Forked from
empyre / empyre
415 commits behind the upstream repository.
-
Jan Caron authorednumcore: deleted c1, c2, c3, name changes in phase_mag_real analytic: name change beta to phi phasemapper: changes to include numcore test_compliance: also checkes .pyx files now setup: now works properly (hopefully on other systems, too) compare_discs: added script to compare discs
Jan Caron authorednumcore: deleted c1, c2, c3, name changes in phase_mag_real analytic: name change beta to phi phasemapper: changes to include numcore test_compliance: also checkes .pyx files now setup: now works properly (hopefully on other systems, too) compare_discs: added script to compare discs
phasemapper.py 12.42 KiB
# -*- coding: utf-8 -*-
"""Create and display a phase map from magnetization data."""
import numpy as np
import numcore
import time
# Physical constants
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
def phase_mag_fourier(res, projection, b_0=1, padding=0):
'''Calculate phasemap from magnetization data (Fourier space approach).
Arguments:
res - the resolution of the grid (grid spacing) in nm
projection - projection of a magnetic distribution (created with pyramid.projector)
b_0 - magnetic induction corresponding to a magnetization Mo in T (default: 1)
padding - factor for zero padding, the default is 0 (no padding), for a factor of n the
number of pixels is increase by (1+n)**2, padded zeros are cropped at the end
Returns:
the phasemap as a 2 dimensional array
'''
v_dim, u_dim = np.shape(projection[0])
v_mag, u_mag = projection
# Create zero padded matrices:
u_pad = u_dim/2 * padding
v_pad = v_dim/2 * padding
u_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
v_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
u_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = u_mag
v_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = v_mag
# Fourier transform of the two components:
u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_big), axes=0)
v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_big), axes=0)
# Calculate the Fourier transform of the phase:
nyq = 1 / res # nyquist frequency
f_u = np.linspace(0, nyq/2, u_mag_fft.shape[1])
f_v = np.linspace(-nyq/2, nyq/2, u_mag_fft.shape[0], endpoint=False)
f_uu, f_vv = np.meshgrid(f_u, f_v)
coeff = (1j*b_0) / (2*PHI_0)
phase_fft = coeff * res * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30)
# Transform to real space and revert padding:
phase_big = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
phase = phase_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
return phase
def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
'''Calculate phasemap from magnetization data (real space approach).
Arguments:
res - the resolution of the grid (grid spacing) in nm
projection - projection of a magnetic distribution (created with pyramid.projector)
method - String, describing the method to use for the pixel field ('slab' or 'disc')
b_0 - magnetic induction corresponding to a magnetization Mo in T (default: 1)
jacobi - matrix in which to save the jacobi matrix (default: None, faster computation)
Returns:
the phasemap as a 2 dimensional array
'''
# Function for creating the lookup-tables:
def phi_lookup(method, n, m, res, b_0):
if method == 'slab':
def F_h(n, m):
a = np.log(res**2 * (n**2 + m**2))
b = np.arctan(n / m)
return n*a - 2*n + 2*m*b
return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
- F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
elif method == 'disc':
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
# Process input parameters:
v_dim, u_dim = np.shape(projection[0])
v_mag, u_mag = projection
coeff = -b_0 * res**2 / (2*PHI_0)
# Create lookup-tables for the phase of one pixel:
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)
phi_u = phi_lookup(method, uu, vv, res, b_0)
phi_v = phi_lookup(method, vv, uu, res, b_0)
# Calculation of the phase:
phase = np.zeros((v_dim, u_dim))
threshold = 0
if jacobi is not None: # With Jacobian matrix (slower)
jacobi[:] = 0 # Jacobi matrix --> zeros
############################### TODO: NUMERICAL CORE #####################################
for j in range(v_dim):
for i in range(u_dim):
phase_u = phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
jacobi[:, i+u_dim*j] = phase_u.reshape(-1)
if abs(u_mag[j, i]) > threshold:
phase += u_mag[j, i] * phase_u
phase_v = phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
jacobi[:, u_dim*v_dim+i+u_dim*j] = -phase_v.reshape(-1)
if abs(v_mag[j, i]) > threshold:
phase -= v_mag[j, i] * phase_v
############################### TODO: NUMERICAL CORE #####################################
else: # Without Jacobi matrix (faster)
numcore.phase_mag_real_helper(v_dim, u_dim, phi_u, phi_v, u_mag, v_mag, phase, threshold)
# Return the phase:
return phase
def phase_mag_real2(res, projection, method, b_0=1, jacobi=None):
'''Calculate phasemap from magnetization data (real space approach).
Arguments:
res - the resolution of the grid (grid spacing) in nm
projection - projection of a magnetic distribution (created with pyramid.projector)
method - String, describing the method to use for the pixel field ('slab' or 'disc')
b_0 - magnetic induction corresponding to a magnetization Mo in T (default: 1)
jacobi - matrix in which to save the jacobi matrix (default: None, faster computation)
Returns:
the phasemap as a 2 dimensional array
'''
# Function for creating the lookup-tables:
def phi_lookup(method, n, m, res, b_0):
if method == 'slab':
def F_h(n, m):
a = np.log(res**2 * (n**2 + m**2))
b = np.arctan(n / m)
return n*a - 2*n + 2*m*b
return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
- F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
elif method == 'disc':
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
# Process input parameters:
v_dim, u_dim = np.shape(projection[0])
v_mag, u_mag = projection
coeff = -b_0 * res**2 / (2*PHI_0)
# Create lookup-tables for the phase of one pixel:
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)
phi_u = phi_lookup(method, uu, vv, res, b_0)
phi_v = phi_lookup(method, vv, uu, res, b_0)
# Calculation of the phase:
phase = np.zeros((v_dim, u_dim))
threshold = 0
if jacobi is not None: # With Jacobian matrix (slower)
jacobi[:] = 0 # Jacobi matrix --> zeros
############################### TODO: NUMERICAL CORE #####################################
for j in range(v_dim):
for i in range(u_dim):
phase_u = phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
jacobi[:, i+u_dim*j] = phase_u.reshape(-1)
if abs(u_mag[j, i]) > threshold:
phase += u_mag[j, i] * phase_u
phase_v = phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
jacobi[:, u_dim*v_dim+i+u_dim*j] = -phase_v.reshape(-1)
if abs(v_mag[j, i]) > threshold:
phase -= v_mag[j, i] * phase_v
############################### TODO: NUMERICAL CORE #####################################
else: # Without Jacobi matrix (faster)
# phasecopy = phase.copy()
# start_time = time.time()
# numcore.phase_mag_real_helper_1(v_dim, u_dim, phi_u, phi_v,
# u_mag, v_mag, phasecopy, threshold)
# print 'with numcore : ', time.time() - start_time
# start_time = time.time()
for j in range(v_dim):
for i in range(u_dim):
if abs(u_mag[j, i]) > threshold:
phase += u_mag[j, i] * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
if abs(v_mag[j, i]) > threshold:
phase -= v_mag[j, i] * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
# print 'without numcore: ', time.time() - start_time
# print 'num. difference: ', ((phase - phasecopy) ** 2).sum()
# Return the phase:
return phase
def phase_mag_real_alt(res, projection, method, b_0=1, jacobi=None): # TODO: Modify
'''Calculate phasemap from magnetization data (real space approach).
Arguments:
res - the resolution of the grid (grid spacing) in nm
projection - projection of a magnetic distribution (created with pyramid.projector)
method - String, describing the method to use for the pixel field ('slab' or 'disc')
b_0 - magnetic induction corresponding to a magnetization Mo in T (default: 1)
jacobi - matrix in which to save the jacobi matrix (default: None, faster computation)
Returns:
the phasemap as a 2 dimensional array
'''
# Function for creating the lookup-tables:
def phi_lookup(method, n, m, res, b_0):
if method == 'slab':
def F_h(n, m):
a = np.log(res**2 * (n**2 + m**2))
b = np.arctan(n / m)
return n*a - 2*n + 2*m*b
return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
- F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
elif method == 'disc':
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
# Function for the phase contribution of one pixel:
def phi_mag(i, j):
return (np.cos(beta[j, i]) * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
- np.sin(beta[j, i]) * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
# Function for the derivative of the phase contribution of one pixel:
def phi_mag_deriv(i, j):
return -(np.sin(beta[j, i]) * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+ np.cos(beta[j, i]) * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
# Process input parameters:
v_dim, u_dim = np.shape(projection[0])
v_mag, u_mag = projection
beta = np.arctan2(v_mag, u_mag)
mag = np.hypot(u_mag, v_mag)
coeff = -b_0 * res**2 / (2*PHI_0)
# Create lookup-tables for the phase of one pixel:
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)
phi_u = phi_lookup(method, uu, vv, res, b_0)
phi_v = phi_lookup(method, vv, uu, res, b_0)
# Calculation of the phase:
phase = np.zeros((v_dim, u_dim))
threshold = 0
if jacobi is not None: # With Jacobian matrix (slower)
jacobi[:] = 0 # Jacobi matrix --> zeros
############################### TODO: NUMERICAL CORE #####################################
for j in range(v_dim):
for i in range(u_dim):
phase_cache = phi_mag(i, j)
jacobi[:, i+u_dim*j] = phase_cache.reshape(-1)
if mag[j, i] > threshold:
phase += mag[j, i]*phase_cache
jacobi[:, u_dim*v_dim+i+u_dim*j] = (mag[j, i]*phi_mag_deriv(i, j)).reshape(-1)
############################### TODO: NUMERICAL CORE #####################################
else: # Without Jacobi matrix (faster)
for j in range(v_dim):
for i in range(u_dim):
if abs(mag[j, i]) > threshold:
phase += mag[j, i] * phi_mag(i, j)
# Return the phase:
return phase
#def phase_elec(mag_data, v_0=0, v_acc=30000):
# # TODO: Docstring
#
# res = mag_data.res
# z_dim, y_dim, x_dim = mag_data.dim
# z_mag, y_mag, x_mag = mag_data.magnitude
#
# phase = np.logical_or(x_mag, y_mag, z_mag)
#
# lam = H_BAR / np.sqrt(2 * M_E * v_acc * (1 + Q_E*v_acc / (2*M_E*C**2)))
#
# Ce = (2*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)))
#
# phase *= res * v_0 * Ce
#
# return phase