KernelCharge has been updated and other some revision has beed done in phasemappercharge.

pyramid/kernel.py

@@ -82,6 +82,7 @@ class Kernel(object):
         self.geometry = geometry
         # 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)

pyramid/kernelcharge.py

+# -*- 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))

pyramid/phasemapper.py

@@ -410,8 +410,30 @@ class PhaseMapperMIP(PhaseMapper):
 
 
 class PhaseMapperCharge(PhaseMapper):
+    """Class representing a phase mapping strategy for the electrostatic charge contribution.
 
-    """"""  # TODO: Write Docstring!
+    The :class:`~.PhaseMapperCharge` class represents a phase mapping strategy for the electrostatic
+    charge contribution to the electron phase shift which results e.g. from the charges in
+    certain samples and which is sensitive to properties of the electron microscope. It directly
+    takes :class:`~.ScalarData` objects and returns :class:`~.PhaseMap` objects.
+
+    Attributes
+    ----------
+    a : float
+        The grid spacing in nm.
+    dim_uv : tuple of int (N=2)
+        Dimensions of the 2-dimensional projected magnetization grid for the kernel setup.
+    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).
+    m: int
+        Size of the image space.
+    n: int
+        Size of the input space.
+
+    """
 
     def __init__(self, a, dim_uv, electrode_vec, v_acc=300000):
         self._log.debug('Calling __init__')
@@ -440,8 +462,6 @@ class PhaseMapperCharge(PhaseMapper):
         """ This is to calculate the phase from many electric dipoles. The model used to avoid singularity is
         the homogeneously-distributed charged metallic sphere.
         The elec_data includes the amount of charge in every grid, unit:electron.
-        The electrode_vec is the  normal vector of the electrode, (elec_a,elec_b), and the distance to the origin is
-        the norm of (elec_a,elec_b).
         R_sam is the sampling rate,pixel/nm."""
 
         R_sam = 1 / self.a
@@ -449,9 +469,7 @@ class PhaseMapperCharge(PhaseMapper):
         field = elec_data.field[0, ...]
         elec_a, elec_b = self.electrode_vec  # electrode vector (orthogonal to the electrode from origin)
         elec_n = elec_a ** 2 + elec_b ** 2
-
         dim_v, dim_u = field.shape
-
         u_cor = range(0, dim_u, 1)
         v_cor = range(0, dim_v, 1)
         v, u = np.meshgrid(v_cor, u_cor)
@@ -465,51 +483,29 @@ class PhaseMapperCharge(PhaseMapper):
         vm = (elec_a ** 2 - elec_b ** 2) / elec_n * vq - 2 * elec_a * elec_b / elec_n * uq + 2 * elec_b
         # Calculate phase contribution for each charge:
         phase = np.zeros((dim_v, dim_u))
-
         for i in range(len(q)):
-
             # The projected distance from the charges or image charges
-
             r1 = np.sqrt((u - uq[i]) ** 2 + (v - vq[i]) ** 2)
-
             r2 = np.sqrt((u - um[i]) ** 2 + (v - vm[i]) ** 2)
-
             # The square height when  the path come across the sphere
-
-            z1 = R ** 2 - r1 ** 2
-
-            z2 = R ** 2 - r2 ** 2
-
+            h1 = R ** 2 - r1 ** 2
+            h2 = R ** 2 - r2 ** 2
             # Phase calculation in 3 different cases
-
             # case 1 totally out of the sphere
-
-            case1 = ((z1 < 0) & (z2 < 0))
-
+            case1 = ((h1 < 0) & (h2 < 0))
             phase[case1] += - q[i] * self.coeff * np.log((r1[case1] ** 2) / (r2[case1] ** 2))
-
             # case 2: inside the charge sphere
-
-            case2 = ((z1 >= 0) & (z2 <= 0))
-
+            case2 = ((h1 >= 0) & (h2 <= 0))
             # The height when the path come across the charge sphere
-
-            z3 = np.sqrt(z1)
-
-            phase[case2] += q[i] * self.coeff * (- np.log((z3[case2] + R) ** 2 / r2[case2] ** 2) +
-                             (2 * z3[case2] / R + 2 * z3[case2] ** 3 / 3 / R ** 3))
-
+            h3 = np.sqrt(h1)
+            phase[case2] += q[i] * self.coeff * (- np.log((h3[case2] + R) ** 2 / r2[case2] ** 2) +
+                                                 (2 * h3[case2] / R + 2 * h3[case2] ** 3 / 3 / R ** 3))
             # case 3 : inside the image charge sphere
-
             case3 = np.logical_not(case1 | case2)
-
             # The height whe the path comes across the image charge sphere
-
-            z4 = np.sqrt(z2)
-
-            phase[case3] += q[i] * self.coeff * (np.log((z4[case3] + R) ** 2 / r1[case3] ** 2) -
-                             (2 * z4[case3] / R + 2 * z4[case3] ** 3 / 3 / R ** 3))
-
+            h4 = np.sqrt(h2)
+            phase[case3] += q[i] * self.coeff * (np.log((h4[case3] + R) ** 2 / r1[case3] ** 2) -
+                                                 (2 * h4[case3] / R + 2 * h4[case3] ** 3 / 3 / R ** 3))
         return PhaseMap(self.a, phase)
 
     def jac_dot(self, vector):

pyramid/tests/test_kernel.py

@@ -20,9 +20,9 @@ class TestCaseKernel(unittest.TestCase):
         self.kernel = None
 
     def test_kernel(self):
-        ref_u = np.load(os.path.join(self.path, 'ref_u.npy'))
+        ref_u = np.load(os.path.join(self.path, 'ref_kc.npy'))
         ref_v = np.load(os.path.join(self.path, 'ref_v.npy'))
-        ref_u_fft = np.load(os.path.join(self.path, 'ref_u_fft.npy'))
+        ref_u_fft = np.load(os.path.join(self.path, 'ref_kc_fft.npy'))
         ref_v_fft = np.load(os.path.join(self.path, 'ref_v_fft.npy'))
         assert_allclose(self.kernel.u, ref_u, err_msg='Unexpected behavior in u')
         assert_allclose(self.kernel.v, ref_v, err_msg='Unexpected behavior in v')

pyramid/tests/test_kernelcharge.py

+# -*- 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