From 27f46b6fdc0cb8872b29fdd0849adc5482e43ba6 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Mon, 18 Nov 2013 14:52:29 +0100
Subject: [PATCH] holoimage, magdata, phasemap: Made plotting more flexible
(now return axes) phasemapper: new function get_kernel() for convenience,
mag_real_fast() for calculation of the phase with pre-computed lookup-tables
projector: single_tilt_projection() now works as intended Added scripts for
the creation of core-shell-structures Added script for charge calculation by
contour integration Added script for creating tilt series of
core-shell-structure for reconstruction Various minor modifications, mostly
nomenclature and streamlining
---
pyramid/holoimage.py | 24 +-
pyramid/magdata.py | 9 +-
pyramid/phasemap.py | 15 +-
pyramid/phasemapper.py | 371 +++++-------------
pyramid/projector.py | 60 +--
scripts/compare methods/compare_methods.py | 49 ++-
.../compare methods/compare_pixel_fields.py | 75 +++-
.../create_core_shell_disc.py | 81 ++++
.../create_core_shell_sphere.py | 81 ++++
scripts/create distributions/create_sample.py | 20 +-
scripts/create distributions/create_vortex.py | 6 +-
scripts/kernel.py | 52 +++
.../ch5-0-evaluation_and_comparison.py | 40 +-
.../ch5-4-comparison_of_methods_new.py | 306 +++++++++++++++
.../reconstruction/reconstruct_core_shell.py | 94 +++++
scripts/robert_filter.py | 58 +++
scripts/test.py | 27 ++
scripts/test_arb_projection.py | 2 +-
scripts/test_methods.py | 23 +-
scripts/test_vadim.py | 186 +++++++++
20 files changed, 1207 insertions(+), 372 deletions(-)
create mode 100644 scripts/create distributions/create_core_shell_disc.py
create mode 100644 scripts/create distributions/create_core_shell_sphere.py
create mode 100644 scripts/kernel.py
create mode 100644 scripts/paper 1/ch5-4-comparison_of_methods_new.py
create mode 100644 scripts/reconstruction/reconstruct_core_shell.py
create mode 100644 scripts/robert_filter.py
create mode 100644 scripts/test.py
create mode 100644 scripts/test_vadim.py
diff --git a/pyramid/holoimage.py b/pyramid/holoimage.py
index 4ec2fad..639c90d 100644
--- a/pyramid/holoimage.py
+++ b/pyramid/holoimage.py
@@ -17,7 +17,7 @@ from numpy import pi
import matplotlib as mpl
import matplotlib.pyplot as plt
-from matplotlib.ticker import NullLocator
+from matplotlib.ticker import NullLocator, MaxNLocator
from PIL import Image
from pyramid.phasemap import PhaseMap
@@ -61,12 +61,10 @@ def holo_image(phase_map, density=1):
'''
assert isinstance(phase_map, PhaseMap), 'phase_map has to be a PhaseMap object!'
- # Extract the phase (considering the unit):
- phase = phase_map.phase
# Calculate the holography image intensity:
- img_holo = (1 + np.cos(density * phase)) / 2
+ img_holo = (1 + np.cos(density * phase_map.phase)) / 2
# Calculate the phase gradients, expressed by magnitude and angle:
- phase_grad_y, phase_grad_x = np.gradient(phase, phase_map.res, phase_map.res)
+ phase_grad_y, phase_grad_x = np.gradient(phase_map.phase, phase_map.res, phase_map.res)
phase_angle = (1 - np.arctan2(phase_grad_y, phase_grad_x)/pi) / 2
phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)
if phase_magnitude.max() != 0:
@@ -126,7 +124,8 @@ def display(holo_image, title='Holographic Contour Map', axis=None, interpolatio
Returns
-------
- None
+ axis: :class:`~matplotlib.axes.AxesSubplot`
+ The axis on which the graph is plotted.
'''
# If no axis is specified, a new figure is created:
@@ -141,6 +140,10 @@ def display(holo_image, title='Holographic Contour Map', axis=None, interpolatio
axis.set_title(title, fontsize=18)
axis.set_xlabel('x-axis [px]', fontsize=15)
axis.set_ylabel('y-axis [px]', fontsize=15)
+ axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+ axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+
+ return axis
def display_combined(phase_map, density=1, title='Combined Plot', interpolation='none'):
@@ -160,7 +163,8 @@ def display_combined(phase_map, density=1, title='Combined Plot', interpolation=
Returns
-------
- None
+ phase_axis, holo_axis: :class:`~matplotlib.axes.AxesSubplot`
+ The axes on which the graphs are plotted.
'''
# Create combined plot and set title:
@@ -169,7 +173,13 @@ def display_combined(phase_map, density=1, title='Combined Plot', interpolation=
# Plot holography image:
holo_axis = fig.add_subplot(1, 2, 1, aspect='equal')
display(holo_image(phase_map, density), axis=holo_axis, interpolation=interpolation)
+ holo_axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+ holo_axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect='equal')
fig.subplots_adjust(right=0.85)
phase_map.display(axis=phase_axis)
+ phase_axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+ phase_axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+
+ return phase_axis, holo_axis
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index d87b8ed..792ca02 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -323,7 +323,8 @@ class MagData:
Returns
-------
- None
+ axis: :class:`~matplotlib.axes.AxesSubplot`
+ The axis on which the graph is plotted.
'''
assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
@@ -361,8 +362,8 @@ class MagData:
axis.set_xlabel(u_label, fontsize=15)
axis.set_ylabel(v_label, fontsize=15)
axis.tick_params(axis='both', which='major', labelsize=14)
- axis.xaxis.set_major_locator(MaxNLocator(nbins=8, integer=True))
- axis.yaxis.set_major_locator(MaxNLocator(nbins=8, integer=True))
+ axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+ axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
plt.show()
def quiver_plot3d(self):
@@ -391,7 +392,7 @@ class MagData:
z_mag = np.reshape(self.magnitude[0], (-1))
# Plot them as vectors:
mlab.figure()
- plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow', scale_factor=10.0)
+ plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow')#, scale_factor=5.0)
mlab.outline(plot)
mlab.axes(plot)
mlab.colorbar()
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index bd7cf74..9062364 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -13,6 +13,7 @@ import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
+from matplotlib.ticker import MaxNLocator, FuncFormatter
import netCDF4
@@ -194,7 +195,8 @@ class PhaseMap:
Returns
-------
- None
+ axis: :class:`~matplotlib.axes.AxesSubplot`
+ The axis on which the graph is plotted.
'''
# Take units into consideration:
@@ -210,11 +212,11 @@ class PhaseMap:
# Set the axes ticks and labels:
axis.set_xlim(0, np.shape(phase)[1])
axis.set_ylim(0, np.shape(phase)[0])
- ticks = (axis.get_xticks()*self.res).astype(int)
- axis.set_xticklabels(ticks)
- ticks = (axis.get_yticks()*self.res).astype(int)
+ axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+ axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+ axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.res)))
+ axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.res)))
axis.tick_params(axis='both', which='major', labelsize=14)
- axis.set_yticklabels(ticks)
axis.set_title(title, fontsize=18)
axis.set_xlabel('x [nm]', fontsize=15)
axis.set_ylabel('y [nm]', fontsize=15)
@@ -241,7 +243,8 @@ class PhaseMap:
Returns
-------
- None
+ axis: :class:`~matplotlib.axes.AxesSubplot`
+ The axis on which the graph is plotted.
'''
# Take units into consideration:
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 8a27ccc..58413ea 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -62,12 +62,12 @@ def phase_mag_fourier(res, projection, padding=0, b_0=1):
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)
+ nyq = 0.5 / res # nyquist frequency
+ f_u = np.linspace(0, nyq, u_mag_fft.shape[1])
+ f_v = np.linspace(-nyq, nyq, 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)
+ phase_fft = coeff * res * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30) #* (8*(res/2)**3)*np.sinc(res/2*f_uu)*np.sinc(res/2*f_vv) * np.exp(2*pi*1j*())
# 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]
@@ -100,54 +100,36 @@ def phase_mag_real(res, projection, method='discs', b_0=1, jacobi=None):
The phase 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])
+ # Process input parameters:
+ dim = np.shape(projection[0])
v_mag, u_mag = projection[:-1]
- 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)
- u_phi = phi_lookup(method, uu, vv, res, b_0)
- v_phi = phi_lookup(method, vv, uu, res, b_0)
+ # Get lookup-tables for the phase of one pixel:
+ u_phi = get_kernel(method, 'u', dim, res, b_0)
+ v_phi = get_kernel(method, 'v', dim, res, b_0)
# Calculation of the phase:
- phase = np.zeros((v_dim, u_dim))
+ phase = np.zeros(dim)
threshold = 0
if jacobi is not None: # With Jacobian matrix (slower)
jacobi[:] = 0 # Jacobi matrix --> zeros
- for j in range(v_dim):
- for i in range(u_dim):
- u_phase = u_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
- jacobi[:, i+u_dim*j] = u_phase.reshape(-1)
+ for j in range(dim[0]):
+ for i in range(dim[1]):
+ u_phase = u_phi[dim[0]-1-j:(2*dim[0]-1)-j, dim[1]-1-i:(2*dim[1]-1)-i]
+ jacobi[:, i+dim[1]*j] = u_phase.reshape(-1)
if abs(u_mag[j, i]) > threshold:
phase += u_mag[j, i] * u_phase
- v_phase = v_phi[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] = -v_phase.reshape(-1)
+ v_phase = v_phi[dim[0]-1-j:(2*dim[0]-1)-j, dim[1]-1-i:(2*dim[1]-1)-i]
+ jacobi[:, dim[1]*dim[0]+i+dim[1]*j] = -v_phase.reshape(-1)
if abs(v_mag[j, i]) > threshold:
phase -= v_mag[j, i] * v_phase
else: # Without Jacobi matrix (faster)
- nc.phase_mag_real_core(v_dim, u_dim, v_phi, u_phi, v_mag, u_mag, phase, threshold)
+ nc.phase_mag_real_core(dim[0], dim[1], v_phi, u_phi, v_mag, u_mag, phase, threshold)
# Return the phase:
return phase
-def phase_mag_real_conv(res, projection, method='disc', b_0=1, jacobi=None):
+def phase_mag_real_conv(res, projection, method='disc', b_0=1):
'''Calculate the magnetic phase from magnetization data (real space approach).
Parameters
@@ -173,116 +155,97 @@ def phase_mag_real_conv(res, projection, method='disc', b_0=1, jacobi=None):
The phase as a 2-dimensional array.
''' # TODO: Docstring!!!
+ # Process input parameters:
+ dim = np.shape(projection[0])
+ v_mag, u_mag = projection[:-1]
- # Function for creating the lookup-tables:
- def phi_lookup(method, n, m, res, b_0):
+ # Get lookup-tables for the phase of one pixel:
+ u_phi = get_kernel(method, 'u', dim, res, b_0)
+ v_phi = get_kernel(method, 'v', dim, res, b_0)
+
+ # Return the phase:
+ result = fftconvolve(u_mag, u_phi, 'same') - fftconvolve(v_mag, v_phi, 'same')
+ return result
+
+
+def phase_mag_real_fast(res, projection, kernels_fourier, b_0=1):
+ '''Calculate the magnetic phase from magnetization data (real space approach).
+
+ Parameters
+ ----------
+ res : float
+ The resolution of the grid (grid spacing) in nm.
+ projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
+ The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
+ of the magnetization and the thickness projection for the resulting 2D-grid.
+ method : {'disc', 'slab'}, optional
+ Specifies the elemental geometry to use for the pixel field.
+ The default is 'disc', because of the smaller computational overhead.
+ b_0 : float, optional
+ The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+ The default is 1.
+ jacobi : :class:`~numpy.ndarray` (N=2), optional
+ Specifies the matrix in which to save the jacobi matrix. The jacobi matrix will not be
+ calculated, if no matrix is specified (default), resulting in a faster computation.
+
+ Returns
+ -------
+ phase : :class:`~numpy.ndarray` (N=2)
+ The phase as a 2-dimensional array.
+
+ ''' # TODO: Docstring!!!
+ # Process input parameters:
+ v_mag, u_mag = projection[:-1]
+ dim = np.shape(u_mag)
+
+ size = 3*np.array(dim) - 1 # dim + (2*dim - 1) magnetisation + kernel
+ fsize = 2 ** np.ceil(np.log2(size)).astype(int)
+ fslice = tuple([slice(0, int(sz)) for sz in size])
+
+ u_mag_f = np.fft.rfftn(u_mag, fsize)
+ v_mag_f = np.fft.rfftn(v_mag, fsize)
+
+ v_kern_f, u_kern_f = kernels_fourier
+
+ fslice = (slice(dim[0]-1, 2*dim[0]-1), slice(dim[1]-1, 2*dim[1]-1))
+ u_phase = np.fft.irfftn(u_mag_f * u_kern_f, fsize)[fslice].copy()
+ v_phase = np.fft.irfftn(v_mag_f * v_kern_f, fsize)[fslice].copy()
+ return u_phase - v_phase
+
+
+def get_kernel_fourier(method, orientation, dim, res, b_0=1):
+
+ kernel = get_kernel(method, orientation, dim, res, b_0)
+
+ size = 3*np.array(dim) - 1 # dim + (2*dim - 1) magnetisation + kernel
+ fsize = 2 ** np.ceil(np.log2(size)).astype(int)
+
+ return np.fft.rfftn(kernel, fsize)
+
+
+def get_kernel(method, orientation, dim, res, b_0=1):
+
+ def get_elementary_phase(method, n, m, res):
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))
+ return 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[:-1]
+ return m / (n**2 + m**2 + 1E-30) * in_or_out
+
coeff = -b_0 * res**2 / (2*PHI_0)
-
- # Create lookup-tables for the phase of one pixel:
+ v_dim, u_dim = dim
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)
- u_phi = phi_lookup(method, uu, vv, res, b_0)
- v_phi = phi_lookup(method, vv, uu, res, b_0)
-
- # Return the phase:
- return fftconvolve(u_mag, u_phi, 'same') - fftconvolve(v_mag, v_phi, 'same')
-
-# def fftconvolve(in1, in2, mode="full"):
-# """Convolve two N-dimensional arrays using FFT.
-#
-# Convolve `in1` and `in2` using the fast Fourier transform method, with
-# the output size determined by the `mode` argument.
-#
-# This is generally much faster than `convolve` for large arrays (n > ~500),
-# but can be slower when only a few output values are needed, and can only
-# output float arrays (int or object array inputs will be cast to float).
-#
-# Parameters
-# ----------
-# in1 : array_like
-# First input.
-# in2 : array_like
-# Second input. Should have the same number of dimensions as `in1`;
-# if sizes of `in1` and `in2` are not equal then `in1` has to be the
-# larger array.
-# mode : str {'full', 'valid', 'same'}, optional
-# A string indicating the size of the output:
-#
-# ``full``
-# The output is the full discrete linear convolution
-# of the inputs. (Default)
-# ``valid``
-# The output consists only of those elements that do not
-# rely on the zero-padding.
-# ``same``
-# The output is the same size as `in1`, centered
-# with respect to the 'full' output.
-#
-# Returns
-# -------
-# out : array
-# An N-dimensional array containing a subset of the discrete linear
-# convolution of `in1` with `in2`.
-#
-# """
-# in1 = asarray(in1)
-# in2 = asarray(in2)
-#
-# if rank(in1) == rank(in2) == 0: # scalar inputs
-# return in1 * in2
-# elif not in1.ndim == in2.ndim:
-# raise ValueError("in1 and in2 should have the same rank")
-# elif in1.size == 0 or in2.size == 0: # empty arrays
-# return array([])
-#
-# s1 = array(in1.shape)
-# s2 = array(in2.shape)
-# complex_result = (np.issubdtype(in1.dtype, np.complex) or
-# np.issubdtype(in2.dtype, np.complex))
-# size = s1 + s2 - 1
-#
-# if mode == "valid":
-# for d1, d2 in zip(s1, s2):
-# if not d1 >= d2:
-# warnings.warn("in1 should have at least as many items as in2 in "
-# "every dimension for 'valid' mode. In scipy "
-# "0.13.0 this will raise an error",
-# DeprecationWarning)
-#
-# # Always use 2**n-sized FFT
-# fsize = 2 ** np.ceil(np.log2(size)).astype(int)
-# print('fsize =', fsize)
-# print('s1 =', s1)
-# print('s2 =', s2)
-# fslice = tuple([slice(0, int(sz)) for sz in size])
-# if not complex_result:
-# ret = irfftn(rfftn(in1, fsize) *
-# rfftn(in2, fsize), fsize)[fslice].copy()
-# ret = ret.real
-# else:
-# ret = ifftn(fftn(in1, fsize) * fftn(in2, fsize))[fslice].copy()
-#
-# if mode == "full":
-# return ret
-# elif mode == "same":
-# return _centered(ret, s1)
-# elif mode == "valid":
-# return _centered(ret, abs(s1 - s2) + 1)
+ if orientation == 'u':
+ return coeff * get_elementary_phase(method, uu, vv, res)
+ elif orientation == 'v':
+ return coeff * get_elementary_phase(method, vv, uu, res)
def phase_elec(res, projection, v_0=1, v_acc=30000):
@@ -311,143 +274,3 @@ def phase_elec(res, projection, v_0=1, v_acc=30000):
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))
# return phase:
return v_0 * Ce * projection[-1] * res*1E-9
-
-
-#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
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 326d895..00d56c8 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -8,11 +8,12 @@ directions with the use of transfer functions (work in progress).
"""
-import time
import numpy as np
from numpy import pi
+import matplotlib.pyplot as plt
+from matplotlib.ticker import MaxNLocator
import itertools
from pyramid.magdata import MagData
@@ -59,26 +60,21 @@ def simple_axis_projection(mag_data, axis='z', threshold=0):
def single_tilt_projection(mag_data, tilt=0, threshold=0):
# TODO: Docstring!!!
+ '''Projection with tilt around the (centered!) y-axis'''
assert isinstance(mag_data, MagData), 'Parameter mag_data has to be a MagData object!'
# assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as string)!'
- # Set starting variables:
- dim = (mag_data.dim[0], mag_data.dim[2])
- z_mag, y_mag, x_mag = mag_data.magnitude
- mask = mag_data.get_mask()
- projection = (np.zeros((mag_data.dim[1], mag_data.dim[2])),
- np.zeros((mag_data.dim[1], mag_data.dim[2])),
- np.zeros((mag_data.dim[1], mag_data.dim[2])))
-
def get_position(p, m, b, size):
- x, y = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
+ y, x = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
+# print 'p:', p
+# print 'position:', (y-m*x-b)/np.sqrt(m**2+1)# + size/2.
return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
def get_impact(pos, r, size):
return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
if 0 <= x < size]
- def get_weight(delta, rho):
+ def get_weight(delta, rho): # use circles to represent the voxels
a, b = delta-rho, delta+rho
if a >= 1 or b <= -1: # TODO: Should not be necessary!
print 'Outside of bounds:', delta
@@ -95,18 +91,27 @@ def single_tilt_projection(mag_data, tilt=0, threshold=0):
w_a = (a*np.sqrt(1-a**2) + np.arctan(a/np.sqrt(1-a**2))) / pi
return w_b - w_a
+ # Set starting variables:
+ # length along projection (proj, z), rotation (rot, y) and perpendicular (perp, x) axis:
+ dim_proj, dim_rot, dim_perp = mag_data.dim
+ z_mag, y_mag, x_mag = mag_data.magnitude
+ mask = mag_data.get_mask()
+ projection = (np.zeros((dim_rot, dim_perp)),
+ np.zeros((dim_rot, dim_perp)),
+ np.zeros((dim_rot, dim_perp)))
# Creating coordinate list of all voxels:
- xi = range(dim[0])
- yj = range(dim[1])
- ii, jj = np.meshgrid(xi, yj)
- voxels = list(itertools.product(yj, xi))
+ voxels = list(itertools.product(range(dim_proj), range(dim_perp)))
+# print 'voxels:', voxels
# Calculate positions along the projected pixel coordinate system:
- direction = (-np.cos(tilt), np.sin(tilt))
- center = (dim[0]/2., dim[1]/2.)
- m = direction[0]/(direction[1]+1E-30)
+ center = (dim_proj/2., dim_perp/2.)
+# direction = (-np.cos(tilt), np.sin(tilt))
+ m = np.where(tilt<=pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))#direction[0]/(direction[1]+1E-30)
+# print 'm:', m
b = center[0] - m * center[1]
- positions = get_position(voxels, m, b, dim[0])
+# print 'b:', b
+ positions = get_position(voxels, m, b, dim_perp)
+# print positions
# Calculate weights:
r = 1/np.sqrt(np.pi) # radius of the voxel circle
@@ -114,19 +119,18 @@ def single_tilt_projection(mag_data, tilt=0, threshold=0):
weights = []
for i, voxel in enumerate(voxels):
voxel_weights = []
- impacts = get_impact(positions[i], r, dim[0])
+ impacts = get_impact(positions[i], r, dim_perp)
for impact in impacts:
distance = np.abs(impact+0.5 - positions[i])
delta = distance / r
voxel_weights.append((impact, get_weight(delta, rho)))
weights.append((voxel, voxel_weights))
-
- # Calculate projection with the calculated weights for the voxels:
- for i, weight in enumerate(weights):
- voxel = weights[i][0]
- voxel_weights = weights[i][1]
- for voxel_weight in voxel_weights:
- pixel, weight = voxel_weight
+
+# print weights
+
+ # Calculate projection with the calculated weights for the voxels (rotation around y-axis):
+ for i, (voxel, voxel_weights) in enumerate(weights):
+ for pixel, weight in voxel_weights:
# Component parallel to tilt axis (':' goes over all slices):
projection[0][:, pixel] += weight * y_mag[voxel[0], :, voxel[1]]
# Component perpendicular to tilt axis:
@@ -134,5 +138,5 @@ def single_tilt_projection(mag_data, tilt=0, threshold=0):
+ z_mag[voxel[0], :, voxel[1]]*np.sin(tilt))
# Thickness profile:
projection[2][:, pixel] += weight * mask[voxel[0], :, voxel[1]]
-
+
return projection
diff --git a/scripts/compare methods/compare_methods.py b/scripts/compare methods/compare_methods.py
index 5fc6801..dd32b50 100644
--- a/scripts/compare methods/compare_methods.py
+++ b/scripts/compare methods/compare_methods.py
@@ -29,14 +29,14 @@ def compare_methods():
'''
# Input parameters:
- b_0 = 1.1 # in T
- res = 10.0 # in nm
+ b_0 = 1.0 # in T
+ res = 1.0 # in nm
phi = pi/4
padding = 3
- density = 1
- dim = (8, 512, 512) # in px (z, y, x)
+ density = 10
+ dim = (128, 128, 128) # in px (z, y, x)
# Create magnetic shape:
- geometry = 'slab'
+ geometry = 'sphere'
if geometry == 'slab':
center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z,y,x) index starts at 0!
width = (dim[0]/2, dim[1]/2, dim[2]/2) # in px (z, y, x)
@@ -49,31 +49,48 @@ def compare_methods():
mag_shape = mc.Shapes.disc(dim, center, radius, height)
phase_ana = an.phase_mag_disc(dim, res, phi, center, radius, height, b_0)
elif geometry == 'sphere':
- center = (50, 50, 50) # in px (z, y, x) index starts with 0!
- radius = 25 # in px
+ center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.50) # in px (z, y, x) index starts with 0!
+ radius = dim[1]/4 # in px
mag_shape = mc.Shapes.sphere(dim, center, radius)
phase_ana = an.phase_mag_sphere(dim, res, phi, center, radius, b_0)
# Project the magnetization data:
mag_data = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
mag_data.quiver_plot(ax_slice=int(center[0]))
projection = pj.simple_axis_projection(mag_data)
+ pj.quiver_plot(projection)
# Construct phase maps:
phase_map_ana = PhaseMap(res, phase_ana)
start_time = time.time()
phase_map_fft = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding, b_0))
- print 'Time for Fourier space approach: ', time.time() - start_time
+ print 'Time for Fourier space approach: ', time.time() - start_time
start_time = time.time()
phase_map_slab = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab', b_0))
- print 'Time for real space approach (Slab):', time.time() - start_time
+ print 'Time for real space approach (Slab): ', time.time() - start_time
start_time = time.time()
phase_map_disc = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc', b_0))
- print 'Time for real space approach (Disc):', time.time() - start_time
+ print 'Time for real space approach (Disc): ', time.time() - start_time
start_time = time.time()
phase_map_slab_conv = PhaseMap(res, pm.phase_mag_real_conv(res, projection, 'slab', b_0))
- print 'Time for real space approach (Slab Convolve):', time.time() - start_time
+ print 'Time for real space approach (Slab Conv):', time.time() - start_time
+ start_time = time.time()
+ phase_map_disc_conv = PhaseMap(res, pm.phase_mag_real_conv(res, projection,'disc', b_0))
+ print 'Time for real space approach (Disc Conv):', time.time() - start_time
+# start_time = time.time()
+ u_kern_f = pm.get_kernel_fourier('slab', 'u', np.shape(projection[0]), res, b_0)
+ v_kern_f = pm.get_kernel_fourier('slab', 'v', np.shape(projection[0]), res, b_0)
+# print 'Time for kernel (Slab):', time.time() - start_time
+ start_time = time.time()
+ phase = pm.phase_mag_real_fast(res, projection, (v_kern_f, u_kern_f), b_0)
+ phase_map_slab_fast = PhaseMap(res, phase)
+ print 'Time for real space approach (Slab Fast):', time.time() - start_time
+# start_time = time.time()
+ u_kern_f = pm.get_kernel_fourier('disc', 'u', np.shape(projection[0]), res, b_0)
+ v_kern_f = pm.get_kernel_fourier('disc', 'v', np.shape(projection[0]), res, b_0)
+# print 'Time for kernel (Disc):', time.time() - start_time
start_time = time.time()
- phase_map_disc_conv = PhaseMap(res, pm.phase_mag_real_conv(res, projection, 'disc', b_0))
- print 'Time for real space approach (Disc Convolve):', time.time() - start_time
+ phase = pm.phase_mag_real_fast(res, projection, (v_kern_f, u_kern_f), b_0)
+ phase_map_disc_fast = PhaseMap(res, phase)
+ print 'Time for real space approach (Disc Fast):', time.time() - start_time
# Display the combinated plots with phasemap and holography image:
hi.display_combined(phase_map_ana, density, 'Analytic Solution')
hi.display_combined(phase_map_fft, density, 'Fourier Space')
@@ -81,6 +98,8 @@ def compare_methods():
hi.display_combined(phase_map_disc, density, 'Real Space (Disc)')
hi.display_combined(phase_map_slab_conv, density, 'Real Space (Slab Convolve)')
hi.display_combined(phase_map_disc_conv, density, 'Real Space (Disc Convolve)')
+ hi.display_combined(phase_map_slab_fast, density, 'Real Space (Slab Fast)')
+ hi.display_combined(phase_map_disc_fast, density, 'Real Space (Disc Disc)')
# Plot differences to the analytic solution:
phase_map_diff_fft = PhaseMap(res, phase_map_ana.phase-phase_map_fft.phase)
phase_map_diff_slab = PhaseMap(res, phase_map_ana.phase-phase_map_slab.phase)
@@ -95,8 +114,8 @@ def compare_methods():
phase_map_diff_fft.display('Fourier Space (RMS = {:3.2e})'.format(RMS_fft))
phase_map_diff_slab.display('Real Space (Slab) (RMS = {:3.2e})'.format(RMS_slab))
phase_map_diff_disc.display('Real Space (Disc) (RMS = {:3.2e})'.format(RMS_disc))
- phase_map_diff_slab_conv.display('Real Space (Disc) (RMS = {:3.2e})'.format(RMS_slab_conv))
- phase_map_diff_disc_conv.display('Real Space (Disc) (RMS = {:3.2e})'.format(RMS_disc_conv))
+ phase_map_diff_slab_conv.display('Real Conv. (Slab) (RMS = {:3.2e})'.format(RMS_slab_conv))
+ phase_map_diff_disc_conv.display('Real Conv. (Disc) (RMS = {:3.2e})'.format(RMS_disc_conv))
if __name__ == "__main__":
diff --git a/scripts/compare methods/compare_pixel_fields.py b/scripts/compare methods/compare_pixel_fields.py
index fd8c8df..414ac02 100644
--- a/scripts/compare methods/compare_pixel_fields.py
+++ b/scripts/compare methods/compare_pixel_fields.py
@@ -16,6 +16,8 @@ import pyramid.phasemapper as pm
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
+from scipy.signal import fftconvolve
+
def compare_pixel_fields():
'''Calculate and display the phase map for different real space approaches.
@@ -43,7 +45,7 @@ def compare_pixel_fields():
f_u = np.linspace(0, nyq/2, dim[2]/2.+1)
f_v = np.linspace(-nyq/2, nyq/2, dim[1], endpoint=False)
f_uu, f_vv = np.meshgrid(f_u, f_v)
- phase_fft = coeff * f_vv / (f_uu**2 + f_vv**2 + 1e-30)
+ phase_fft = coeff * f_vv / (f_uu**2 + f_vv**2 + 1e-30) #* (8*(res/2)**3)*np.sinc(res/2*f_uu)*np.sinc(res/2*f_vv)
# Transform to real space and revert padding:
phase_fft = np.fft.ifftshift(phase_fft, axes=0)
phase_fft_kernel = np.fft.fftshift(np.fft.irfft2(phase_fft), axes=(0, 1))
@@ -62,30 +64,75 @@ def compare_pixel_fields():
# Kernel of the Fourier method:
phase_map_fft = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding=0), 'mrad')
phase_map_fft.display('Phase of one Pixel (Fourier)', limit=limit)
+
# Kernel of the Fourier method, calculated directly:
phase_map_fft_kernel = PhaseMap(res, get_fourier_kernel(), 'mrad')
phase_map_fft_kernel.display('Phase of one Pixel (Fourier Kernel)', limit=limit)
# Kernel differences:
print 'Fourier Kernel, direct and indirect method are identical:', \
np.all(phase_map_fft_kernel.phase - phase_map_fft.phase) == 0
- phase_map_diff = PhaseMap(res, phase_map_fft.phase - phase_map_disc.phase, 'mrad')
+ phase_map_diff = PhaseMap(res, phase_map_disc.phase - phase_map_fft.phase, 'mrad')
phase_map_diff.display('Phase difference of one Pixel (Disc - Fourier)')
- phase_inv_fft_r = np.real(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))
- phase_map_inv_fft = PhaseMap(res, phase_inv_fft_r)
- phase_map_inv_fft.display('FT of the Phase of one Pixel (Fourier, real)')
+ phase_inv_fft = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))**2
+ phase_map_inv_fft = PhaseMap(res, phase_inv_fft)
+ phase_map_inv_fft.display('FT of the Phase of one Pixel (Fourier, Power)')
+
+ phase_inv_disc = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))**2
+ phase_map_inv_disc = PhaseMap(res, phase_inv_disc)
+ phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, Power)')
+
+ import matplotlib.pyplot as plt
+ from matplotlib.ticker import IndexLocator
+
+ fig = plt.figure()
+ axis = fig.add_subplot(1, 1, 1)
+ x = range(dim[1])
+ y_ft = phase_map_fft.phase[:, dim[1]/2]
+ y_re = phase_map_disc.phase[:, dim[1]/2]
+ axis.axhline(0, color='k')
+ axis.plot(x, y_re, label='Real space method')
+ axis.plot(x, y_ft, label='Fourier method')
+# axis.plot(x, y_re-y_ft, label='Difference')
+ axis.grid()
+ axis.legend()
+ axis.set_title('Real Space Kernel')
+ axis.set_xlim(0, dim[1]-1)
+ axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
+
+ fig = plt.figure()
+ axis = fig.add_subplot(1, 1, 1)
+ x = range(dim[1])
+ k_re = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))[:, 0]**2
+ k_ft = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))[:, 0]**2
+
+
+
+ phase = np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0)
+ nyq = 0.5 / res
+ f_u = np.linspace(0, nyq, dim[2]/2.+1)
+ f_v = np.linspace(-nyq, nyq, dim[1], endpoint=False)
+ f_uu, f_vv = np.meshgrid(f_u, f_v)
+ kern = 4*(res/2)**2*np.sinc(f_uu*res/2)*np.sinc(f_vv*res/2)
+ phase_mod = fftconvolve(phase, kern, 'same')
+ k_ft_mod = np.abs(phase_mod)[:, 0]**2
+# sincy = np.sinc((np.asarray(x)-(dim[1]/2-1))*dim[1])
+# quark
- phase_inv_fft_i = np.imag(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))
- phase_map_inv_fft = PhaseMap(res, phase_inv_fft_i)
- phase_map_inv_fft.display('FT of the Phase of one Pixel (Fourier, imag)')
- phase_inv_disc_r = np.real(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))
- phase_map_inv_disc = PhaseMap(res, phase_inv_disc_r)
- phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, real)')
- phase_inv_disc_i = np.imag(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))
- phase_map_inv_disc = PhaseMap(res, phase_inv_disc_i)
- phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, imag)')
+ axis.axhline(0, color='k')
+ axis.plot(x, k_re, label='Real space method')
+ axis.plot(x, k_ft, label='Fourier method')
+# axis.plot(x, k_ft_mod, label='Fourier modified')
+ axis.grid()
+ axis.legend()
+ axis.set_title('Fourier Space Kernel')
+ axis.set_xlim(0, dim[1]-1)
+ axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
+
+# # Convolve:
+# phase_fft = fftconvolve(phase_fft, np.sinc(f_uu/(nyq/2))*np.sinc(f_vv/nyq), 'same')
if __name__ == "__main__":
diff --git a/scripts/create distributions/create_core_shell_disc.py b/scripts/create distributions/create_core_shell_disc.py
new file mode 100644
index 0000000..0242008
--- /dev/null
+++ b/scripts/create distributions/create_core_shell_disc.py
@@ -0,0 +1,81 @@
+# -*- coding: utf-8 -*-
+"""Create a core-shell disc."""
+
+
+import pdb
+import traceback
+import sys
+import os
+
+import numpy as np
+
+import matplotlib.pyplot as plt
+
+import pyramid.magcreator as mc
+import pyramid.projector as pj
+import pyramid.phasemapper as pm
+import pyramid.holoimage as hi
+from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
+
+
+def create_core_shell_disc():
+ '''Calculate and display the Pyramid-Logo.
+ Arguments:
+ None
+ Returns:
+ None
+
+ '''
+ directory = '../../output/magnetic distributions'
+ if not os.path.exists(directory):
+ os.makedirs(directory)
+ # Input parameters:
+ filename = directory + '/mag_dist_core_shell_disc.txt'
+ res = 1.0 # in nm
+ density = 1
+ dim = (64, 64, 64)
+ center = (dim[0]/2-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
+ radius_core = dim[1]/8
+ radius_shell = dim[1]/4
+ height = dim[0]/2
+ # Create magnetic shape:
+ mag_shape_core = mc.Shapes.disc(dim, center, radius_core, height)
+ mag_shape_outer = mc.Shapes.disc(dim, center, radius_shell, height)
+ mag_shape_shell = np.logical_xor(mag_shape_outer, mag_shape_core)
+ mag_data = MagData(res, mc.create_mag_dist_vortex(mag_shape_shell, magnitude=0.75))
+ mag_data.add_magnitude(mc.create_mag_dist_homog(mag_shape_core, phi=0, theta=0))
+ mag_data.quiver_plot('z-projection', proj_axis='z')
+ mag_data.quiver_plot('x-projection', proj_axis='x')
+ mag_data.quiver_plot3d()
+ mag_data.save_to_llg(filename)
+ projection_z = pj.simple_axis_projection(mag_data, axis='z')
+ projection_x = pj.simple_axis_projection(mag_data, axis='x')
+ phase_map_z = PhaseMap(res, pm.phase_mag_real(res, projection_z, 'slab'))
+ phase_map_x = PhaseMap(res, pm.phase_mag_real(res, projection_x, 'slab'))
+ hi.display_combined(phase_map_z, density, 'Core-Shell structure (z-projection)')
+ phase_axis, holo_axis = hi.display_combined(phase_map_x, density,
+ 'Core-Shell structure (x-projection)')
+ phase_axis.set_xlabel('y [nm]')
+ phase_axis.set_ylabel('z [nm]')
+ holo_axis.set_xlabel('y-axis [px]')
+ holo_axis.set_ylabel('z-axis [px]')
+ phase_slice_z = phase_map_z.phase[dim[1]/2, :]
+ phase_slice_x = phase_map_x.phase[dim[0]/2, :]
+ fig = plt.figure()
+ fig.add_subplot(1, 1, 1)
+ plt.plot(range(dim[2]), phase_slice_z)
+ plt.title('Phase slice along x for z-projection')
+ fig = plt.figure()
+ fig.add_subplot(1, 1, 1)
+ plt.plot(range(dim[1]), phase_slice_x)
+ plt.title('Phase slice along y for x-projection')
+
+
+if __name__ == "__main__":
+ try:
+ create_core_shell_disc()
+ except:
+ type, value, tb = sys.exc_info()
+ traceback.print_exc()
+ pdb.post_mortem(tb)
diff --git a/scripts/create distributions/create_core_shell_sphere.py b/scripts/create distributions/create_core_shell_sphere.py
new file mode 100644
index 0000000..dd6c0aa
--- /dev/null
+++ b/scripts/create distributions/create_core_shell_sphere.py
@@ -0,0 +1,81 @@
+# -*- coding: utf-8 -*-
+"""Create a core-shell disc."""
+
+
+import pdb
+import traceback
+import sys
+import os
+
+import numpy as np
+
+import matplotlib.pyplot as plt
+
+import pyramid.magcreator as mc
+import pyramid.projector as pj
+import pyramid.phasemapper as pm
+import pyramid.holoimage as hi
+from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
+
+
+def create_core_shell_sphere():
+ '''Calculate and display the Pyramid-Logo.
+ Arguments:
+ None
+ Returns:
+ None
+
+ '''
+ directory = '../../output/magnetic distributions'
+ if not os.path.exists(directory):
+ os.makedirs(directory)
+ # Input parameters:
+ filename = directory + '/mag_dist_core_shell_sphere.txt'
+ res = 1.0 # in nm
+ density = 1
+ dim = (32, 32, 32)
+ center = (dim[0]/2-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
+ radius_core = dim[1]/8
+ radius_shell = dim[1]/4
+ # Create magnetic shape:
+ mag_shape_sphere = mc.Shapes.sphere(dim, center, radius_shell)
+ mag_shape_disc = mc.Shapes.disc(dim, center, radius_core, height=dim[0])
+ mag_shape_core = np.logical_and(mag_shape_sphere, mag_shape_disc)
+ mag_shape_shell = np.logical_and(mag_shape_sphere, np.logical_not(mag_shape_core))
+ mag_data = MagData(res, mc.create_mag_dist_vortex(mag_shape_shell, magnitude=0.75))
+ mag_data.add_magnitude(mc.create_mag_dist_homog(mag_shape_core, phi=0, theta=0))
+ mag_data.quiver_plot('z-projection', proj_axis='z')
+ mag_data.quiver_plot('x-projection', proj_axis='x')
+ mag_data.quiver_plot3d()
+ mag_data.save_to_llg(filename)
+ projection_z = pj.simple_axis_projection(mag_data, axis='z')
+ projection_x = pj.simple_axis_projection(mag_data, axis='x')
+ phase_map_z = PhaseMap(res, pm.phase_mag_real(res, projection_z, 'slab'))
+ phase_map_x = PhaseMap(res, pm.phase_mag_real(res, projection_x, 'slab'))
+ hi.display_combined(phase_map_z, density, 'Core-Shell structure (z-projection)')
+ phase_axis, holo_axis = hi.display_combined(phase_map_x, density,
+ 'Core-Shell structure (x-projection)')
+ phase_axis.set_xlabel('y [nm]')
+ phase_axis.set_ylabel('z [nm]')
+ holo_axis.set_xlabel('y-axis [px]')
+ holo_axis.set_ylabel('z-axis [px]')
+ phase_slice_z = phase_map_z.phase[dim[1]/2, :]
+ phase_slice_x = phase_map_x.phase[dim[0]/2, :]
+ fig = plt.figure()
+ fig.add_subplot(1, 1, 1)
+ plt.plot(range(dim[2]), phase_slice_z)
+ plt.title('Phase slice along x for z-projection')
+ fig = plt.figure()
+ fig.add_subplot(1, 1, 1)
+ plt.plot(range(dim[1]), phase_slice_x)
+ plt.title('Phase slice along y for x-projection')
+
+
+if __name__ == "__main__":
+ try:
+ create_core_shell_sphere()
+ except:
+ type, value, tb = sys.exc_info()
+ traceback.print_exc()
+ pdb.post_mortem(tb)
diff --git a/scripts/create distributions/create_sample.py b/scripts/create distributions/create_sample.py
index 2941f57..cfc625a 100644
--- a/scripts/create distributions/create_sample.py
+++ b/scripts/create distributions/create_sample.py
@@ -26,18 +26,18 @@ def create_sample():
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
- key = 'sphere'
+ key = 'pixel'
filename = directory + '/mag_dist_' + key + '.txt'
- dim = (128, 128, 128) # in px (z, y, x)
- res = 10.0 # in nm
- phi = pi/4
+ dim = (1, 1, 5) # in px (z, y, x)
+ res = 1.0 # in nm
+ phi = pi/2
# Geometry parameters:
- center = (64, 64, 64) # in px (z, y, x), index starts with 0!
- width = (1, 50, 50) # in px (z, y, x)
- radius = 25 # in px
- height = 1 # in px
- pos = (0, 63) # in px (tuple of length 2)
- pixel = (0, 63, 63) # in px (z, y, x), index starts with 0!
+ center = (dim[0]/2-0.5, dim[1]/2-0.5, dim[2]/2-0.5) # in px (z, y, x), index starts with 0!
+ width = (dim[0]/2, dim[1]/2, dim[2]/2) # in px (z, y, x)
+ radius = dim[2]/4 # in px
+ height = dim[0]/2 # in px
+ pos = (0, dim[1]/2) # in px (tuple of length 2)
+ pixel = (0, 0, 1) # in px (z, y, x), index starts with 0!
# Determine the magnetic shape:
if key == 'slab':
mag_shape = mc.Shapes.slab(dim, center, width)
diff --git a/scripts/create distributions/create_vortex.py b/scripts/create distributions/create_vortex.py
index e96ea13..794a02b 100644
--- a/scripts/create distributions/create_vortex.py
+++ b/scripts/create distributions/create_vortex.py
@@ -32,10 +32,10 @@ def create_vortex():
filename = directory + '/mag_dist_vortex.txt'
res = 10.0 # in nm
density = 1
- dim = (1, 128, 128)
- center = (0, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
+ dim = (64, 64, 64)
+ center = (int(dim[0]/2)-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
radius = 0.25 * dim[1]
- height = 1
+ height = dim[0]/4
# Create magnetic shape:
mag_shape = mc.Shapes.disc(dim, center, radius, height)
mag_data = MagData(res, mc.create_mag_dist_vortex(mag_shape))
diff --git a/scripts/kernel.py b/scripts/kernel.py
new file mode 100644
index 0000000..6fc4bd0
--- /dev/null
+++ b/scripts/kernel.py
@@ -0,0 +1,52 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Sep 27 14:35:46 2013
+
+@author: Jan
+"""
+
+import numpy as np
+from numpy import pi
+
+import matplotlib.pyplot as plt
+from matplotlib.ticker import IndexLocator
+
+N = 16
+T = 1
+fs = 1. / T
+fn = fs / 2.
+
+x = np.linspace(-N/2, N/2, N, endpoint=False)
+y = np.abs(x) / (x + 1E-30)**3
+y[N/2] = 0
+Y_ft = np.fft.fftshift(np.fft.rfft(y))
+
+f = np.linspace(0, fn, N/2+1)
+Y = -4 * pi * 1j / (f + 1E-30)
+Y[0] = 0
+y_ft = np.fft.fftshift(np.fft.irfft(Y))
+
+#np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))[:, 0]**2
+
+
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+axis.axhline(0, color='k')
+axis.plot(x, y, label='Real space method')
+axis.plot(x, y_ft, label='Real space method')
+axis.grid()
+axis.legend()
+axis.set_title('Real Space Kernel')
+axis.set_xlim(-N/2, N/2)
+axis.xaxis.set_major_locator(IndexLocator(base=N/8, offset=0))
+
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+axis.axhline(0, color='k')
+#axis.plot(f, np.abs(Y)**2, label='Real space method')
+axis.plot(f, np.abs(Y_ft)**2, label='Real space method')
+axis.grid()
+axis.legend()
+axis.set_title('Fourier Space Kernel')
+#axis.set_xlim(0, dim[1]-1)
+#axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
\ No newline at end of file
diff --git a/scripts/paper 1/ch5-0-evaluation_and_comparison.py b/scripts/paper 1/ch5-0-evaluation_and_comparison.py
index 7c12cb5..1f53da2 100644
--- a/scripts/paper 1/ch5-0-evaluation_and_comparison.py
+++ b/scripts/paper 1/ch5-0-evaluation_and_comparison.py
@@ -11,6 +11,7 @@ import traceback
import pdb
import os
+import numpy as np
from numpy import pi
import shelve
@@ -20,6 +21,8 @@ from pyramid.magdata import MagData
import matplotlib.pyplot as plt
+from matplotlib.ticker import FixedFormatter, IndexLocator, LinearLocator
+
def run():
@@ -30,6 +33,23 @@ def run():
os.makedirs(directory)
data_shelve = shelve.open(directory + '/paper_1_shelve')
+ ###############################################################################################
+ print 'CH5-0 KERNELS'
+
+ x = np.linspace(-5, 5, 1000)
+
+ y_r = x/np.abs(x)**3
+ y_k = x/np.abs(x)**2
+ fig = plt.figure()
+ axis = fig.add_subplot(1, 1, 1, aspect='equal')
+ axis.plot(x,y_r, 'r', label=r'$r/|r|^3$')
+ axis.plot(x,y_k, 'b', label=r'$k/|k|^2$')
+ axis.set_xlim(-5, 5)
+ axis.set_ylim(-5, 5)
+ axis.axvline(0, linewidth=2, color='k')
+ axis.axhline(0, linewidth=2, color='k')
+ axis.legend()
+
###############################################################################################
print 'CH5-0 MAGNETIC DISTRIBUTIONS'
@@ -40,9 +60,9 @@ def run():
else:
print '--CREATE MAGNETIC DISTRIBUTIONS'
# Input parameters:
- res = 0.5 # in nm
+ res = 1.0 # in nm
phi = pi/2
- dim = (32, 256, 256) # in px (z, y, x)
+ dim = (16, 128, 128) # in px (z, y, x)
# Create magnetic shape:
center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts at 0!
radius = dim[1]/4 # in px
@@ -50,10 +70,10 @@ def run():
mag_shape = mc.Shapes.disc(dim, center, radius, height)
print '--CREATE MAGN. DISTR. OF HOMOG. MAG. DISC'
mag_data_disc = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
- mag_data_disc.scale_down(3)
+ mag_data_disc.scale_down(2)
print '--CREATE MAGN. DISTR. OF VORTEX STATE DISC'
mag_data_vort = MagData(res, mc.create_mag_dist_vortex(mag_shape, center))
- mag_data_vort.scale_down(3)
+ mag_data_vort.scale_down(2)
# Mayavi-Plots:
mag_data_disc.quiver_plot3d()
mag_data_vort.quiver_plot3d()
@@ -66,9 +86,21 @@ def run():
# Plot MagData (Disc):
mag_data_disc.quiver_plot('Homog. magn. disc', axis=axes[0])
axes[0].set_aspect('equal')
+ axes[0].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
+ axes[0].yaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
+ axes[0].xaxis.set_major_formatter(FixedFormatter([16*i for i in range(10)]))
+ axes[0].yaxis.set_major_formatter(FixedFormatter([16*i for i in range(10)]))
+ axes[0].set_xlabel('x [nm]', fontsize=15)
+ axes[0].set_ylabel('x [ym]', fontsize=15)
# Plot MagData (Disc):
mag_data_vort.quiver_plot('Vortex state disc', axis=axes[1])
axes[1].set_aspect('equal')
+ axes[1].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
+ axes[1].yaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
+ axes[1].xaxis.set_major_formatter(FixedFormatter([16*i for i in range(10)]))
+ axes[1].yaxis.set_major_formatter(FixedFormatter([16*i for i in range(10)]))
+ axes[1].set_xlabel('x [nm]', fontsize=15)
+ axes[1].set_ylabel('x [ym]', fontsize=15)
# Save Plots:
plt.figtext(0.15, 0.15, 'a)', fontsize=30)
plt.figtext(0.57, 0.15, 'b)', fontsize=30)
diff --git a/scripts/paper 1/ch5-4-comparison_of_methods_new.py b/scripts/paper 1/ch5-4-comparison_of_methods_new.py
new file mode 100644
index 0000000..b99f764
--- /dev/null
+++ b/scripts/paper 1/ch5-4-comparison_of_methods_new.py
@@ -0,0 +1,306 @@
+#! python
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jul 26 14:37:30 2013
+
+@author: Jan
+"""
+
+
+import time
+import pdb
+import traceback
+import sys
+import os
+
+import numpy as np
+from numpy import pi
+
+import matplotlib.pyplot as plt
+from matplotlib.ticker import NullLocator, LogLocator, LogFormatter
+
+import pyramid.magcreator as mc
+import pyramid.projector as pj
+import pyramid.phasemapper as pm
+import pyramid.analytic as an
+from pyramid.magdata import MagData
+import shelve
+
+
+def run():
+
+ print '\nACCESS SHELVE'
+ # Create / Open databank:
+ directory = '../../output/paper 1'
+ if not os.path.exists(directory):
+ os.makedirs(directory)
+ data_shelve = shelve.open(directory + '/paper_1_shelve')
+
+ ###############################################################################################
+ print 'CH5-4 METHOD COMPARISON'
+
+ key = 'ch5-4-method_comparison_new'
+ if key in data_shelve:
+ print '--LOAD METHOD DATA'
+ (data_disc_fourier0, data_vort_fourier0,
+ data_disc_fourier3, data_vort_fourier3,
+ data_disc_fourier10, data_vort_fourier10,
+ data_disc_real_s, data_vort_real_s,
+ data_disc_real_d, data_vort_real_d) = data_shelve[key]
+ else:
+ # Input parameters:
+ steps = 5
+ res = 0.25 # in nm
+ phi = pi/2
+ dim = (64, 512, 512) # in px (z, y, x)
+ center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts at 0!
+ radius = dim[1]/4 # in px
+ height = dim[0]/2 # in px
+
+ print '--CREATE MAGNETIC SHAPE'
+ mag_shape = mc.Shapes.disc(dim, center, radius, height)
+ # Create MagData (4 times the size):
+ print '--CREATE MAG. DIST. HOMOG. MAGN. DISC'
+ mag_data_disc = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
+ print '--CREATE MAG. DIST. VORTEX STATE DISC'
+ mag_data_vort = MagData(res, mc.create_mag_dist_vortex(mag_shape, center))
+
+ # Create Data Arrays:
+ dim_list = [dim[2]/2**i for i in np.linspace(1, steps, steps)]
+ data_disc_fourier0 = np.vstack((dim_list, np.zeros((2, steps))))
+ data_vort_fourier0 = np.vstack((dim_list, np.zeros((2, steps))))
+ data_disc_fourier3 = np.vstack((dim_list, np.zeros((2, steps))))
+ data_vort_fourier3 = np.vstack((dim_list, np.zeros((2, steps))))
+ data_disc_fourier10 = np.vstack((dim_list, np.zeros((2, steps))))
+ data_vort_fourier10 = np.vstack((dim_list, np.zeros((2, steps))))
+ data_disc_real_s = np.vstack((dim_list, np.zeros((2, steps))))
+ data_vort_real_s = np.vstack((dim_list, np.zeros((2, steps))))
+ data_disc_real_d = np.vstack((dim_list, np.zeros((2, steps))))
+ data_vort_real_d = np.vstack((dim_list, np.zeros((2, steps))))
+
+ for i in range(steps):
+ # Scale mag_data, resolution and dimensions:
+ mag_data_disc.scale_down()
+ mag_data_vort.scale_down()
+ dim = mag_data_disc.dim
+ res = mag_data_disc.res
+ center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px, index starts at 0!
+ radius = dim[1]/4 # in px
+ height = dim[0]/2 # in px
+
+ print '--res =', res, 'nm', 'dim =', dim
+
+ print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
+ # Create projections along z-axis:
+ projection_disc = pj.simple_axis_projection(mag_data_disc)
+ # Analytic solution:
+ phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
+ # Fourier unpadded:
+ padding = 0
+ start_time = time.clock()
+ phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+ data_disc_fourier0[2, i] = time.clock() - start_time
+ print '------time (disc, fourier0) =', data_disc_fourier0[2, i]
+ phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
+ data_disc_fourier0[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+ # Fourier padding 3:
+ padding = 3
+ start_time = time.clock()
+ phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+ data_disc_fourier3[2, i] = time.clock() - start_time
+ print '------time (disc, fourier3) =', data_disc_fourier3[2, i]
+ phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
+ data_disc_fourier3[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+ # Fourier padding 10:
+ padding = 10
+ start_time = time.clock()
+ phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+ data_disc_fourier10[2, i] = time.clock() - start_time
+ print '------time (disc, fourier10) =', data_disc_fourier10[2, i]
+ phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
+ data_disc_fourier10[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+ # Real space slab:
+ u_kern_f = pm.get_kernel_fourier('slab', 'u', np.shape(projection_disc[0]), res)
+ v_kern_f = pm.get_kernel_fourier('slab', 'v', np.shape(projection_disc[0]), res)
+ start_time = time.clock()
+ phase_num_disc = pm.phase_mag_real_fast(res, projection_disc, (v_kern_f, u_kern_f))
+ data_disc_real_s[2, i] = time.clock() - start_time
+ print '------time (disc, real slab) =', data_disc_real_s[2, i]
+ phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
+ data_disc_real_s[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+ # Real space disc:
+ u_kern_f = pm.get_kernel_fourier('disc', 'u', np.shape(projection_disc[0]), res)
+ v_kern_f = pm.get_kernel_fourier('disc', 'v', np.shape(projection_disc[0]), res)
+ start_time = time.clock()
+ phase_num_disc = pm.phase_mag_real_fast(res, projection_disc, (v_kern_f, u_kern_f))
+ data_disc_real_d[2, i] = time.clock() - start_time
+ print '------time (disc, real disc) =', data_disc_real_d[2, i]
+ phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
+ data_disc_real_d[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
+
+ print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
+ # Create projections along z-axis:
+ projection_vort = pj.simple_axis_projection(mag_data_vort)
+ # Analytic solution:
+ phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
+ # Fourier unpadded:
+ padding = 0
+ start_time = time.clock()
+ phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+ data_vort_fourier0[2, i] = time.clock() - start_time
+ print '------time (vortex, fourier0) =', data_vort_fourier0[2, i]
+ phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
+ data_vort_fourier0[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+ # Fourier padding 3:
+ padding = 3
+ start_time = time.clock()
+ phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+ data_vort_fourier3[2, i] = time.clock() - start_time
+ print '------time (vortex, fourier3) =', data_vort_fourier3[2, i]
+ phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
+ data_vort_fourier3[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+ # Fourier padding 10:
+ padding = 10
+ start_time = time.clock()
+ phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+ data_vort_fourier10[2, i] = time.clock() - start_time
+ print '------time (vortex, fourier10) =', data_vort_fourier10[2, i]
+ phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
+ data_vort_fourier10[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+ # Real space slab:
+ u_kern_f = pm.get_kernel_fourier('slab', 'u', np.shape(projection_vort[0]), res)
+ v_kern_f = pm.get_kernel_fourier('slab', 'v', np.shape(projection_vort[0]), res)
+ start_time = time.clock()
+ phase_num_vort = pm.phase_mag_real_fast(res, projection_vort, (v_kern_f, u_kern_f))
+ data_vort_real_s[2, i] = time.clock() - start_time
+ print '------time (vortex, real slab) =', data_vort_real_s[2, i]
+ phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
+ data_vort_real_s[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+ # Real space disc:
+ u_kern_f = pm.get_kernel_fourier('disc', 'u', np.shape(projection_vort[0]), res)
+ v_kern_f = pm.get_kernel_fourier('disc', 'v', np.shape(projection_vort[0]), res)
+ start_time = time.clock()
+ phase_num_vort = pm.phase_mag_real_fast(res, projection_vort, (v_kern_f, u_kern_f))
+ data_vort_real_d[2, i] = time.clock() - start_time
+ print '------time (vortex, real disc) =', data_vort_real_d[2, i]
+ phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
+ data_vort_real_d[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
+
+ print '--SHELVE METHOD DATA'
+ data_shelve[key] = (data_disc_fourier0, data_vort_fourier0,
+ data_disc_fourier3, data_vort_fourier3,
+ data_disc_fourier10, data_vort_fourier10,
+ data_disc_real_s, data_vort_real_s,
+ data_disc_real_d, data_vort_real_d)
+
+ print '--PLOT/SAVE METHOD DATA'
+
+ # Plot using shared rows and colums:
+ fig, axes = plt.subplots(2, 2, sharex='col', sharey='row', figsize=(12, 8))
+ fig.tight_layout(rect=(0.05, 0.05, 0.95, 0.95))
+ fig.suptitle('Method Comparison', fontsize=20)
+
+ # Plot duration against res (disc) [top/left]:
+ axes[1, 0].set_xscale('log')
+ axes[1, 0].set_yscale('log')
+ axes[1, 0].plot(data_disc_fourier0[0], data_disc_fourier0[1], ':bs')
+ axes[1, 0].plot(data_disc_fourier3[0], data_disc_fourier3[1], ':bo')
+ axes[1, 0].plot(data_disc_fourier10[0], data_disc_fourier10[1], ':b^')
+ axes[1, 0].plot(data_disc_real_s[0], data_disc_real_s[1], '--rs')
+ axes[1, 0].plot(data_disc_real_d[0], data_disc_real_d[1], '--ro')
+ axes[1, 0].set_xlabel('grid size [px]', fontsize=15)
+ axes[1, 0].set_ylabel('RMS [mrad]', fontsize=15)
+ axes[1, 0].set_xlim(12, 350)
+ axes[1, 0].tick_params(axis='both', which='major', labelsize=14)
+ axes[1, 0].xaxis.set_major_locator(LogLocator(base=2))
+ axes[1, 0].xaxis.set_major_formatter(LogFormatter(base=2))
+ axes[1, 0].xaxis.set_minor_locator(NullLocator())
+ axes[1, 0].grid()
+
+ # Plot RMS against res (disc) [bottom/left]:
+ plt.tick_params(axis='both', which='major', labelsize=14)
+ axes[0, 0].set_xscale('log')
+ axes[0, 0].set_yscale('log')
+ axes[0, 0].plot(data_disc_fourier0[0], data_disc_fourier0[2], ':bs')
+ axes[0, 0].plot(data_disc_fourier3[0], data_disc_fourier3[2], ':bo')
+ axes[0, 0].plot(data_disc_fourier10[0], data_disc_fourier10[2], ':b^')
+ axes[0, 0].plot(data_disc_real_s[0], data_disc_real_s[2], '--rs')
+ axes[0, 0].plot(data_disc_real_d[0], data_disc_real_d[2], '--ro')
+ axes[0, 0].set_title('Homog. magn. disc', fontsize=18)
+ axes[0, 0].set_ylabel('duration [s]', fontsize=15)
+ axes[0, 0].set_xlim(12, 350)
+ axes[0, 0].tick_params(axis='both', which='major', labelsize=14)
+ axes[0, 0].xaxis.set_major_locator(LogLocator(base=2))
+ axes[0, 0].xaxis.set_major_formatter(LogFormatter(base=2))
+ axes[0, 0].xaxis.set_minor_locator(NullLocator())
+ axes[0, 0].grid()
+
+ # Plot duration against res (vortex) [top/right]:
+ axes[1, 1].set_xscale('log')
+ axes[1, 1].set_yscale('log')
+ axes[1, 1].plot(data_vort_fourier0[0], data_vort_fourier0[1], ':bs',
+ label='Fourier padding=0')
+ axes[1, 1].plot(data_vort_fourier3[0], data_vort_fourier3[1], ':bo',
+ label='Fourier padding=3')
+ axes[1, 1].plot(data_vort_fourier10[0], data_vort_fourier10[1], ':b^',
+ label='Fourier padding=10')
+ axes[1, 1].plot(data_vort_real_s[0], data_vort_real_s[1], '--rs',
+ label='Real space (slab)')
+ axes[1, 1].plot(data_vort_real_d[0], data_vort_real_d[1], '--ro',
+ label='Real space (disc)')
+ axes[1, 1].set_xlabel('grid size [px]', fontsize=15)
+ axes[1, 1].set_xlim(12, 350)
+ axes[1, 1].tick_params(axis='both', which='major', labelsize=14)
+ axes[1, 1].xaxis.set_major_locator(LogLocator(base=2))
+ axes[1, 1].xaxis.set_major_formatter(LogFormatter(base=2))
+ axes[1, 1].xaxis.set_minor_locator(NullLocator())
+ axes[1, 1].grid()
+
+ # Plot RMS against res (vortex) [bottom/right]:
+ axes[0, 1].set_xscale('log')
+ axes[0, 1].set_yscale('log')
+ axes[0, 1].plot(data_vort_fourier0[0], data_vort_fourier0[2], ':bs',
+ label='Fourier padding=0')
+ axes[0, 1].plot(data_vort_fourier3[0], data_vort_fourier3[2], ':bo',
+ label='Fourier padding=3')
+ axes[0, 1].plot(data_vort_fourier10[0], data_vort_fourier10[2], ':b^',
+ label='Fourier padding=10')
+ axes[0, 1].plot(data_vort_real_s[0], data_vort_real_s[2], '--rs',
+ label='Real space (slab)')
+ axes[0, 1].plot(data_vort_real_d[0], data_vort_real_d[2], '--ro',
+ label='Real space (disc)')
+ axes[0, 1].set_title('Vortex state disc', fontsize=18)
+ axes[0, 1].set_xlim(12, 350)
+ axes[0, 1].tick_params(axis='both', which='major', labelsize=14)
+ axes[0, 1].xaxis.set_major_locator(LogLocator(base=2))
+ axes[0, 1].xaxis.set_major_formatter(LogFormatter(base=2))
+ axes[0, 1].xaxis.set_minor_locator(NullLocator())
+ axes[0, 1].grid()
+
+ # Add legend:
+ axes[1, 1].legend(bbox_to_anchor=(0, 0, 0.955, 0.615), bbox_transform=fig.transFigure,
+ prop={'size':12})
+
+ # Save figure as .png:
+ plt.show()
+ plt.figtext(0.12, 0.85, 'a)', fontsize=30)
+ plt.figtext(0.57, 0.85, 'b)', fontsize=30)
+ plt.figtext(0.12, 0.15, 'c)', fontsize=30)
+ plt.figtext(0.57, 0.15, 'd)', fontsize=30)
+ plt.savefig(directory + '/ch5-3-method comparison.png', bbox_inches='tight')
+
+ ###############################################################################################
+ print 'CLOSING SHELVE\n'
+ # Close shelve:
+ data_shelve.close()
+
+ ###############################################################################################
+
+
+if __name__ == "__main__":
+ try:
+ run()
+ except:
+ type, value, tb = sys.exc_info()
+ traceback.print_exc()
+ pdb.post_mortem(tb)
diff --git a/scripts/reconstruction/reconstruct_core_shell.py b/scripts/reconstruction/reconstruct_core_shell.py
new file mode 100644
index 0000000..1cb4c37
--- /dev/null
+++ b/scripts/reconstruction/reconstruct_core_shell.py
@@ -0,0 +1,94 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Nov 13 10:25:53 2013
+
+@author: Jan
+"""
+
+import sys
+import traceback
+import pdb
+import os
+
+import numpy as np
+from numpy import pi
+
+import shelve
+
+import pyramid.phasemapper as pm
+import pyramid.projector as pj
+import pyramid.holoimage as hi
+from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
+
+import matplotlib.pyplot as plt
+from matplotlib.animation import ArtistAnimation
+
+from matplotlib.ticker import FixedFormatter, IndexLocator, LinearLocator
+
+
+print '\nACCESS SHELVE'
+# Create / Open databank:
+directory = '../../output'
+if not os.path.exists(directory):
+ os.makedirs(directory)
+data_shelve = shelve.open(directory + '/reconstruction/reconstruct_shelf')
+
+key = 'reconstruct_core_shell'
+if key in data_shelve:
+ (projections, phase_maps, mag_data) = data_shelve[key]
+else:
+ mag_data = MagData.load_from_llg(directory +
+ '/magnetic distributions/mag_dist_core_shell_disc.txt')
+ mag_data.quiver_plot()
+ res = mag_data.res
+ dim = mag_data.dim
+ slices = 30
+ angles = np.linspace(0, 2*pi, slices, endpoint=True)
+ projections = []
+ phase_maps = []
+ for tilt_angle in angles:
+ print '{:.2f} pi'.format(tilt_angle/pi)
+ # Project along tilt_angle
+ projection = pj.single_tilt_projection(mag_data, tilt_angle)
+ projections.append(projection)
+ # Construct phase maps:
+ phase_map = PhaseMap(res, pm.phase_mag_real_conv(res, projection), unit='rad')
+ phase_maps.append(phase_map)
+ # Display the combinated plots with phasemap and holography image:
+ title = 'Phase at {:.2f} pi'.format(tilt_angle/pi)
+ # phase_map.display(title, limit=None)
+ # hi.display_combined(phase_map, 1, title=title)
+ # plt.savefig(directory+'/reconstruction/tilt_series_comb_'+
+ # '{:.2f}'.format(tilt_angle/pi)+'pi.png')
+ # Save data
+ data_shelve[key] = (projections, phase_maps, mag_data)
+
+# Close shelve:
+data_shelve.close()
+
+# Plot Animation:
+fig = plt.figure()
+limit = None
+images_ph = []
+for i in range(len(phase_maps)):
+ images_ph.append([plt.imshow(phase_maps[i].phase, cmap='RdBu')])
+ani_ph = ArtistAnimation(fig, images_ph, interval=50, blit=True)
+plt.show()
+fig_phase = plt.figure()
+#fig_projection = plt.figure()
+images_pr = []
+for i in range(len(phase_maps)):
+ images_pr.append([plt.quiver(projections[i][1], projections[i][0], pivot='middle',
+ angles='xy', scale_units='xy', headwidth=6, headlength=7)])
+ani_pr = ArtistAnimation(fig_phase, images_pr, interval=50, blit=True)
+plt.show()
+
+# RECONSTRUCTION:
+mask = mag_data.get_mask()
+dim = mag_data.dim
+res = mag_data.res
+mag_data_reconstruct = MagData(res, (np.zeros(dim),)*3)
+print np.shape(mag_data_reconstruct.magnitude)
+
+mag_vector = mag_data_reconstruct.get_vector(mask)
diff --git a/scripts/robert_filter.py b/scripts/robert_filter.py
new file mode 100644
index 0000000..8c0423d
--- /dev/null
+++ b/scripts/robert_filter.py
@@ -0,0 +1,58 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Oct 29 15:05:39 2013
+
+@author: Jan
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+from numpy import pi
+
+w = 8
+l = 108
+x = np.linspace(0, 20*l, 20*l)
+
+
+# Kastenfunktion:
+y = np.where(x <= l, 1, 0)
+
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+axis.plot(x, y)
+axis.set_xlim(0, 2*l)
+axis.set_ylim(0, 1.2)
+axis.set_title('Kastenfunction --> FFT --> Sinc')
+
+y_fft = np.fft.fft(y)
+
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+axis.plot(x, np.abs(y_fft), 'k', label='abs')
+axis.plot(x, np.real(y_fft), 'r', label='real')
+axis.plot(x, np.imag(y_fft), 'b', label='imag')
+axis.set_xlim(0, 200)
+axis.legend()
+axis.set_title('Kastenfunction --> FFT --> Sinc')
+
+
+# Kastenfunktion mit Cos-"Softness":
+y = np.where(x <= l, 1, np.where(x >= l+w, 0, 0.5*(1+np.cos(pi*(x-l)/w))))
+
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+axis.plot(x, y)
+axis.set_xlim(0, 2*l)
+axis.set_ylim(0, 1.2)
+axis.set_title('Kastenfunction mit Cos-Kanten')
+
+y_fft = np.fft.fft(y)
+
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+axis.plot(x, np.abs(y_fft), 'k', label='abs')
+axis.plot(x, np.real(y_fft), 'r', label='real')
+axis.plot(x, np.imag(y_fft), 'b', label='imag')
+axis.set_xlim(0, 200)
+axis.legend()
+axis.set_title('Kastenfunction mit Cos-Kanten')
\ No newline at end of file
diff --git a/scripts/test.py b/scripts/test.py
new file mode 100644
index 0000000..cfbfb91
--- /dev/null
+++ b/scripts/test.py
@@ -0,0 +1,27 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Nov 05 14:52:48 2013
+
+@author: Jan
+"""
+
+import numpy as np
+from numpy import pi
+import matplotlib.pyplot as plt
+from matplotlib.ticker import MaxNLocator, FuncFormatter
+
+x = np.linspace(0, 2*pi, 1000)
+y = np.where(x<pi, -1/np.tan(x+1E-30), 1/np.tan(x+1E-30))#20*np.arctan(np.exp(x))
+
+def func(x, pos):
+ return '{}'.format(x*10)
+
+formatter = FuncFormatter(lambda x, pos: '{:g}'.format(x))
+
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+axis.plot(x, y)
+axis.set_ylim(-10, 10)
+axis.yaxis.set_major_formatter(formatter)
+
+plt.show()
\ No newline at end of file
diff --git a/scripts/test_arb_projection.py b/scripts/test_arb_projection.py
index a802772..fa71062 100644
--- a/scripts/test_arb_projection.py
+++ b/scripts/test_arb_projection.py
@@ -14,7 +14,7 @@ import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator, NullLocator
-dim = (4, 4)
+dim = (8, 8)
offset = (dim[0]/2., dim[1]/2.)
phi = 1.000000001*pi/4
field = np.zeros(dim)
diff --git a/scripts/test_methods.py b/scripts/test_methods.py
index ec3cdf3..be89108 100644
--- a/scripts/test_methods.py
+++ b/scripts/test_methods.py
@@ -27,10 +27,12 @@ def compare_methods():
# Input parameters:
res = 10.0 # in nm
phi = 0
- density = 1
+ theta = pi/2
+ tilt = pi/4
+ density = 0.25
dim = (128, 128, 128) # in px (z, y, x)
# Create magnetic shape:
- geometry = 'slab'
+ geometry = 'sphere'
if geometry == 'slab':
center = (dim[0]/2-0.5, dim[1]/2-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts at 0!
width = (dim[0]/2, dim[1]/2., dim[2]/2.) # in px (z, y, x)
@@ -45,17 +47,18 @@ def compare_methods():
radius = dim[0]/4 # in px
mag_shape = mc.Shapes.sphere(dim, center, radius)
# Project the magnetization data:
- mag_data = MagData(res, mc.create_mag_dist_homog(mag_shape, phi, theta=0))
+ mag_data = MagData(res, mc.create_mag_dist_homog(mag_shape, phi, theta))
# density = 0.3
# mag_data = MagData.load_from_llg('../output/magnetic distributions/mag_dist_sphere.txt')
res = mag_data.res
- mag_data.quiver_plot()
+ mag_data.quiver_plot(ax_slice=dim[1]/2)
# mag_data.quiver_plot3d()
import time
start = time.time()
- projection = pj.single_tilt_projection(mag_data, tilt=pi/4)
+ projection = pj.single_tilt_projection(mag_data, tilt)
+ pj.quiver_plot(projection)
print 'Total projection time:', time.time() - start
# Construct phase maps:
phase_map_mag = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding=1))
@@ -66,7 +69,15 @@ def compare_methods():
phase_map = PhaseMap(res, phase_map_mag.phase+phase_map_elec.phase)
hi.display_combined(phase_map, density)
-
+
+ import matplotlib.pyplot as plt
+ x = range(dim[2])
+ y = phase_map_elec.phase[dim[1]/2, :]
+ fig = plt.figure()
+ axis = fig.add_subplot(1, 1, 1)
+ axis.plot(x, y, label='Real space method')
+ axis.grid()
+
if __name__ == "__main__":
try:
diff --git a/scripts/test_vadim.py b/scripts/test_vadim.py
new file mode 100644
index 0000000..55d6e8d
--- /dev/null
+++ b/scripts/test_vadim.py
@@ -0,0 +1,186 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Nov 07 16:47:52 2013
+
+@author: Jan
+"""
+
+import numpy as np
+from numpy import pi
+
+import os
+
+import pyramid.magcreator as mc
+import pyramid.projector as pj
+import pyramid.phasemapper as pm
+import pyramid.holoimage as hi
+from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
+
+import matplotlib.pyplot as plt
+
+Q_E = 1#1.602176565e-19
+EPSILON_0 = 1#8.8e-9
+C_E = 1
+
+
+def calculate_charge_batch(phase_map):
+ def calculate_charge(grad_y, grad_x, l, r, t, b):
+ left = np.sum(-grad_x[b:t, l])
+ top = np.sum(grad_y[t, l:r])
+ right = np.sum(grad_x[b:t, r])
+ bottom = np.sum(-grad_y[b, l:r])
+ integral = left + top + right + bottom
+ Q = -EPSILON_0/C_E * integral
+ return Q
+ grad_y, grad_x = np.gradient(phase_map.phase, phase_map.res, phase_map.res)
+ xi, yi = np.zeros(t-b), np.zeros(t-b)
+ for i, ti in enumerate(np.arange(b, t)):
+ xi[i] = ti
+ yi[i] = calculate_charge(grad_y, grad_x, l, r, ti, b)
+ return xi, yi
+
+directory = '../output/vadim/'
+if not os.path.exists(directory):
+ os.makedirs(directory)
+
+res = 1.0 # in nm
+phi = pi/4
+density = 30
+dim = (64, 256, 256) # in px (z, y, x)
+
+l = dim[2]/4.
+r = dim[2]/4. + dim[2]/2.
+b = dim[1]/4.
+t = dim[1]/4. + dim[1]/2.
+
+# Create magnetic shape:
+center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts at 0!
+radius = dim[1]/8 # in px
+height = dim[0]/4 # in px
+width = (dim[0]/4, dim[1]/4, dim[2]/4)
+mag_shape_sphere = mc.Shapes.sphere(dim, center, radius)
+mag_shape_disc = mc.Shapes.disc(dim, center, radius, height)
+mag_shape_slab = mc.Shapes.slab(dim, center, (height, radius*2, radius*2))
+# Create magnetization distributions:
+mag_data_sphere = MagData(res, mc.create_mag_dist_homog(mag_shape_sphere, phi))
+mag_data_disc = MagData(res, mc.create_mag_dist_homog(mag_shape_disc, phi))
+mag_data_vortex = MagData(res, mc.create_mag_dist_vortex(mag_shape_disc))
+mag_data_slab = MagData(res, mc.create_mag_dist_homog(mag_shape_slab, phi))
+# Project the magnetization data:
+projection_sphere = pj.simple_axis_projection(mag_data_sphere)
+projection_disc = pj.simple_axis_projection(mag_data_disc)
+projection_vortex = pj.simple_axis_projection(mag_data_vortex)
+projection_slab = pj.simple_axis_projection(mag_data_slab)
+# Magnetic phase map of a homogeneously magnetized disc:
+phase_map_mag_disc = PhaseMap(res, pm.phase_mag_real_conv(res, projection_disc))
+phase_map_mag_disc.save_to_txt(directory+'phase_map_mag_disc.txt')
+axis, _ = hi.display_combined(phase_map_mag_disc, density)
+axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+plt.savefig(directory+'phase_map_mag_disc.png')
+x, y = calculate_charge_batch(phase_map_mag_disc)
+np.savetxt(directory+'charge_integral_mag_disc.txt', np.array([x, y]).T)
+plt.figure()
+plt.plot(x, y)
+plt.savefig(directory+'charge_integral_mag_disc.png')
+# Magnetic phase map of a vortex state disc:
+phase_map_mag_vortex = PhaseMap(res, pm.phase_mag_real_conv(res, projection_vortex))
+phase_map_mag_vortex.save_to_txt(directory+'phase_map_mag_vortex.txt')
+axis, _ = hi.display_combined(phase_map_mag_vortex, density)
+axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+plt.savefig(directory+'phase_map_mag_vortex.png')
+x, y = calculate_charge_batch(phase_map_mag_vortex)
+np.savetxt(directory+'charge_integral_mag_vortex.txt', np.array([x, y]).T)
+plt.figure()
+plt.plot(x, y)
+plt.savefig(directory+'charge_integral_mag_vortex.png')
+# MIP phase of a slab:
+phase_map_mip_slab = PhaseMap(res, pm.phase_elec(res, projection_slab, v_0=1, v_acc=300000))
+phase_map_mip_slab.save_to_txt(directory+'phase_map_mip_slab.txt')
+axis, _ = hi.display_combined(phase_map_mip_slab, density)
+axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+plt.savefig(directory+'phase_map_mip_slab.png')
+x, y = calculate_charge_batch(phase_map_mip_slab)
+np.savetxt(directory+'charge_integral_mip_slab.txt', np.array([x, y]).T)
+plt.figure()
+plt.plot(x, y)
+plt.savefig(directory+'charge_integral_mip_slab.png')
+# MIP phase of a disc:
+phase_map_mip_disc = PhaseMap(res, pm.phase_elec(res, projection_disc, v_0=1, v_acc=300000))
+phase_map_mip_disc.save_to_txt(directory+'phase_map_mip_disc.txt')
+axis, _ = hi.display_combined(phase_map_mip_disc, density)
+axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+plt.savefig(directory+'phase_map_mip_disc.png')
+x, y = calculate_charge_batch(phase_map_mip_disc)
+np.savetxt(directory+'charge_integral_mip_disc.txt', np.array([x, y]).T)
+plt.figure()
+plt.plot(x, y)
+plt.savefig(directory+'charge_integral_mip_disc.png')
+# MIP phase of a sphere:
+phase_map_mip_sphere = PhaseMap(res, pm.phase_elec(res, projection_sphere, v_0=1, v_acc=300000))
+phase_map_mip_sphere.save_to_txt(directory+'phase_map_mip_sphere.txt')
+axis, _ = hi.display_combined(phase_map_mip_sphere, density)
+axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+plt.savefig(directory+'phase_map_mip_sphere.png')
+x, y = calculate_charge_batch(phase_map_mip_sphere)
+np.savetxt(directory+'charge_integral_mip_sphere.txt', np.array([x, y]).T)
+plt.figure()
+plt.plot(x, y)
+plt.savefig(directory+'charge_integral_mip_sphere.png')
+
+
+
+
+
+# Display phasemap and holography image:
+
+#y, x = np.mgrid[0:dim[1], 0:dim[2]]
+
+#phase_charge = 1/((xx-center[2])**2+(yy-center[1])**2)
+#phase_map.phase = phase_charge
+
+#q_e = 1#1.602176565e-19
+#epsilon0 = 1#8.8e-9
+#C_e = 1
+#
+#dist = np.sqrt((y-center[1])**2+(x-center[2])**2)
+#V = q_e/((4*pi*epsilon0)*dist)
+#thickness = 10000
+#phase = C_e*q_e/(4*pi*epsilon0) * 1/2*np.arcsinh(thickness/dist)
+
+#phase_map = PhaseMap(res, phase)#C_e * res * V.sum(axis=0)
+
+
+## USER INPUT:
+#phase_map.display()
+#point1, point2 = plt.ginput(n=2, timeout=0)
+#plt.close()
+#l = np.round(min(point1[0], point2[0]))
+#r = np.round(max(point1[0], point2[0]))
+#b = np.round(min(point1[1], point2[1]))
+#t = np.round(max(point1[1], point2[1]))
--
GitLab