diff --git a/pyramid/kernel.py b/pyramid/kernel.py index 7bb72410d3b30f16150cbd885e13f8efc6774207..accd876e5e5f96dbbfbce9de129210da5d9270c5 100644 --- a/pyramid/kernel.py +++ b/pyramid/kernel.py @@ -14,6 +14,11 @@ from jutil import fft __all__ = ['Kernel', 'PHI_0'] PHI_0 = 2067.83 # magnetic flux in T*nm² +H_BAR = 6.626E-34 # Planck constant in J*s +M_E = 9.109E-31 # electron mass in kg +Q_E = 1.602E-19 # electron charge in C +C = 2.998E8 # speed of light in m/s +EPS_0 = 8.8542E-12 # electrical field constant class Kernel(object): @@ -163,3 +168,130 @@ class Kernel(object): print('Grid spacing : {} nm'.format(self.a)) print('Geometry :', self.geometry) print('PRW vector : {} T'.format(self.prw_vec)) + + +class KernelCharge(object): + """Class for calculating kernel matrices for the phase calculation. + + Represents the phase of a single charged pixel, which can be accessed via the corresponding attributes. + During the construction, a few attributes are calculated that are used in + the convolution during phase calculation in the different :class:`~Phasemapper` classes. + + Attributes + ---------- + a : float + The grid spacing in nm. + v_acc : float, optional + The acceleration voltage of the electron microscope in V. The default is 300000. + electrode_vec : tuple of float (N=2) + The norm vector of the counter electrode, (elec_a,elec_b), and the distance to the origin is + the norm of (elec_a,elec_b). + dim_uv : tuple of int (N=2), optional + Dimensions of the 2-dimensional electrostatic charge grid from which the phase should + be calculated. + dim_kern : tuple of int (N=2) + Dimensions of the kernel, which is ``2N-1`` for both axes compared to `dim_uv`. + dim_pad : tuple of int (N=2) + Dimensions of the padded FOV, which is ``2N`` (if FFTW is used) or the next highest power + of 2 (for numpy-FFT). + dim_fft : tuple of int (N=2) + Dimensions of the grid, which is used for the FFT, taking into account that a RFFT should + be used (one axis is halved in comparison to `dim_pad`). + kc : :class:`~numpy.ndarray` (N=3) + The phase contribution of one charged pixel. + kc_fft : :class:`~numpy.ndarray` (N=3) + The real FFT of the phase contribution of one charged pixel. + slice_phase : tuple (N=2) of :class:`slice` + A tuple of :class:`slice` objects to extract the original FOV from the increased one with + size `dim_pad` for the elementary kernel phase. The kernel is shifted, thus the center is + not at (0, 0), which also shifts the slicing compared to `slice_mag`. + slice_mag : tuple (N=2) of :class:`slice` + A tuple of :class:`slice` objects to extract the original FOV from the increased one with + size `dim_pad` for the projected charge distribution. + prw_vec: tuple of 2 int, optional + A two-component vector describing the displacement of the reference wave to include + perturbation of this reference by the object itself (via fringing fields), (y, x). + dtype: numpy dtype, optional + Data type of the kernel. Default is np.float32. + + """ + + _log = logging.getLogger(__name__ + '.KernelCharge') + + def __init__(self, a, dim_uv, electrode_vec, v_acc=300000., prw_vec=None, dtype=np.float32): + self._log.debug('Calling __init__') + # Set basic properties: + self.prw_vec = prw_vec + self.dim_uv = dim_uv # Dimensions of the FOV + self.dim_kern = tuple(2 * np.array(dim_uv) - 1) # Dimensions of the kernel + self.a = a + self.electrode_vec = electrode_vec + self.v_acc = v_acc + lam = H_BAR / np.sqrt(2 * M_E * Q_E * v_acc * (1 + Q_E * v_acc / (2 * M_E * C ** 2))) + c_e = 2 * np.pi * Q_E / lam * (Q_E * v_acc + M_E * C ** 2) / ( + Q_E * v_acc * (Q_E * v_acc + 2 * M_E * C ** 2)) + # Set up FFT: + if fft.HAVE_FFTW: + self.dim_pad = tuple(2 * np.array(dim_uv)) # is at least even (not nec. power of 2) + else: + self.dim_pad = tuple(2 ** np.ceil(np.log2(2 * np.array(dim_uv))).astype(int)) # pow(2) + self.dim_fft = (self.dim_pad[0], self.dim_pad[1] // 2 + 1) # last axis is real + self.slice_phase = (slice(dim_uv[0] - 1, self.dim_kern[0]), # Shift because kernel center + slice(dim_uv[1] - 1, self.dim_kern[1])) # is not at (0, 0)! + self.slice_mag = (slice(0, dim_uv[0]), # Magnetization is padded on the far end! + slice(0, dim_uv[1])) # (Phase cutout is shifted as listed above) + # Calculate kernel (single pixel phase): + coeff = c_e * Q_E / (4 * np.pi * EPS_0) # Minus is gone because of negative z-direction + v_dim, u_dim = dim_uv + u = np.linspace(-(u_dim - 1), u_dim - 1, num=2 * u_dim - 1) + v = np.linspace(-(v_dim - 1), v_dim - 1, num=2 * v_dim - 1) + uu, vv = np.meshgrid(u, v) + self.kc = np.empty(self.dim_kern, dtype=dtype) + self.kc[...] = coeff * self._get_elementary_phase(electrode_vec, uu, vv, a) + # Include perturbed reference wave: + if prw_vec is not None: + uu += prw_vec[1] + vv += prw_vec[0] + self.kc[...] -= coeff * self._get_elementary_phase(electrode_vec, uu, vv, a) + # Calculate Fourier transform of kernel: + self.kc_fft = fft.rfftn(self.kc, self.dim_pad) + self._log.debug('Created ' + str(self)) + + def __repr__(self): + self._log.debug('Calling __repr__') + return '%s(a=%r, dim_uv=%r, electrode_vec=%r,prw_vec=%r)' % \ + (self.__class__, self.a, self.dim_uv, self.electrode_vec, self.prw_vec) + + def __str__(self): + self._log.debug('Calling __str__') + return 'Kernel(a=%s, dim_uv=%s, electrode_vec=%s, prw_vec=%s,)' % \ + (self.a, self.dim_uv, self.electrode_vec, self.prw_vec) + + def _get_elementary_phase(self, electrode_vec, n, m, a): + self._log.debug('Calling _get_elementary_phase') + elec_a, elec_b = electrode_vec + n_img = 2 * elec_a + m_img = 2 * elec_b + in_or_out1 = ~ np.logical_and(n == 0, m == 0) + in_or_out2 = ~ np.logical_and((n - n_img) == 0, (m - m_img) == 0) + return (1. / np.sqrt(a ** 2 * (n ** 2 + m ** 2 + 1E-30))) * in_or_out1 - \ + (1. / np.sqrt(a ** 2 * ((n - n_img) ** 2 + (m - m_img) ** 2 + 1E-30))) * in_or_out2 + + def print_info(self): + """Print information about the kernel. + + Returns + ------- + None + + """ + self._log.debug('Calling log_info') + print('Shape of the FOV :', self.dim_uv) + print('Shape of the Kernel :', self.dim_kern) + print('Zero-padded shape :', self.dim_pad) + print('Shape of the FFT :', self.dim_fft) + print('Slice for the phase :', self.slice_phase) + print('Slice for the magn. :', self.slice_mag) + print('Grid spacing : {} nm'.format(self.a)) + print('PRW vector : {} T'.format(self.prw_vec)) + print('Electrode vector : {} T'.format(self.electrode_vec)) diff --git a/pyramid/kernelcharge.py b/pyramid/kernelcharge.py deleted file mode 100644 index 5625c2dfe489d036f35352672d52f05ae2311c7a..0000000000000000000000000000000000000000 --- a/pyramid/kernelcharge.py +++ /dev/null @@ -1,151 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright 2014 by Forschungszentrum Juelich GmbH -# Author: J. Caron -# -"""This module provides the :class:`~.KernelCharge` class, representing the phase contribution of one -single Charged pixel.""" - -import logging - -import numpy as np - -from jutil import fft - -__all__ = ['Kernel', 'PHI_0'] # TODO rewrite! - -# PHI_0 = 2067.83 # magnetic flux in T*nm² -H_BAR = 6.626E-34 # Planck constant in J*s -M_E = 9.109E-31 # electron mass in kg -Q_E = 1.602E-19 # electron charge in C -C = 2.998E8 # speed of light in m/s -EPS_0 = 8.8542E-12 # electrical field constant - - -class KernelCharge(object): - """Class for calculating kernel matrices for the phase calculation. - - Represents the phase of a single charged pixel, which can be accessed via the corresponding attributes. - During the construction, a few attributes are calculated that are used in - the convolution during phase calculation in the different :class:`~Phasemapper` classes. - - Attributes - ---------- - a : float - The grid spacing in nm. - v_acc : float, optional - The acceleration voltage of the electron microscope in V. The default is 300000. - electrode_vec : tuple of float (N=2) - The norm vector of the counter electrode, (elec_a,elec_b), and the distance to the origin is - the norm of (elec_a,elec_b). - dim_uv : tuple of int (N=2), optional - Dimensions of the 2-dimensional projected magnetization grid from which the phase should - be calculated. - dim_kern : tuple of int (N=2) - Dimensions of the kernel, which is ``2N-1`` for both axes compared to `dim_uv`. - dim_pad : tuple of int (N=2) - Dimensions of the padded FOV, which is ``2N`` (if FFTW is used) or the next highest power - of 2 (for numpy-FFT). - dim_fft : tuple of int (N=2) - Dimensions of the grid, which is used for the FFT, taking into account that a RFFT should - be used (one axis is halved in comparison to `dim_pad`). - kc : :class:`~numpy.ndarray` (N=3) - The phase contribution of one charged pixel. - kc_fft : :class:`~numpy.ndarray` (N=3) - The real FFT of the phase contribution of one charged pixel. - slice_phase : tuple (N=2) of :class:`slice` - A tuple of :class:`slice` objects to extract the original FOV from the increased one with - size `dim_pad` for the elementary kernel phase. The kernel is shifted, thus the center is - not at (0, 0), which also shifts the slicing compared to `slice_mag`. - slice_mag : tuple (N=2) of :class:`slice` - A tuple of :class:`slice` objects to extract the original FOV from the increased one with - size `dim_pad` for the projected magnetization distribution. - prw_vec: tuple of 2 int, optional - A two-component vector describing the displacement of the reference wave to include - perturbation of this reference by the object itself (via fringing fields), (y, x). - dtype: numpy dtype, optional - Data type of the kernel. Default is np.float32. - - """ - - _log = logging.getLogger(__name__ + '.Kernel') # TODO 'KernelCharge' - - def __init__(self, a, dim_uv, electrode_vec, v_acc=300000., prw_vec=None, dtype=np.float32): - self._log.debug('Calling __init__') - # Set basic properties: - self.prw_vec = prw_vec - self.dim_uv = dim_uv # Dimensions of the FOV - self.dim_kern = tuple(2 * np.array(dim_uv) - 1) # Dimensions of the kernel - self.a = a - self.electrode_vec = electrode_vec - self.v_acc = v_acc - lam = H_BAR / np.sqrt(2 * M_E * Q_E * v_acc * (1 + Q_E * v_acc / (2 * M_E * C ** 2))) - c_e = 2 * np.pi * Q_E / lam * (Q_E * v_acc + M_E * C ** 2) / ( - Q_E * v_acc * (Q_E * v_acc + 2 * M_E * C ** 2)) - - # Set up FFT: - if fft.HAVE_FFTW: - - self.dim_pad = tuple(2 * np.array(dim_uv)) # is at least even (not nec. power of 2) - else: - self.dim_pad = tuple(2 ** np.ceil(np.log2(2 * np.array(dim_uv))).astype(int)) # pow(2) - - self.dim_fft = (self.dim_pad[0], self.dim_pad[1] // 2 + 1) # last axis is real - self.slice_phase = (slice(dim_uv[0] - 1, self.dim_kern[0]), # Shift because kernel center - slice(dim_uv[1] - 1, self.dim_kern[1])) # is not at (0, 0)! - self.slice_mag = (slice(0, dim_uv[0]), # Magnetization is padded on the far end! - slice(0, dim_uv[1])) # (Phase cutout is shifted as listed above) - # Calculate kernel (single pixel phase): - coeff = c_e * Q_E / (4 * np.pi * EPS_0) # Minus is gone because of negative z-direction - v_dim, u_dim = dim_uv - u = np.linspace(-(u_dim - 1), u_dim - 1, num=2 * u_dim - 1) - v = np.linspace(-(v_dim - 1), v_dim - 1, num=2 * v_dim - 1) - uu, vv = np.meshgrid(u, v) - self.kc = np.empty(self.dim_kern, dtype=dtype) - self.kc[...] = coeff * self._get_elementary_phase(electrode_vec, uu, vv, a) - # Include perturbed reference wave: - if prw_vec is not None: - uu += prw_vec[1] - vv += prw_vec[0] - self.kc[...] -= coeff * self._get_elementary_phase(electrode_vec, uu, vv, a) - # Calculate Fourier transform of kernel: - self.kc_fft = fft.rfftn(self.kc, self.dim_pad) - self._log.debug('Created ' + str(self)) - - def __repr__(self): - self._log.debug('Calling __repr__') - return '%s(a=%r, dim_uv=%r, electrode_vec=%r,prw_vec=%r)' % \ - (self.__class__, self.a, self.dim_uv, self.electrode_vec, self.prw_vec) - - def __str__(self): - self._log.debug('Calling __str__') - return 'Kernel(a=%s, dim_uv=%s, electrode_vec=%s, prw_vec=%s,)' % \ - (self.a, self.dim_uv, self.electrode_vec, self.prw_vec) - - def _get_elementary_phase(self, electrode_vec, n, m, a): - self._log.debug('Calling _get_elementary_phase') - elec_a, elec_b = electrode_vec - n_img = 2 * elec_a - m_img = 2 * elec_b - in_or_out1 = ~ np.logical_and(n == 0, m == 0) - in_or_out2 = ~ np.logical_and((n - n_img) == 0, (m - m_img) == 0) - return (1. / np.sqrt(a ** 2 * (n ** 2 + m ** 2 + 1E-30))) * in_or_out1 - \ - (1. / np.sqrt(a ** 2 * ((n - n_img) ** 2 + (m - m_img) ** 2 + 1E-30))) * in_or_out2 - - def print_info(self): - """Print information about the kernel. - - Returns - ------- - None - - """ - self._log.debug('Calling log_info') - print('Shape of the FOV :', self.dim_uv) - print('Shape of the Kernel :', self.dim_kern) - print('Zero-padded shape :', self.dim_pad) - print('Shape of the FFT :', self.dim_fft) - print('Slice for the phase :', self.slice_phase) - print('Slice for the magn. :', self.slice_mag) - print('Grid spacing : {} nm'.format(self.a)) - print('PRW vector : {} T'.format(self.prw_vec)) - print('Electrode vector : {} T'.format(self.electrode_vec)) diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py index 27d1a6751818cdeddc7e953e174613653bc6b04f..a8852a8e8b60848644aaabd5b6974c0bf39fbac8 100644 --- a/pyramid/phasemapper.py +++ b/pyramid/phasemapper.py @@ -434,6 +434,7 @@ class PhaseMapperCharge(PhaseMapper): Size of the input space. """ + _log = logging.getLogger(__name__ + '.PhaseMapperCharge') def __init__(self, a, dim_uv, electrode_vec, v_acc=300000): self._log.debug('Calling __init__') diff --git a/pyramid/tests/.DS_Store b/pyramid/tests/.DS_Store deleted file mode 100644 index 5008ddfcf53c02e82d7eee2e57c38e5672ef89f6..0000000000000000000000000000000000000000 Binary files a/pyramid/tests/.DS_Store and /dev/null differ diff --git a/pyramid/tests/test_kernelcharge.py b/pyramid/tests/test_kernelcharge.py deleted file mode 100644 index 7b36bb96456f602ebcb4ef020b988029a7221df8..0000000000000000000000000000000000000000 --- a/pyramid/tests/test_kernelcharge.py +++ /dev/null @@ -1,27 +0,0 @@ -# -*- coding: utf-8 -*- -"""Testcase for the magdata module.""" - -import os -import unittest - -import numpy as np -from numpy.testing import assert_allclose - -from pyramid.kernelcharge import KernelCharge - - -class TestCaseKernelCharge(unittest.TestCase): - def setUp(self): - self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_kernelcharge') - self.kernel = KernelCharge(1., dim_uv=(8, 8), electrode_vec=(3, 3)) - - def tearDown(self): - self.path = None - self.kernel = None - - def test_kernelcharge(self): - ref_kc = np.load(os.path.join(self.path, 'ref_kc.npy')) - ref_kc_fft = np.load(os.path.join(self.path, 'ref_kc_fft.npy')) - assert_allclose(self.kernel.kc, ref_kc, err_msg='Unexpected behavior in kc') - assert_allclose(self.kernel.kc_fft, ref_kc_fft, atol=1E-7, - err_msg='Unexpected behavior in kc_fft') \ No newline at end of file diff --git a/tests/test_kernel.py b/tests/test_kernel.py index b4273e5651d9198f7ba250a8ea153bc947a2fb18..4fdd2daca705bed85e130c8a153f455a59f492ec 100644 --- a/tests/test_kernel.py +++ b/tests/test_kernel.py @@ -7,7 +7,7 @@ import unittest import numpy as np from numpy.testing import assert_allclose -from pyramid.kernel import Kernel +from pyramid.kernel import Kernel, KernelCharge class TestCaseKernel(unittest.TestCase): @@ -30,3 +30,20 @@ class TestCaseKernel(unittest.TestCase): err_msg='Unexpected behavior in u_fft') assert_allclose(self.kernel.v_fft, ref_v_fft, atol=1E-7, err_msg='Unexpected behavior in v_fft') + + +class TestCaseKernelCharge(unittest.TestCase): + def setUp(self): + self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_kernel') + self.kernel = KernelCharge(1., dim_uv=(8, 8), electrode_vec=(3, 3)) + + def tearDown(self): + self.path = None + self.kernel = None + + def test_kernelcharge(self): + ref_kc = np.load(os.path.join(self.path, 'ref_kc.npy')) + ref_kc_fft = np.load(os.path.join(self.path, 'ref_kc_fft.npy')) + assert_allclose(self.kernel.kc, ref_kc, err_msg='Unexpected behavior in kc') + assert_allclose(self.kernel.kc_fft, ref_kc_fft, atol=1E-7, + err_msg='Unexpected behavior in kc_fft') diff --git a/pyramid/tests/test_kernelcharge/ref_kc.npy b/tests/test_kernel/ref_kc.npy similarity index 100% rename from pyramid/tests/test_kernelcharge/ref_kc.npy rename to tests/test_kernel/ref_kc.npy diff --git a/pyramid/tests/test_kernelcharge/ref_kc_fft.npy b/tests/test_kernel/ref_kc_fft.npy similarity index 100% rename from pyramid/tests/test_kernelcharge/ref_kc_fft.npy rename to tests/test_kernel/ref_kc_fft.npy