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

PhaseMapperCharge has been updated.

parent 8bf61ba6
No related branches found
No related tags found
No related merge requests found
Showing
with 153 additions and 33 deletions
Metadata-Version: 1.1
Name: pyramid
Version: 0.1.0.dev0
Summary: PYthon based Reconstruction Algorithm for MagnetIc Distributions
Home-page: UNKNOWN
Author: Jan Caron
Author-email: j.caron@fz-juelich.de
License: UNKNOWN
Description: long description (TODO!)
Platform: UNKNOWN
Requires: numpy
Requires: scipy
Requires: matplotlib
Requires: Pillow
Requires: mayavi
Requires: pyfftw
Requires: hyperspy
Requires: nose
Requires: jutil
setup.py
pyramid/__init__.py
pyramid/analytic.py
pyramid/colors.py
pyramid/costfunction.py
pyramid/dataset.py
pyramid/diagnostics.py
pyramid/fft.py
pyramid/fieldconverter.py
pyramid/fielddata.py
pyramid/forwardmodel.py
pyramid/kernel.py
pyramid/phasemap.py
pyramid/phasemapper.py
pyramid/plottools.py
pyramid/projector.py
pyramid/quaternion.py
pyramid/ramp.py
pyramid/reconstruction.py
pyramid/regularisator.py
pyramid/version.py
pyramid.egg-info/PKG-INFO
pyramid.egg-info/SOURCES.txt
pyramid.egg-info/dependency_links.txt
pyramid.egg-info/top_level.txt
pyramid/file_io/__init__.py
pyramid/file_io/io_dataset.py
pyramid/file_io/io_phasemap.py
pyramid/file_io/io_projector.py
pyramid/file_io/io_scalardata.py
pyramid/file_io/io_vectordata.py
pyramid/magcreator/__init__.py
pyramid/magcreator/examples.py
pyramid/magcreator/magcreator.py
pyramid/magcreator/shapes.py
pyramid/utils/__init__.py
pyramid/utils/mag_slicer.py
pyramid/utils/phasemap_creator.py
pyramid/utils/pm.py
pyramid/utils/reconstruction_2d_from_phasemap.py
pyramid/utils/reconstruction_3d_from_magdata.py
\ No newline at end of file
pyramid
File added
File added
File added
File added
File added
File added
File added
File added
...@@ -411,13 +411,13 @@ class PhaseMapperMIP(PhaseMapper): ...@@ -411,13 +411,13 @@ class PhaseMapperMIP(PhaseMapper):
class PhaseMapperCharge(PhaseMapper): class PhaseMapperCharge(PhaseMapper):
"""""" #This part will be updated. """"""
def __init__(self, a, dim_uv, biprism_vec, v_acc=300000): def __init__(self, a, dim_uv, electrode_vec, v_acc=300000):
self._log.debug('Calling __init__') self._log.debug('Calling __init__')
self.a = a self.a = a
self.dim_uv = dim_uv self.dim_uv = dim_uv
self.biprism_vec = biprism_vec self.electrode_vec = electrode_vec
self.v_acc = v_acc 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))) 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) / (
...@@ -437,43 +437,76 @@ class PhaseMapperCharge(PhaseMapper): ...@@ -437,43 +437,76 @@ class PhaseMapperCharge(PhaseMapper):
(self.a, self.dim_uv, self.v_acc) (self.a, self.dim_uv, self.v_acc)
def __call__(self, elec_data): def __call__(self, elec_data):
""" phase_dipoles() is to caculate the phase from many electric dipoles. """ phase_dipoles() is to calculate the phase from many electric dipoles. The model used to avoid singularity is
field include the amount of charge in every grid, unit:electron. the homogeneously-distributed charged metallic sphere.
The normal vector of the electrode is (a,b), and the distance to the origin is the norm The elec_data includes the amount of charge in every grid, unit:electron.
of (a,b). R is the sampling rate,pixel/nm.""" The electrode_vec is the normal vector of the electrode, (elec_a,elec_b), and the distance to the origin is
R = 1 / self.a the norm of (elec_a,elec_b).
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, ...] field = elec_data.field[0, ...]
bu, bv = self.biprism_vec # biprism vector (orthogonal to biprism from origin) elec_a, elec_b = self.electrode_vec # electrode vector (orthogonal to the electrode from origin)
bn = bu ** 2 + bv ** 2 # norm of the biprism vector elec_n = elec_a ** 2 + elec_b ** 2
dim_v, dim_u = field.shape dim_v, dim_u = field.shape
v, u = np.meshgrid(dim_v, dim_u)
# Find list of charge locations and charge values: # Find list of charge locations and charge values:
vq, uq = np.nonzero(field) vq, uq = np.nonzero(field)
q = field[vq, uq] q = field[vq, uq]
# Calculate locations of image charges: # Calculate locations of image charges:
vm = (bu ** 2 - bv ** 2) / bn * vq - 2 * bu * bv / bn * uq + 2 * bv # 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 = (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: # Calculate phase contribution for each charge:
phase = np.zeros((dim_v, dim_u)) phase = np.zeros((dim_v, dim_u))
for i in range(len(q)): for i in range(len(q)):
for u in range(dim_u):
for v in range(dim_v): # The projected distance from the charges or image charges
rq = np.sqrt((u - uq[i]) ** 2 + (v - vq[i]) ** 2) # charge distance
rm = np.sqrt((u - um[i]) ** 2 + (v - vm[i]) ** 2) # mirror distance r1 = np.sqrt((u - uq[i]) ** 2 + (v - vq[i]) ** 2)
# Distance:
z1 = (R / 2) ** 2 - rq ** 2 r2 = np.sqrt((u - um[i]) ** 2 + (v - vm[i]) ** 2)
z2 = (R / 2) ** 2 - rm ** 2
if z1 < 0 and z2 < 0: # The square height when the path come across the sphere
phase[v, u] += -q[i] * self.coeff * np.log((rq ** 2) / (rm ** 2))
elif z1 >= 0 >= z2: z1 = R ** 2 - r1 ** 2
z3 = np.sqrt(z1)
phase[v, u] += (-q[i] * self.coeff * z2 = R ** 2 - r2 ** 2
np.log(((z3 + R / 2) ** 2) / (rm ** 2))
+ 2 * q[i] * Q_E * z3 / 2 / np.pi / EPS_0 / R) # Phase calculation in 3 different cases
else:
z4 = np.sqrt(z2) # case 1 totally out of the sphere
phase[v, u] += (-q[i] * self.coeff *
np.log((rq ** 2) / ((z4 + R / 2) ** 2)) case1 = ((z1 < 0) & (z2 < 0))
- 2 * q[i] * Q_E * z4 / 2 / np.pi / EPS_0 / R)
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))
# 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))
# 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))
return PhaseMap(self.a, phase) return PhaseMap(self.a, phase)
def jac_dot(self, vector): def jac_dot(self, vector):
......
...@@ -8,7 +8,7 @@ import numpy as np ...@@ -8,7 +8,7 @@ import numpy as np
from numpy.testing import assert_allclose from numpy.testing import assert_allclose
from pyramid.kernel import Kernel from pyramid.kernel import Kernel
from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC, PhaseMapperMIP from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC, PhaseMapperMIP, PhaseMapperCharge
from pyramid import load_phasemap, load_vectordata, load_scalardata from pyramid import load_phasemap, load_vectordata, load_scalardata
...@@ -177,6 +177,31 @@ class TestCasePhaseMapperMIP(unittest.TestCase): ...@@ -177,6 +177,31 @@ class TestCasePhaseMapperMIP(unittest.TestCase):
self.assertRaises(NotImplementedError, self.mapper.jac_T_dot, None) self.assertRaises(NotImplementedError, self.mapper.jac_T_dot, None)
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:])
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'))
phasemap = self.mapper(self.Charge_proj)
assert_allclose(phasemap.phase, Charge_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__()!')
def test_jac_dot(self):
self.assertRaises(NotImplementedError, self.mapper.jac_dot, None)
def test_jac_T_dot(self):
self.assertRaises(NotImplementedError, self.mapper.jac_T_dot, None)
if __name__ == '__main__': if __name__ == '__main__':
import nose import nose
nose.run(defaultTest=__name__) nose.run(defaultTest=__name__)
File added
File added
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
""""This file is generated automatically by the Pyramid `setup.py`""" """"This file is generated automatically by the Pyramid `setup.py`"""
version = "0.1.0-dev" version = "0.1.0-dev"
hg_revision = "???" # TODO: Now uses git! Maybe delete alltogether? See setup.py! hg_revision = "???"
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment