Skip to content
Snippets Groups Projects
Commit 04691a22 authored by Fengshan Zheng's avatar Fengshan Zheng
Browse files

KernelCharge updated but not finished.

parent 08b4e196
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
pyramid/kernel.py
+ 35
− 11
View file @ 04691a22
......@@ -182,7 +182,7 @@ class KernelCharge(object):
a : float
The grid spacing in nm.
v_acc : float, optional
The acceleration voltage of the electron microscope in V. The default is 300000.
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).
......@@ -228,7 +228,7 @@ class KernelCharge(object):
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) / (
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:
......@@ -241,18 +241,18 @@ class KernelCharge(object):
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
self.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)
self.kc[...] = self.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)
self.kc[...] -= self.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))
......@@ -269,13 +269,37 @@ class KernelCharge(object):
def _get_elementary_phase(self, electrode_vec, n, m, a):
self._log.debug('Calling _get_elementary_phase')
phase = np.zeros((self.dim_uv[1], self.dim_uv[0]))
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
u_img = 2 * elec_a
v_img = 2 * elec_b
# in_or_out1 = ~ np.logical_and(n == 0, m == 0)
# in_or_out2 = ~ np.logical_and((n - u_img) == 0, (m - v_img) == 0)
# The projected distance from the charges or image charges
r1 = np.sqrt(n ** 2 + m ** 2)
r2 = np.sqrt((n - u_img) ** 2 + (n - v_img) ** 2)
# The square height when the path come across the sphere
R = a / 2 # the radius of the pixel
h1 = R ** 2 - r1 ** 2
h2 = R ** 2 - r2 ** 2
# Phase calculation in 3 different cases
# case 1 totally out of the sphere
case1 = ((h1 < 0) & (h2 < 0))
phase[case1] += - self.coeff * np.log((r1[case1] ** 2) / (r2[case1] ** 2))
# case 2: inside the charge sphere
case2 = ((h1 >= 0) & (h2 <= 0))
# The height when the path come across the charge sphere
h3 = np.sqrt(h1)
phase[case2] += 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
h4 = np.sqrt(h2)
phase[case3] += self.coeff * (np.log((h4[case3] + R) ** 2 / r1[case3] ** 2) -
(2 * h4[case3] / R + 2 * h4[case3] ** 3 / 3 / R ** 3))
return phase
def print_info(self):
"""Print information about the kernel.
......
pyramid/phasemapper.py
+ 26
− 72
View file @ 04691a22
......@@ -419,15 +419,8 @@ class PhaseMapperCharge(PhaseMapper):
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).
kernelcharge : :class:`~pyramid.KernelCharge`
Convolution kernel, representing the phase contribution of one single charge pixel.
m: int
Size of the image space.
n: int
......@@ -436,78 +429,39 @@ class PhaseMapperCharge(PhaseMapper):
"""
_log = logging.getLogger(__name__ + '.PhaseMapperCharge')
def __init__(self, a, dim_uv, electrode_vec, v_acc=300000):
def __init__(self, kernelcharge):
self._log.debug('Calling __init__')
self.a = a
self.dim_uv = dim_uv
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))
self.coeff = C_e * Q_E / (4 * np.pi * EPS_0)
self.kernelcharge = kernelcharge
self.m = np.prod(kernelcharge.dim_uv)
self.n = 2 * self.m
self.c = np.zeros(kernelcharge.dim_pad, dtype=kernelcharge.u.dtype)
self.phase_adj = np.zeros(kernelcharge.dim_pad, dtype=kernelcharge.u.dtype)
self._log.debug('Created ' + str(self))
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, v_acc=%r)' % \
(self.__class__, self.a, self.dim_uv, self.v_acc)
return '%s(kernelcharge=%r)' % (self.__class__, self.kernelcharge)
def __str__(self):
self._log.debug('Calling __str__')
return 'PhaseMapperCharge(a=%s, dim_uv=%s, v_acc=%s)' % \
(self.a, self.dim_uv, self.v_acc)
return 'PhaseMapperCharge(kernelcharge=%s)' % self.kernelcharge
def __call__(self, elec_data):
""" 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.
R_sam is the sampling rate,pixel/nm."""
R_sam = 1 / self.a
R = R_sam / 2 # The radius of the charged sphere
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)
# Find list of charge locations and charge values:
vq, uq = np.nonzero(field)
q = field[vq, uq]
# Calculate locations of image charges:
# vm = (bu ** 2 - bv ** 2) / bn * vq - 2 * bu * bv / bn * uq + 2 * bv
# um = (bv ** 2 - bu ** 2) / bn * uq - 2 * bu * bv / bn * vq + 2 * bu
um = (elec_b ** 2 - elec_a ** 2) / elec_n * uq - 2 * elec_a * elec_b / elec_n * vq + 2 * elec_a
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
h1 = R ** 2 - r1 ** 2
h2 = R ** 2 - r2 ** 2
# Phase calculation in 3 different cases
# case 1 totally out of the sphere
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 = ((h1 >= 0) & (h2 <= 0))
# The height when the path come across the charge sphere
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
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 __call__(self, elecdata):
assert isinstance(elecdata, ScalarData), 'Only ScalarData objects can be mapped!'
assert elecdata.a == self.kernelcharge.a, 'Grid spacing has to match!'
assert elecdata.dim[0] == 1, 'Charge distribution must be 2-dimensional!'
assert elecdata.dim[1:3] == self.kernelcharge.dim_uv, 'Dimensions do not match!'
# Process input parameters:
self.c[self.kernelcharge.slice_mag] = elecdata.field[0, 0, ...]
return PhaseMap(elecdata.a, self._convolve())
def _convolve(self):
# Fourier transform the projected magnetisation:
self.c_fft = fft.rfftn(self.c)
# Convolve the magnetization with the kernel in Fourier space:
self.phase_fft = self.c_fft * self.kernelcharge.kc_fft
# Return the result:
return fft.irfftn(self.phase_fft)[self.kernelcharge.slice_phase]
def jac_dot(self, vector):
"""Calculate the product of the Jacobi matrix with a given `vector`.
......
tests/test_phasemapper.py
+ 23
− 11
View file @ 04691a22
......@@ -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
from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC, PhaseMapperMIP, PhaseMapperCharge
from pyramid import load_phasemap, load_vectordata, load_scalardata
......@@ -181,23 +181,35 @@ class TestCasePhaseMapperCharge(unittest.TestCase):
def setUp(self):
self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_phasemapper')
self.charge_proj = load_scalardata(os.path.join(self.path, 'charge_proj.hdf5'))
self.mapper = PhaseMapperCharge(self.charge_proj.a, self.charge_proj.dim[1:],
electrode_vec=[8, 8], v_acc=300000)
self.mapper = PhaseMapperRDFC(KernelCharge(self.charge_proj.a, self.charge_proj.dim[1:]))
def tearDown(self):
self.path = None
self.charge_proj = None
self.mapper = None
def test_call(self):
charge_phase_ref = load_phasemap(os.path.join(self.path, 'charge_phase_ref.hdf5'))
def test_PhaseMapperCharge_call(self):
phase_ref = load_phasemap(os.path.join(self.path, 'phasemap.hdf5'))
phasemap = self.mapper(self.charge_proj)
assert_allclose(phasemap.phase, charge_phase_ref.phase, atol=1E-7,
assert_allclose(phasemap.phase, phase_ref.phase, atol=1E-7,
err_msg='Unexpected behavior in __call__()!')
assert_allclose(phasemap.a, charge_phase_ref.a, err_msg='Unexpected behavior in __call__()!')
assert_allclose(phasemap.a, phase_ref.a, err_msg='Unexpected behavior in __call__()!')
def test_jac_dot(self):
self.assertRaises(NotImplementedError, self.mapper.jac_dot, None)
def test_PhaseMapperCharge_jac_dot(self):
phase = self.mapper(self.charge_proj).phase
charge_proj_scalar = self.charge_proj.field[:2, ...].ravel()
phase_jac = self.mapper.jac_dot(charge_proj_scalar).reshape(self.mapper.kernelcharge.dim_uv)
assert_allclose(phase, phase_jac, atol=1E-7,
err_msg='Inconsistency between __call__() and jac_dot()!')
n = self.mapper.n
jac = np.array([self.mapper.jac_dot(np.eye(n)[:, i]) for i in range(n)]).T
jac_ref = np.load(os.path.join(self.path, 'jac.npy'))
assert_allclose(jac, jac_ref, atol=1E-7,
err_msg='Unexpected behaviour in the the jacobi matrix!')
def test_jac_T_dot(self):
self.assertRaises(NotImplementedError, self.mapper.jac_T_dot, None)
def test_PhaseMapperCharge_jac_T_dot(self):
m = self.mapper.m
jac_T = np.array([self.mapper.jac_T_dot(np.eye(m)[:, i]) for i in range(m)]).T
jac_T_ref = np.load(os.path.join(self.path, 'jac.npy')).T
assert_allclose(jac_T, jac_T_ref, atol=1E-7,
err_msg='Unexpected behaviour in the the transposed jacobi matrix!')
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment