diff --git a/docs/pyramid.rst b/docs/pyramid.rst
index b4ade82b76e0c3727715973488714e50c216339e..3491e66ec3ab49baf2a06d77b1b654c0f9670ee2 100644
--- a/docs/pyramid.rst
+++ b/docs/pyramid.rst
@@ -19,6 +19,15 @@ pyramid Package
:show-inheritance:
:special-members:
+:mod:`kernel` Module
+-----------------------
+
+.. automodule:: pyramid.kernel
+ :members:
+ :undoc-members:
+ :show-inheritance:
+ :special-members:
+
:mod:`magcreator` Module
------------------------
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index 63782bf9bfa87bde56ac6e8719755cb19563eb06..d33e5c92653888980594ae4fa3f488e6269d3082 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -9,6 +9,8 @@ magdata
Class for the storage of magnetizatin data.
projector
Create projections of a given magnetization distribution.
+kernel
+ Class for the kernel matrix representing one magnetized pixel.
phasemapper
Create magnetic and electric phase maps from magnetization data.
phasemap
diff --git a/pyramid/analytic.py b/pyramid/analytic.py
index 6035de75c17f040bb87f45e8a252e0a8b91e239a..5445990a923ede7dccd77d4d1232bd90aeff8c27 100644
--- a/pyramid/analytic.py
+++ b/pyramid/analytic.py
@@ -15,15 +15,15 @@ from numpy import pi
PHI_0 = -2067.83 # magnetic flux in T*nm²
-def phase_mag_slab(dim, res, phi, center, width, b_0=1):
+def phase_mag_slab(dim, a, phi, center, width, b_0=1):
'''Calculate the analytic magnetic phase for a homogeneously magnetized slab.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The grid spacing in nm.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
@@ -57,27 +57,27 @@ def phase_mag_slab(dim, res, phi, center, width, b_0=1):
+ F_0(y-y0+Ly/2, x-x0+Lx/2)))
# Process input parameters:
z_dim, y_dim, x_dim = dim
- y0 = res * (center[1] + 0.5) # y0, x0 define the center of a pixel,
- x0 = res * (center[2] + 0.5) # hence: (cellindex + 0.5) * resolution
- Lz, Ly, Lx = res * width[0], res * width[1], res * width[2]
+ y0 = a * (center[1] + 0.5) # y0, x0 define the center of a pixel,
+ x0 = a * (center[2] + 0.5) # hence: (cellindex + 0.5) * grid spacing
+ Lz, Ly, Lx = a * width[0], a * width[1], a * width[2]
coeff = b_0 / (4*PHI_0)
# Create grid:
- x = np.linspace(res/2, x_dim*res-res/2, num=x_dim)
- y = np.linspace(res/2, y_dim*res-res/2, num=y_dim)
+ x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
+ y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return phi_mag(xx, yy)
-def phase_mag_disc(dim, res, phi, center, radius, height, b_0=1):
+def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
'''Calculate the analytic magnetic phase for a homogeneously magnetized disc.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The grid spacing in nm.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
@@ -105,28 +105,28 @@ def phase_mag_disc(dim, res, phi, center, radius, height, b_0=1):
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
- y0 = res * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
- x0 = res * (center[2] + 0.5) # hence: cellindex + 0.5
- Lz = res * height
- R = res * radius
+ y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
+ x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
+ Lz = a * height
+ R = a * radius
coeff = - pi * b_0 / (2*PHI_0)
# Create grid:
- x = np.linspace(res/2, x_dim*res-res/2, num=x_dim)
- y = np.linspace(res/2, y_dim*res-res/2, num=y_dim)
+ x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
+ y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return phi_mag(xx, yy)
-def phase_mag_sphere(dim, res, phi, center, radius, b_0=1):
+def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
'''Calculate the analytic magnetic phase for a homogeneously magnetized sphere.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The grid spacing in nm.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
@@ -152,27 +152,27 @@ def phase_mag_sphere(dim, res, phi, center, radius, b_0=1):
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
- y0 = res * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
- x0 = res * (center[2] + 0.5) # hence: cellindex + 0.5
- R = res * radius
+ y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
+ x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
+ R = a * radius
coeff = - 2./3. * pi * b_0 / PHI_0
# Create grid:
- x = np.linspace(res / 2, x_dim * res - res / 2, num=x_dim)
- y = np.linspace(res / 2, y_dim * res - res / 2, num=y_dim)
+ x = np.linspace(a / 2, x_dim * a - a / 2, num=x_dim)
+ y = np.linspace(a / 2, y_dim * a - a / 2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return phi_mag(xx, yy)
-def phase_mag_vortex(dim, res, center, radius, height, b_0=1):
+def phase_mag_vortex(dim, a, center, radius, height, b_0=1):
'''Calculate the analytic magnetic phase for a vortex state disc.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The grid spacing in nm.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`, which is also the vortex center.
radius : float
@@ -196,14 +196,14 @@ def phase_mag_vortex(dim, res, center, radius, height, b_0=1):
return result
# Process input parameters:
z_dim, y_dim, x_dim = dim
- y0 = res * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
- x0 = res * (center[2] + 0.5) # hence: cellindex + 0.5
- Lz = res * height
- R = res * radius
+ y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
+ x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
+ Lz = a * height
+ R = a * radius
coeff = pi * b_0 * Lz / PHI_0
# Create grid:
- x = np.linspace(res/2, x_dim*res-res/2, num=x_dim)
- y = np.linspace(res/2, y_dim*res-res/2, num=y_dim)
+ x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
+ y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
xx, yy = np.meshgrid(x, y)
# Return phase:
return phi_mag(xx, yy)
diff --git a/pyramid/holoimage.py b/pyramid/holoimage.py
index 639c90d2b8608399cbb4cc038c29b3ee2108a827..f48aa61327458ddbc47a9d57b75dc4d637dd9f2f 100644
--- a/pyramid/holoimage.py
+++ b/pyramid/holoimage.py
@@ -64,7 +64,7 @@ def holo_image(phase_map, density=1):
# Calculate the holography image intensity:
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_map.phase, phase_map.res, phase_map.res)
+ phase_grad_y, phase_grad_x = np.gradient(phase_map.phase, phase_map.a, phase_map.a)
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:
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
new file mode 100644
index 0000000000000000000000000000000000000000..51bf156067655ca2d4ccef6badb7d550167b1a62
--- /dev/null
+++ b/pyramid/kernel.py
@@ -0,0 +1,206 @@
+# -*- coding: utf-8 -*-
+"""Class for the calculation and storage of kernel.
+
+This module provides the :class:`~.Kernel` class whose instances can be used to calculate and
+store the kernel matrix representing the phase of a single pixel for the convolution used in the
+phase calculation. The phasemap of a single pixel for two orthogonal directions (`u` and `v`) are
+stored seperately as 2-dimensional matrices. The Jacobi matrix of the phasemapping just depends
+on the kernel and can be calculated via the :func:`~.get_jacobi` function. Storing the Jacobi
+matrix uses much memory, thus it is also possible to directly get the multiplication of a given
+vector with the (transposed) Jacobi matrix without explicit calculation of the latter.
+It is possible to load data from and save them to NetCDF4 files. See :class:`~.Kernel` for further
+information.
+
+"""
+
+
+import numpy as np
+
+
+PHI_0 = -2067.83 # magnetic flux in T*nm²
+
+
+class Kernel:
+ '''Class for calculating kernel matrices for the phase calculation.
+
+ Represents the phase of a single magnetized pixel for two orthogonal directions (`u` and `v`),
+ which can be accessed via the corresponding attributes. The default elementary geometry is
+ `disc`, but can also be specified as the phase of a `slab` representation of a single
+ magnetized pixel. During the construction, a few attributes are calculated that are used in
+ the convolution during phase calculation.
+
+ Attributes
+ ----------
+ dim : tuple (N=2)
+ Dimensions of the projected magnetization grid.
+ a : float
+ The grid spacing in nm.
+ geometry : {'disc', 'slab'}, optional
+ The elementary geometry of the single magnetized pixel.
+ b_0 : float, optional
+ The saturation magnetic induction. Default is 1.
+ u : :class:`~numpy.ndarray` (N=3)
+ The phase contribution of one pixel magnetized in u-direction.
+ v : :class:`~numpy.ndarray` (N=3)
+ The phase contribution of one pixel magnetized in v-direction.
+ u_fft : :class:`~numpy.ndarray` (N=3)
+ The real FFT of the phase contribution of one pixel magnetized in u-direction.
+ v_fft : :class:`~numpy.ndarray` (N=3)
+ The real FFT of the phase contribution of one pixel magnetized in v-direction.
+ dim_fft : tuple (N=2)
+ Dimensions of the grid, which is used for the FFT. Calculated by adding the dimensions
+ `dim` of the magnetization grid and the dimensions of the kernel (given by ``2*dim-1``)
+ and increasing to the next multiple of 2 (for faster FFT).
+ slice_fft : tuple (N=2) of :class:`slice`
+ A tuple of :class:`slice` objects to extract the original field of view from the increased
+ size (size_fft) of the grid for the FFT-convolution.
+
+ '''
+
+ def __init__(self, dim, a, b_0=1, geometry='disc'):
+ '''Constructor for a :class:`~.Kernel` object for representing a kernel matrix.
+
+ Parameters
+ ----------
+ dim : tuple (N=2)
+ Dimensions of the projected magnetization grid.
+ a : float
+ The grid spacing in nm.
+ b_0 : float, optional
+ The saturation magnetic induction. Default is 1.
+ geometry : {'disc', 'slab'}, optional
+ The elementary geometry of the single magnetized pixel.
+
+ '''
+ # Function for the phase of an elementary geometry:
+ def get_elementary_phase(geometry, n, m, res):
+ if geometry == 'disc':
+ in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
+ return m / (n**2 + m**2 + 1E-30) * in_or_out
+ elif geometry == '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 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))
+ # Set basic properties:
+ self.dim = dim # !!! size of the FOV, not the kernel (kernel is bigger)!
+ self.a = a
+ self.geometry = geometry
+ self.b_0 = b_0
+ # Calculate kernel (single pixel phase):
+ coeff = -a**2 / (2*PHI_0)
+ 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)
+ self.u = coeff * get_elementary_phase(geometry, uu, vv, a)
+ self.v = coeff * get_elementary_phase(geometry, vv, uu, a)
+ # Calculate Fourier trafo of kernel components:
+ dim_combined = 3*np.array(dim) - 1 # dim + (2*dim - 1) magnetisation + kernel
+ self.dim_fft = 2 ** np.ceil(np.log2(dim_combined)).astype(int) # next multiple of 2
+ self.slice_fft = (slice(dim[0]-1, 2*dim[0]-1), slice(dim[1]-1, 2*dim[1]-1))
+ self.u_fft = np.fft.rfftn(self.u, self.dim_fft)
+ self.v_fft = np.fft.rfftn(self.v, self.dim_fft)
+
+ def get_jacobi(self):
+ '''Calculate the Jacobi matrix for the phase calculation from a projected magnetization.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ jacobi : :class:`~numpy.ndarray` (N=2)
+ Jacobi matrix containing the derivatives of the phase at every pixel with respect to
+ the projected magetization. Has `N` columns for the `u`-component of the magnetization
+ and `N` columns for the `v`-component (from left to right) and ``N**2`` rows for the
+ phase at every pixel.
+
+ Notes
+ -----
+ Just use for small dimensions, Jacobi Matrix scales with order of ``N**4``.
+
+ '''
+ v_dim, u_dim = self.dim
+ jacobi = np.zeros((v_dim*u_dim, 2*v_dim*u_dim))
+# nc.get_jacobi_core(dim[0], dim[1], v_phi, u_phi, jacobi)
+# return jacobi
+ for j in range(v_dim):
+ for i in range(u_dim):
+ u_column = i + u_dim*j
+ v_column = i + u_dim*j + u_dim*v_dim
+ u_min = (u_dim-1) - i
+ u_max = (2*u_dim-1) - i
+ v_min = (v_dim-1) - j
+ v_max = (2*v_dim-1) - j
+ # u_dim*v_dim columns for the u-component:
+ jacobi[:, u_column] = self.u[v_min:v_max, u_min:u_max].reshape(-1)
+ # u_dim*v_dim columns for the v-component (note the minus!):
+ jacobi[:, v_column] = -self.v[v_min:v_max, u_min:u_max].reshape(-1)
+ return jacobi
+
+ def multiply_jacobi(self, vector):
+ '''Calculate the product of the Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
+ (row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
+ ``N**2`` elements to the `v`-component of the magnetization.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
+
+ '''
+ v_dim, u_dim = self.dim
+ size = v_dim * u_dim
+ assert len(vector) == 2*size, 'vector size not compatible!'
+ result = np.zeros(size)
+ for s in range(size): # column-wise (two columns at a time, u- and v-component)
+ i = s % u_dim
+ j = int(s/u_dim)
+ u_min = (u_dim-1) - i
+ u_max = (2*u_dim-1) - i
+ v_min = (v_dim-1) - j
+ v_max = (2*v_dim-1) - j
+ result += vector[s]*self.u[v_min:v_max, u_min:u_max].reshape(-1) # u
+ result += vector[s+size]*-self.v[v_min:v_max, u_min:u_max].reshape(-1) # v
+ return result
+
+ def multiply_jacobi_T(self, vector):
+ '''Calculate the product of the transposed Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
+ (row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
+ ``N**2`` elements to the `v`-component of the magnetization.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Product of the transposed Jacobi matrix (which is not explicitely calculated) with
+ the vector.
+
+ '''
+ v_dim, u_dim = self.dim
+ size = v_dim * u_dim
+ assert len(vector) == size, 'vector size not compatible!'
+ result = np.zeros(2*size)
+ for s in range(size): # row-wise (two rows at a time, u- and v-component)
+ i = s % u_dim
+ j = int(s/u_dim)
+ u_min = (u_dim-1) - i
+ u_max = (2*u_dim-1) - i
+ v_min = (v_dim-1) - j
+ v_max = (2*v_dim-1) - j
+ result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1)) # u
+ result[s+size] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1)) # v
+ return result
diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py
index 11eea245d188c37f7cb195dc55deb2353aa42c2d..4d8f6b4340a1c61e90ccb9abf1c51f0dc74bf354 100644
--- a/pyramid/magcreator.py
+++ b/pyramid/magcreator.py
@@ -22,7 +22,7 @@ class Shapes:
'''Class containing functions for generating magnetic shapes.
- The :class:`~.Shapes` class is a collection of some mehtods that return a 3-dimensional
+ The :class:`~.Shapes` class is a collection of some methods that return a 3-dimensional
matrix that represents the magnetized volume and consists of values between 0 and 1.
This matrix is used in the functions of the :mod:`~.magcreator` module to create
:class:`~pyramid.magdata.MagData` objects which store the magnetic informations.
@@ -61,7 +61,7 @@ class Shapes:
@classmethod
def disc(cls, dim, center, radius, height, axis='z'):
- '''Create the shape of a zylindrical disc in x-, y-, or z-direction.
+ '''Create the shape of a cylindrical disc in x-, y-, or z-direction.
Parameters
----------
@@ -139,7 +139,7 @@ class Shapes:
@classmethod
def ellipsoid(cls, dim, center, width):
- '''Create the shape of a sphere.
+ '''Create the shape of an ellipsoid.
Parameters
----------
@@ -147,8 +147,8 @@ class Shapes:
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the sphere in pixel coordinates `(z, y, x)`.
- radius : float
- The radius of the sphere in pixel coordinates.
+ width : tuple (N=3)
+ The width of the ellipsoid `(z, y, x)`.
Returns
-------
@@ -158,7 +158,7 @@ class Shapes:
'''
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
-# assert radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!' # TODO: change
+ assert np.shape(width) == (3,), 'Parameter width has to be a a tuple of length 3!'
mag_shape = np.array([[[np.sqrt((x-center[2])**2/(width[2]/2)**2
+ (y-center[1])**2/(width[1]/2)**2
+ (z-center[0])**2/(width[0]/2)**2) <= 1
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 792ca026a4b1d12a43bd97dd3634763733da1b84..6ac5834f774c99c573da69c2e8c3dffdbc92a28a 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -2,9 +2,9 @@
"""Class for the storage of magnetizatin data.
This module provides the :class:`~.MagData` class whose instances can be used to store
-magnetization data for 3 components for a 3-dimensional grid. It is possible to load data from
-NetCDF4 or LLG (.txt) files or to save the data in these formats. Also plotting methods are
-provided. See :class:`~.MagData` for further information.
+magnetization distributions with 3 components for a 3-dimensional grid. It is possible to load
+data from NetCDF4 or LLG (.txt) files or to save the data in these formats. Plotting methods
+are also provided. See :class:`~.MagData` for further information.
"""
@@ -29,12 +29,13 @@ class MagData:
This is useful, if the `magnitude` is more complex and more than one magnetized object should
be represented by the :class:`~.MagData` object, which can be added one after another by the
:func:`~.add_magnitude` function. The dimensions `dim` of the grid will be set as soon as the
- magnitude is specified. However, `res` has to be always specified at construction time.
+ magnitude is specified. However, the grid spacding `a` has to be always specified at
+ construction time.
Attributes
----------
- res : float
- The resolution of the grid (grid spacing in nm)
+ a : float
+ The grid spacing in nm.
dim : tuple (N=3)
Dimensions of the grid.
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
@@ -43,22 +44,17 @@ class MagData:
'''
- def __init__(self, res, magnitude=None):
+ def __init__(self, a, magnitude=None):
'''Constructor for a :class:`~.MagData` object for storing magnetization data.
Parameters
----------
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The grid spacing in nm.
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3), optional
The `z`-, `y`- and `x`-component of the magnetization vector for every 3D-gridpoint
as a tuple. Is zero everywhere if not specified.
- Returns
- -------
- mag_data : :class:`~.MagData`
- The 3D magnetic distribution as a :class:`~.MagData` object.
-
'''
if magnitude is not None:
dim = np.shape(magnitude[0])
@@ -70,7 +66,7 @@ class MagData:
else:
self.magnitude = None
self.dim = None
- self.res = res
+ self.a = a
def add_magnitude(self, magnitude):
'''Add a given magnitude to the magnitude of the :class:`~.MagData`.
@@ -175,7 +171,7 @@ class MagData:
Notes
-----
- Acts in place and changes dimensions and resolution accordingly.
+ Acts in place and changes dimensions and grid spacing accordingly.
Only possible, if each axis length is a power of 2!
'''
@@ -188,43 +184,20 @@ class MagData:
y_mag = self.magnitude[1].reshape(self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2)
x_mag = self.magnitude[2].reshape(self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2)
self.dim = (self.dim[0]/2, self.dim[1]/2, self.dim[2]/2)
- self.res = self.res * 2
+ self.a = self.a * 2
self.magnitude = (z_mag.mean(axis=5).mean(axis=3).mean(axis=1),
y_mag.mean(axis=5).mean(axis=3).mean(axis=1),
x_mag.mean(axis=5).mean(axis=3).mean(axis=1))
- @classmethod
- def load_from_llg(cls, filename):
- '''Construct :class:`~.MagData` object from LLG-file.
-
- Parameters
- ----------
- filename : string
- The name of the LLG-file from which to load the data.
-
- Returns
- -------
- mag_data: :class:`~.MagData`
- A :class:`~.MagData` object containing the loaded data.
- '''
- SCALE = 1.0E-9 / 1.0E-2 # From cm to nm
- data = np.genfromtxt(filename, skip_header=2)
- x_dim, y_dim, z_dim = np.genfromtxt(filename, dtype=int, skip_header=1,
- skip_footer=len(data[:, 0]))
- res = (data[1, 0] - data[0, 0]) / SCALE
- # Reshape in Python and Igor is different, Python fills rows first, Igor columns:
- x_mag, y_mag, z_mag = [data[:, i].reshape(z_dim, y_dim, x_dim) for i in range(3, 6)]
- return MagData(res, (z_mag, y_mag, x_mag))
-
- def save_to_llg(self, filename='magdata_output.txt'):
+ def save_to_llg(self, filename='..\output\magdata_output.txt'):
'''Save magnetization data in a file with LLG-format.
Parameters
----------
filename : string, optional
The name of the LLG-file in which to store the magnetization data.
- The default is 'magdata_output.txt'.
+ The default is '..\output\magdata_output.txt'.
Returns
-------
@@ -232,11 +205,11 @@ class MagData:
'''
dim = self.dim
- res = self.res * 1.0E-9 / 1.0E-2 # from nm to cm
+ a = self.a * 1.0E-9 / 1.0E-2 # from nm to cm
# Create 3D meshgrid and reshape it and the magnetization into a list where x varies first:
- zz, yy, xx = np.mgrid[res/2:(dim[0]*res-res/2):dim[0]*1j,
- res/2:(dim[1]*res-res/2):dim[1]*1j,
- res/2:(dim[2]*res-res/2):dim[2]*1j]
+ zz, yy, xx = np.mgrid[a/2:(dim[0]*a-a/2):dim[0]*1j,
+ a/2:(dim[1]*a-a/2):dim[1]*1j,
+ a/2:(dim[2]*a-a/2):dim[2]*1j]
xx = xx.reshape(-1)
yy = yy.reshape(-1)
zz = zz.reshape(-1)
@@ -252,13 +225,13 @@ class MagData:
for cell in row) for row in data))
@classmethod
- def load_from_netcdf4(cls, filename):
- '''Construct :class:`~.DataMag` object from NetCDF4-file.
+ def load_from_llg(cls, filename):
+ '''Construct :class:`~.MagData` object from LLG-file.
Parameters
----------
filename : string
- The name of the NetCDF4-file from which to load the data. Standard format is '\*.nc'.
+ The name of the LLG-file from which to load the data.
Returns
-------
@@ -266,13 +239,14 @@ class MagData:
A :class:`~.MagData` object containing the loaded data.
'''
- mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
- res = mag_file.res
- z_mag = mag_file.variables['z_mag'][:]
- y_mag = mag_file.variables['y_mag'][:]
- x_mag = mag_file.variables['x_mag'][:]
- mag_file.close()
- return MagData(res, (z_mag, y_mag, x_mag))
+ SCALE = 1.0E-9 / 1.0E-2 # From cm to nm
+ data = np.genfromtxt(filename, skip_header=2)
+ x_dim, y_dim, z_dim = np.genfromtxt(filename, dtype=int, skip_header=1,
+ skip_footer=len(data[:, 0]))
+ a = (data[1, 0] - data[0, 0]) / SCALE
+ # Reshape in Python and Igor is different, Python fills rows first, Igor columns:
+ x_mag, y_mag, z_mag = [data[:, i].reshape(z_dim, y_dim, x_dim) for i in range(3, 6)]
+ return MagData(a, (z_mag, y_mag, x_mag))
def save_to_netcdf4(self, filename='..\output\magdata_output.nc'):
'''Save magnetization data in a file with NetCDF4-format.
@@ -281,7 +255,7 @@ class MagData:
----------
filename : string, optional
The name of the NetCDF4-file in which to store the magnetization data.
- The default is 'magdata_output.nc'.
+ The default is '..\output\magdata_output.nc'.
Returns
-------
@@ -289,7 +263,7 @@ class MagData:
'''
mag_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
- mag_file.res = self.res
+ mag_file.a = self.a
mag_file.createDimension('z_dim', self.dim[0])
mag_file.createDimension('y_dim', self.dim[1])
mag_file.createDimension('x_dim', self.dim[2])
@@ -299,9 +273,31 @@ class MagData:
z_mag[:] = self.magnitude[0]
y_mag[:] = self.magnitude[1]
x_mag[:] = self.magnitude[2]
- print mag_file
mag_file.close()
+ @classmethod
+ def load_from_netcdf4(cls, filename):
+ '''Construct :class:`~.DataMag` object from NetCDF4-file.
+
+ Parameters
+ ----------
+ filename : string
+ The name of the NetCDF4-file from which to load the data. Standard format is '\*.nc'.
+
+ Returns
+ -------
+ mag_data: :class:`~.MagData`
+ A :class:`~.MagData` object containing the loaded data.
+
+ '''
+ mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
+ a = mag_file.a
+ z_mag = mag_file.variables['z_mag'][:]
+ y_mag = mag_file.variables['y_mag'][:]
+ x_mag = mag_file.variables['x_mag'][:]
+ mag_file.close()
+ return MagData(a, (z_mag, y_mag, x_mag))
+
def quiver_plot(self, title='Magnetic Distribution', filename=None, axis=None,
proj_axis='z', ax_slice=None):
'''Plot a slice of the magnetization as a quiver plot.
@@ -365,6 +361,7 @@ class MagData:
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
plt.show()
+ return axis
def quiver_plot3d(self):
'''Plot the magnetization as 3D-vectors in a quiverplot.
@@ -378,12 +375,12 @@ class MagData:
None
'''
- res = self.res
+ a = self.a
dim = self.dim
# Create points and vector components as lists:
- zz, yy, xx = np.mgrid[res/2:(dim[0]*res-res/2):dim[0]*1j,
- res/2:(dim[1]*res-res/2):dim[1]*1j,
- res/2:(dim[2]*res-res/2):dim[2]*1j]
+ zz, yy, xx = np.mgrid[a/2:(dim[0]*a-a/2):dim[0]*1j,
+ a/2:(dim[1]*a-a/2):dim[1]*1j,
+ a/2:(dim[2]*a-a/2):dim[2]*1j]
xx = xx.reshape(-1)
yy = yy.reshape(-1)
zz = zz.reshape(-1)
diff --git a/pyramid/numcore/phase_mag_real.pyx b/pyramid/numcore/phase_mag_real.pyx
index d67993367b216a854df217f5d9fa186f391979c1..d3a23568466b482c51e04087ef92196bd865a93a 100644
--- a/pyramid/numcore/phase_mag_real.pyx
+++ b/pyramid/numcore/phase_mag_real.pyx
@@ -66,8 +66,22 @@ def get_jacobi_core(
unsigned int v_dim, unsigned int u_dim,
double[:, :] v_phi, double[:, :] u_phi,
double[:, :] jacobi):
- '''DOCSTRING!'''
- # TODO: Docstring!!!
+ '''Numerical core routine for the jacobi matrix calculation.
+
+ Parameters
+ ----------
+ v_dim, u_dim : int
+ Dimensions of the projection along the two major axes.
+ v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
+ Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
+ jacobi : :class:`~numpy.ndarray` (N=2)
+ Jacobi matrix which is filled by this routine.
+
+ Returns
+ -------
+ None
+
+ '''
cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
for j in range(v_dim):
for i in range(u_dim):
@@ -86,22 +100,3 @@ def get_jacobi_core(
jacobi[row, u_column] = u_phi[q, p]
# v-component (note the minus!):
jacobi[row, v_column] = -v_phi[q, p]
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def get_jacobi_core(
- unsigned int v_dim, unsigned int u_dim,
- np.ndarray v_phi, np.ndarray u_phi,
- np.ndarray jacobi):
- '''DOCSTRING!'''
- # TODO: Docstring!!!
- cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
- for j in range(v_dim):
- for i in range(u_dim):
- # v_dim*u_dim columns for the u-component:
- jacobi[:, i+u_dim*j] = \
- np.reshape(u_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)
- # v_dim*u_dim columns for the v-component (note the minus!):
- jacobi[:, u_dim*v_dim+i+u_dim*j] = \
- -np.reshape(v_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index 90623640ff7c23e5824d111b289a5a98f3b6b493..2fd0e669d39eb3af758878deb04609769932d800 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -2,8 +2,8 @@
"""Class for the storage of phase data.
This module provides the :class:`~.PhaseMap` class whose instances can be used to store
-phase data for a 2-dimensional grid. It is possible to load data from NetCDF4 or LLG (.txt) files
-or to save the data in these formats. Also plotting methods are provided. See :class:`~.PhaseMap`
+phase data for a 2-dimensional grid. It is possible to load data from NetCDF4 or textfiles or to
+save the data in these formats. Also plotting methods are provided. See :class:`~.PhaseMap`
for further information.
"""
@@ -23,13 +23,13 @@ class PhaseMap:
'''Class for storing phase map data.
Represents 2-dimensional phase maps. the phase information itself is stored in `phase`.
- The dimensions `dim` of the grid with resolution `res` will be calculated at construction
- time, but `res` has to be specified.
+ The dimensions `dim` of the grid with grid spacing `a` will be calculated at construction
+ time, but `a` has to be specified.
Attributes
----------
- res : float
- The resolution of the grid (grid spacing in nm)
+ a : float
+ The grid spacing in nm.
dim : tuple (N=2)
Dimensions of the grid.
phase : :class:`~numpy.ndarray` (N=2)
@@ -45,13 +45,13 @@ class PhaseMap:
'mrad': 1E3,
'µrad': 1E6}
- def __init__(self, res, phase, unit='rad'):
+ def __init__(self, a, phase, unit='rad'):
'''Constructor for a :class:`~.PhaseMap` object for storing phase data.
Parameters
----------
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The grid spacing in nm.
phase : :class:`~numpy.ndarray` (N=2)
Matrix containing the phase shift.
unit : {'rad', 'mrad', 'µrad'}, optional
@@ -59,15 +59,10 @@ class PhaseMap:
because the phase is scaled accordingly. Does not change the phase itself, which is
always in `rad`.
- Returns
- -------
- phase_map : :class:`~.PhaseMap`
- The 2D phase shift as a :class:`~.PhaseMap` object.
-
'''
dim = np.shape(phase)
assert len(dim) == 2, 'Phasemap has to be 2-dimensional!'
- self.res = res
+ self.a = a
self.dim = dim
self.unit = unit
self.phase = phase
@@ -90,27 +85,6 @@ class PhaseMap:
assert unit in ['rad', 'mrad']
self.unit = unit
- @classmethod
- def load_from_txt(cls, filename):
- '''Construct :class:`~.PhaseMap` object from a human readable txt-file.
-
- Parameters
- ----------
- filename : string
- The name of the file from which to load the data.
-
- Returns
- -------
- phase_map : :class:`~.PhaseMap`
- A :class:`~.PhaseMap` object containing the loaded data.
-
- '''
- with open(filename, 'r') as phase_file:
- phase_file.readline() # Headerline is not used
- res = float(phase_file.readline()[13:-4])
- phase = np.loadtxt(filename, delimiter='\t', skiprows=2)
- return PhaseMap(res, phase)
-
def save_to_txt(self, filename='..\output\phasemap_output.txt'):
'''Save :class:`~.PhaseMap` data in a file with txt-format.
@@ -118,7 +92,7 @@ class PhaseMap:
----------
filename : string
The name of the file in which to store the phase map data.
- The default is 'phasemap_output.txt'.
+ The default is '..\output\phasemap_output.txt'.
Returns
-------
@@ -127,29 +101,29 @@ class PhaseMap:
'''
with open(filename, 'w') as phase_file:
phase_file.write('{}\n'.format(filename.replace('.txt', '')))
- phase_file.write('resolution = {} nm\n'.format(self.res))
+ phase_file.write('grid spacing = {} nm\n'.format(self.a))
np.savetxt(phase_file, self.phase, fmt='%7.6e', delimiter='\t')
@classmethod
- def load_from_netcdf4(cls, filename):
- '''Construct :class:`~.PhaseMap` object from NetCDF4-file.
+ def load_from_txt(cls, filename):
+ '''Construct :class:`~.PhaseMap` object from a human readable txt-file.
Parameters
----------
filename : string
- The name of the NetCDF4-file from which to load the data. Standard format is '\*.nc'.
+ The name of the file from which to load the data.
Returns
-------
- phase_map: :class:`~.PhaseMap`
+ phase_map : :class:`~.PhaseMap`
A :class:`~.PhaseMap` object containing the loaded data.
'''
- phase_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
- res = phase_file.res
- phase = phase_file.variables['phase'][:]
- phase_file.close()
- return PhaseMap(res, phase)
+ with open(filename, 'r') as phase_file:
+ phase_file.readline() # Headerline is not used
+ a = float(phase_file.readline()[15:-4])
+ phase = np.loadtxt(filename, delimiter='\t', skiprows=2)
+ return PhaseMap(a, phase)
def save_to_netcdf4(self, filename='..\output\phasemap_output.nc'):
'''Save :class:`~.PhaseMap` data in a file with NetCDF4-format.
@@ -158,7 +132,7 @@ class PhaseMap:
----------
filename : string, optional
The name of the NetCDF4-file in which to store the phase data.
- The default is 'phasemap_output.nc'.
+ The default is '..\output\phasemap_output.nc'.
Returns
-------
@@ -166,14 +140,34 @@ class PhaseMap:
'''
phase_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
- phase_file.res = self.res
+ phase_file.a = self.a
phase_file.createDimension('v_dim', self.dim[0])
phase_file.createDimension('u_dim', self.dim[1])
phase = phase_file.createVariable('phase', 'f', ('v_dim', 'u_dim'))
phase[:] = self.phase
- print phase_file
phase_file.close()
+ @classmethod
+ def load_from_netcdf4(cls, filename):
+ '''Construct :class:`~.PhaseMap` object from NetCDF4-file.
+
+ Parameters
+ ----------
+ filename : string
+ The name of the NetCDF4-file from which to load the data. Standard format is '\*.nc'.
+
+ Returns
+ -------
+ phase_map: :class:`~.PhaseMap`
+ A :class:`~.PhaseMap` object containing the loaded data.
+
+ '''
+ phase_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
+ a = phase_file.a
+ phase = phase_file.variables['phase'][:]
+ phase_file.close()
+ return PhaseMap(a, phase)
+
def display(self, title='Phase Map', cmap='RdBu', limit=None, norm=None, axis=None):
'''Display the phasemap as a colormesh.
@@ -214,8 +208,8 @@ class PhaseMap:
axis.set_ylim(0, np.shape(phase)[0])
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.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
+ axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlabel('x [nm]', fontsize=15)
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 8ddf90be580397e29c19657c99326174ec0abb7b..d69f262a2deaf273e855435044153fcc5e6833d9 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -4,9 +4,7 @@
This module executes several forward models to calculate the magnetic or electric phase map from
a given projection of a 3-dimensional magnetic distribution (see :mod:`~pyramid.projector`).
For the magnetic phase map, an approach using real space and one using Fourier space is provided.
-The real space approach also calculates the Jacobi matrix, which is used for the reconstruction in
-the :mod:`~pyramid.reconstructor` module. The electrostatic contribution is calculated by using
-the assumption of a mean inner potentail (MIP).
+The electrostatic contribution is calculated by using the assumption of a mean inner potential.
"""
@@ -15,8 +13,7 @@ import numpy as np
from numpy import pi
import pyramid.numcore as nc
-
-from scipy.signal import fftconvolve
+from pyramid.kernel import Kernel
PHI_0 = -2067.83 # magnetic flux in T*nm²
@@ -26,117 +23,22 @@ Q_E = 1.602E-19 # electron charge in C
C = 2.998E8 # speed of light in m/s
-# TODO: Kernel class? Creator, transform to fourier, get_jacobi?
-
-class Kernel:
- # TODO: Docstrings!!!
-
- def __init__(self, dim, res, geometry='disc', b_0=1):
- # TODO: Docstring!!!
- def get_elementary_phase(geometry, n, m, res):
- if geometry == '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 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 geometry == 'disc':
- in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
- return m / (n**2 + m**2 + 1E-30) * in_or_out
- self.dim = dim # !!! size of the FOV, kernel is bigger!
- self.res = res
- self.geometry = geometry
- self.b_0 = b_0
- coeff = -res**2 / (2*PHI_0)
- 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)
- self.u = coeff * get_elementary_phase(geometry, uu, vv, res)
- self.v = coeff * get_elementary_phase(geometry, vv, uu, res)
- size = 3*np.array(dim) - 1 # dim + (2*dim - 1) magnetisation + kernel
- fsize = 2 ** np.ceil(np.log2(size)).astype(int) # next multiple of 2
- self.u_fft = np.fft.rfftn(self.u, fsize)
- self.v_fft = np.fft.rfftn(self.v, fsize)
-
- def get_jacobi(self):
- # TODO: Docstring!!!
- '''CAUTIOUS! Just use for small dimensions or computer will explode violently!'''
- v_dim, u_dim = self.dim
- jacobi = np.zeros((v_dim*u_dim, 2*v_dim*u_dim))
- # nc.get_jacobi_core(dim[0], dim[1], v_phi, u_phi, jacobi)
- # return jacobi
- for j in range(v_dim):
- for i in range(u_dim):
- u_column = i + u_dim*j
- v_column = i + u_dim*j + u_dim*v_dim
- u_min = (u_dim-1) - i
- u_max = (2*u_dim-1) - i
- v_min = (v_dim-1) - j
- v_max = (2*v_dim-1) - j
- # u_dim*v_dim columns for the u-component:
- jacobi[:, u_column] = self.u[v_min:v_max, u_min:u_max].reshape(-1)
- # u_dim*v_dim columns for the v-component (note the minus!):
- jacobi[:, v_column] = -self.v[v_min:v_max, u_min:u_max].reshape(-1)
- return jacobi
-
- def multiply_jacobi(self, vector):
- # TODO: Docstring!!!
- # vector: v_dim*u_dim elements for u_mag and v_dim*u_dim elements for v_mag
- v_dim, u_dim = self.dim
- size = v_dim * u_dim
- assert len(vector) == 2*size, 'vector size not compatible!'
- result = np.zeros(size)
- for s in range(size): # column-wise (two columns at a time, u- and v-component)
- i = s % u_dim
- j = int(s/u_dim)
- u_min = (u_dim-1) - i
- u_max = (2*u_dim-1) - i
- v_min = (v_dim-1) - j
- v_max = (2*v_dim-1) - j
- result += vector[s]*self.u[v_min:v_max, u_min:u_max].reshape(-1) # u
- result += vector[s+size]*-self.v[v_min:v_max, u_min:u_max].reshape(-1) # v
- return result
-
- def multiply_jacobi_T(self, vector):
- # TODO: Docstring!!!
- # vector: v_dim*u_dim elements for u_mag and v_dim*u_dim elements for v_mag
- v_dim, u_dim = self.dim
- size = v_dim * u_dim
- assert len(vector) == size, 'vector size not compatible!'
- result = np.zeros(2*size)
- for s in range(size): # row-wise (two rows at a time, u- and v-component)
- i = s % u_dim
- j = int(s/u_dim)
- u_min = (u_dim-1) - i
- u_max = (2*u_dim-1) - i
- v_min = (v_dim-1) - j
- v_max = (2*v_dim-1) - j
- result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1)) # u
- result[s+size] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1)) # v
- return result
-
-
-def phase_mag_real(res, projection, geometry='disc', b_0=1, jacobi=None):
- '''Calculate the magnetic phase from magnetization data (real space approach).
+def phase_mag(a, projection, b_0=1, kernel=None):
+ '''Calculate the magnetic phase from magnetization data.
Parameters
----------
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The 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.
- geometry : {'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.
+ kernel : :class:`~pyramid.kernel.Kernel`, optional
+ Specifies the kernel for the convolution with the magnetization data. If none is specified,
+ one will be created with `disc` as the default geometry.
Returns
-------
@@ -144,86 +46,65 @@ def phase_mag_real(res, projection, geometry='disc', b_0=1, jacobi=None):
The phase as a 2-dimensional array.
'''
- # Process input parameters:
- dim = np.shape(projection[0])
+ # Process input parameters:
v_mag, u_mag = projection[:-1]
+ dim = np.shape(u_mag)
- # Create kernel (lookup-tables for the phase of one pixel):
- kernel = Kernel(dim, res, geometry, b_0)
- u_phi = kernel.u
- v_phi = kernel.v
+ # Create kernel (lookup-tables for the phase of one pixel) if none is given:
+ if kernel is None:
+ kernel = Kernel(dim, a, b_0)
- # Calculation of the phase:
- phase = np.zeros(dim)
- threshold = 0
- if jacobi is not None: # With Jacobian matrix (slower)
- jacobi[:] = 0 # Jacobi matrix --> zeros
- 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[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(dim[0], dim[1], v_phi, u_phi, v_mag, u_mag, phase, threshold)
- # Return the phase:
- return phase
+ # Fourier transform the projected magnetisation:
+ u_mag_fft = np.fft.rfftn(u_mag, kernel.dim_fft)
+ v_mag_fft = np.fft.rfftn(v_mag, kernel.dim_fft)
+ # Convolve the magnetization with the kernel in Fourier space:
+ u_phase = np.fft.irfftn(u_mag_fft * kernel.u_fft, kernel.dim_fft)[kernel.slice_fft].copy()
+ v_phase = np.fft.irfftn(v_mag_fft * kernel.v_fft, kernel.dim_fft)[kernel.slice_fft].copy()
-def phase_mag_real_conv(res, projection, method='disc', b_0=1):
- '''Calculate the magnetic phase from magnetization data (real space approach).
+ # Return the result:
+ return u_phase - v_phase
+
+
+def phase_elec(a, projection, v_0=1, v_acc=30000):
+ '''Calculate the electric phase from magnetization distributions.
Parameters
----------
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The 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
+ 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.
+ v_0 : float, optional
+ The mean inner potential of the specimen in V. The default is 1.
+ v_acc : float, optional
+ The acceleration voltage of the electron microscope in V. The default is 30000.
Returns
-------
phase : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
- ''' # TODO: Docstring!!!
+ '''
# Process input parameters:
- dim = np.shape(projection[0])
- v_mag, u_mag = projection[:-1]
-
- # 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
+ lam = H_BAR / np.sqrt(2 * M_E * Q_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))
+ # return phase:
+ return v_0 * Ce * projection[-1] * a*1E-9
-def phase_mag_real_fast(res, projection, kernels_fourier, b_0=1):
- '''Calculate the magnetic phase from magnetization data (real space approach).
+def phase_mag_real(a, projection, b_0=1, geometry='disc', jacobi=None):
+ '''Calculate the magnetic phase from magnetization data (pure real space, no FFT-convolution).
Parameters
----------
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The 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
+ geometry : {'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
@@ -238,33 +119,44 @@ def phase_mag_real_fast(res, projection, kernels_fourier, b_0=1):
phase : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
- ''' # TODO: Docstring!!!
- # Process input parameters:
+ '''
+ # Process input parameters:
+ dim = np.shape(projection[0])
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
+ # Create kernel (lookup-tables for the phase of one pixel):
+ kernel = Kernel(dim, a, b_0, geometry)
+ u_phi = kernel.u
+ v_phi = kernel.v
- 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
+ # Calculation of the phase:
+ phase = np.zeros(dim)
+ threshold = 0
+ if jacobi is not None: # With Jacobian matrix (slower)
+ jacobi[:] = 0 # Jacobi matrix --> zeros
+ 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[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(dim[0], dim[1], v_phi, u_phi, v_mag, u_mag, phase, threshold)
+ # Return the phase:
+ return phase
-def phase_mag_fourier(res, projection, padding=0, b_0=1):
+def phase_mag_fourier(a, projection, padding=0, b_0=1):
'''Calculate the magnetic phase from magnetization data (Fourier space approach).
Parameters
----------
- res : float
- The resolution of the grid (grid spacing) in nm.
+ a : float
+ The 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.
@@ -294,86 +186,13 @@ 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 = 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_nyq = 0.5 / a # nyquist frequency
+ f_u = np.linspace(0, f_nyq, u_mag_fft.shape[1])
+ f_v = np.linspace(-f_nyq, f_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) #* (8*(res/2)**3)*np.sinc(res/2*f_uu)*np.sinc(res/2*f_vv) * np.exp(2*pi*1j*())
+ phase_fft = coeff * a * (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(res, projection, kernel, 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 phase_elec(res, projection, v_0=1, v_acc=30000):
- '''Calculate the electric phase from magnetization distributions.
-
- 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.
- v_0 : float, optional
- The mean inner potential of the specimen in V. The default is 1.
- v_acc : float, optional
- The acceleration voltage of the electron microscope in V. The default is 30000.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- '''
- # Process input parameters:
- lam = H_BAR / np.sqrt(2 * M_E * Q_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))
- # return phase:
- return v_0 * Ce * projection[-1] * res*1E-9
diff --git a/pyramid/projection.py b/pyramid/projection.py
new file mode 100644
index 0000000000000000000000000000000000000000..01aa038965abdd8d36ab8ac2f25e8f686bb581ae
--- /dev/null
+++ b/pyramid/projection.py
@@ -0,0 +1,7 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Nov 29 12:07:44 2013
+
+@author: Jan
+"""
+
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 00d56c82b364d6b93e65d4485ebab99f08820de5..9fa94d265d7414b04044f63dbc3098597da6a9c2 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -3,8 +3,8 @@
This module creates 2-dimensional projections from 3-dimensional magnetic distributions, which
are stored in :class:`~pyramid.magdata.MagData` objects. Either simple projections along the
-major axes are possible (:func:`~.simple_axis_projection`), or projections along arbitrary
-directions with the use of transfer functions (work in progress).
+major axes are possible (:func:`~.simple_axis_projection`), or projections with a tilt around
+the y-axis. The thickness profile is also calculated and can be used for electric phase maps.
"""
@@ -12,8 +12,6 @@ directions with the use of transfer functions (work in progress).
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
@@ -38,7 +36,7 @@ def simple_axis_projection(mag_data, axis='z', threshold=0):
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. The latter
- has to be multiplied with the resolution for a value in nm.
+ has to be multiplied with the grid spacing for a value in nm.
'''
assert isinstance(mag_data, MagData), 'Parameter mag_data has to be a MagData object!'
@@ -58,16 +56,30 @@ def simple_axis_projection(mag_data, axis='z', threshold=0):
return projection
-def single_tilt_projection(mag_data, tilt=0, threshold=0):
- # TODO: Docstring!!!
- '''Projection with tilt around the (centered!) y-axis'''
+def single_tilt_projection(mag_data, tilt=0):
+ '''
+ Project a magnetization distribution which is tilted around the centered y-axis.
+
+ Parameters
+ ----------
+ mag_data : :class:`~pyramid.magdata.MagData`
+ A :class:`~pyramid.magdata.MagData` object storing the magnetization distribution,
+ which should be projected.
+ tilt : float, optional
+ The counter-clockwise tilt angle around the y-axis. Default is 1.
+
+ Returns
+ -------
+ 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. The latter
+ has to be multiplied with the grid spacing for a value in nm.
+
+ '''
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)!'
def get_position(p, m, b, size):
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):
@@ -75,21 +87,18 @@ def single_tilt_projection(mag_data, tilt=0, threshold=0):
if 0 <= x < size]
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
- return 0
+ lo, up = delta-rho, delta+rho
# Upper boundary:
- if b >= 1:
- w_b = 0.5
+ if up >= 1:
+ w_up = 0.5
else:
- w_b = (b*np.sqrt(1-b**2) + np.arctan(b/np.sqrt(1-b**2))) / pi
+ w_up = (up*np.sqrt(1-up**2) + np.arctan(up/np.sqrt(1-up**2))) / pi
# Lower boundary:
- if a <= -1:
- w_a = -0.5
+ if lo <= -1:
+ w_lo = -0.5
else:
- w_a = (a*np.sqrt(1-a**2) + np.arctan(a/np.sqrt(1-a**2))) / pi
- return w_b - w_a
+ w_lo = (lo*np.sqrt(1-lo**2) + np.arctan(lo/np.sqrt(1-lo**2))) / pi
+ return w_up - w_lo
# Set starting variables:
# length along projection (proj, z), rotation (rot, y) and perpendicular (perp, x) axis:
@@ -101,22 +110,17 @@ def single_tilt_projection(mag_data, tilt=0, threshold=0):
np.zeros((dim_rot, dim_perp)))
# Creating coordinate list of all voxels:
voxels = list(itertools.product(range(dim_proj), range(dim_perp)))
-# print 'voxels:', voxels
# Calculate positions along the projected pixel coordinate system:
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
+ m = np.where(tilt<=pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
b = center[0] - m * center[1]
-# 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
- rho = 0.5 / r # TODO: ratio of radii
- weights = []
+ rho = 0.5 / r
+ weights = {}
for i, voxel in enumerate(voxels):
voxel_weights = []
impacts = get_impact(positions[i], r, dim_perp)
@@ -124,12 +128,11 @@ def single_tilt_projection(mag_data, tilt=0, threshold=0):
distance = np.abs(impact+0.5 - positions[i])
delta = distance / r
voxel_weights.append((impact, get_weight(delta, rho)))
- weights.append((voxel, voxel_weights))
-
-# print weights
+ weights[voxel] = voxel_weights
# Calculate projection with the calculated weights for the voxels (rotation around y-axis):
- for i, (voxel, voxel_weights) in enumerate(weights):
+ for i, voxel in enumerate(weights):
+ voxel_weights = weights[voxel]
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]]
@@ -140,3 +143,199 @@ def single_tilt_projection(mag_data, tilt=0, threshold=0):
projection[2][:, pixel] += weight * mask[voxel[0], :, voxel[1]]
return projection
+
+class Projection:
+
+ '''Class for calculating kernel matrices for the phase calculation.
+
+ Represents the phase of a single magnetized pixel for two orthogonal directions (`u` and `v`),
+ which can be accessed via the corresponding attributes. The default elementary geometry is
+ `disc`, but can also be specified as the phase of a `slab` representation of a single
+ magnetized pixel. During the construction, a few attributes are calculated that are used in
+ the convolution during phase calculation.
+
+ Attributes
+ ----------
+ dim : tuple (N=2)
+ Dimensions of the projected magnetization grid.
+ a : float
+ The grid spacing in nm.
+ geometry : {'disc', 'slab'}, optional
+ The elementary geometry of the single magnetized pixel.
+ b_0 : float, optional
+ The saturation magnetic induction. Default is 1.
+ u : :class:`~numpy.ndarray` (N=3)
+ The phase contribution of one pixel magnetized in u-direction.
+ v : :class:`~numpy.ndarray` (N=3)
+ The phase contribution of one pixel magnetized in v-direction.
+ u_fft : :class:`~numpy.ndarray` (N=3)
+ The real FFT of the phase contribution of one pixel magnetized in u-direction.
+ v_fft : :class:`~numpy.ndarray` (N=3)
+ The real FFT of the phase contribution of one pixel magnetized in v-direction.
+ dim_fft : tuple (N=2)
+ Dimensions of the grid, which is used for the FFT. Calculated by adding the dimensions
+ `dim` of the magnetization grid and the dimensions of the kernel (given by ``2*dim-1``)
+ and increasing to the next multiple of 2 (for faster FFT).
+ slice_fft : tuple (N=2) of :class:`slice`
+ A tuple of :class:`slice` objects to extract the original field of view from the increased
+ size (size_fft) of the grid for the FFT-convolution.
+
+ ''' # TODO: Docstring!
+
+ def __init__(self, u_mag, v_mag, thickness, weights):
+ '''Constructor for a :class:`~.Kernel` object for representing a kernel matrix.
+
+ Parameters
+ ----------
+ dim : tuple (N=2)
+ Dimensions of the projected magnetization grid.
+ a : float
+ The grid spacing in nm.
+ b_0 : float, optional
+ The saturation magnetic induction. Default is 1.
+ geometry : {'disc', 'slab'}, optional
+ The elementary geometry of the single magnetized pixel.
+
+ ''' # TODO: Docstring!
+ self.dim = np.shape(u_mag)
+# self.a = a
+# self.b_0 = b_0
+ self.u = u_mag
+ self.v = v_mag
+ self.weights = weights
+
+ def get_jacobi(self):
+ '''Calculate the Jacobi matrix for the phase calculation from a projected magnetization.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ jacobi : :class:`~numpy.ndarray` (N=2)
+ Jacobi matrix containing the derivatives of the phase at every pixel with respect to
+ the projected magetization. Has `N` columns for the `u`-component of the magnetization
+ and `N` columns for the `v`-component (from left to right) and ``N**2`` rows for the
+ phase at every pixel.
+
+ Notes
+ -----
+ Just use for small dimensions, Jacobi Matrix scales with order of ``N**4``.
+
+ ''' # TODO: Docstring!
+ v_dim, u_dim = self.dim
+ size = v_dim * u_dim
+ jacobi = np.zeros((v_dim*u_dim, 2*v_dim*u_dim))
+ weight_matrix = np.zeros(dim) # Check if correct (probably not, x and z instead of x and y)
+ voxels = list(itertools.product(range(dim_proj), range(dim_perp)))
+ for i, voxel in enumerate(self.weights):
+ voxel_weights = self.weights[voxel]
+ 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:
+ projection[1][:, pixel] += weight * (x_mag[voxel[0], :, voxel[1]]*np.cos(tilt)
+ + z_mag[voxel[0], :, voxel[1]]*np.sin(tilt))
+ # Thickness profile:
+ projection[2][:, pixel] += weight * mask[voxel[0], :, voxel[1]]
+
+ for i, vx in enumerate(self.weights):
+ vx_pos = vx[1]*3 + vx[0]
+ for px, weight in self.weights[voxel]:
+ px_pos = (px[2]*3 + px[1])*3 + px[0]
+ weight_matrix[vx_pos, px_pos] = weight
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ v_dim, u_dim = self.dim
+ jacobi = np.zeros((v_dim*u_dim, 2*v_dim*u_dim))
+# nc.get_jacobi_core(dim[0], dim[1], v_phi, u_phi, jacobi)
+# return jacobi
+ for j in range(v_dim):
+ for i in range(u_dim):
+ u_column = i + u_dim*j
+ v_column = i + u_dim*j + u_dim*v_dim
+ u_min = (u_dim-1) - i
+ u_max = (2*u_dim-1) - i
+ v_min = (v_dim-1) - j
+ v_max = (2*v_dim-1) - j
+ # u_dim*v_dim columns for the u-component:
+ jacobi[:, u_column] = self.u[v_min:v_max, u_min:u_max].reshape(-1)
+ # u_dim*v_dim columns for the v-component (note the minus!):
+ jacobi[:, v_column] = -self.v[v_min:v_max, u_min:u_max].reshape(-1)
+ return jacobi
+
+ def multiply_jacobi(self, vector):
+ '''Calculate the product of the Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
+ (row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
+ ``N**2`` elements to the `v`-component of the magnetization.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
+
+ ''' # TODO: Docstring!
+ v_dim, u_dim = self.dim
+ size = v_dim * u_dim
+ assert len(vector) == 2*size, 'vector size not compatible!'
+ result = np.zeros(size)
+ for s in range(size): # column-wise (two columns at a time, u- and v-component)
+ i = s % u_dim
+ j = int(s/u_dim)
+ u_min = (u_dim-1) - i
+ u_max = (2*u_dim-1) - i
+ v_min = (v_dim-1) - j
+ v_max = (2*v_dim-1) - j
+ result += vector[s]*self.u[v_min:v_max, u_min:u_max].reshape(-1) # u
+ result += vector[s+size]*-self.v[v_min:v_max, u_min:u_max].reshape(-1) # v
+ return result
+
+ def multiply_jacobi_T(self, vector):
+ '''Calculate the product of the transposed Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
+ (row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
+ ``N**2`` elements to the `v`-component of the magnetization.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Product of the transposed Jacobi matrix (which is not explicitely calculated) with
+ the vector.
+
+ ''' # TODO: Docstring!
+ v_dim, u_dim = self.dim
+ size = v_dim * u_dim
+ assert len(vector) == size, 'vector size not compatible!'
+ result = np.zeros(2*size)
+ for s in range(size): # row-wise (two rows at a time, u- and v-component)
+ i = s % u_dim
+ j = int(s/u_dim)
+ u_min = (u_dim-1) - i
+ u_max = (2*u_dim-1) - i
+ v_min = (v_dim-1) - j
+ v_max = (2*v_dim-1) - j
+ result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1)) # u
+ result[s+size] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1)) # v
+ return result
diff --git a/pyramid/reconstructor.py b/pyramid/reconstructor.py
index b31d3e5f2312ccd86fb53cdba909c74d212f362b..7c995dc7d1eb54bc4d5e1ea3755cfa6c4b8fd7cb 100644
--- a/pyramid/reconstructor.py
+++ b/pyramid/reconstructor.py
@@ -49,17 +49,17 @@ def reconstruct_simple_leastsq(phase_map, mask, b_0=1):
'''
# Read in parameters:
y_m = phase_map.phase.reshape(-1) # Measured phase map as a vector
- res = phase_map.res # Resolution
+ a = phase_map.a # Grid spacing
dim = mask.shape # Dimensions of the mag. distr.
count = mask.sum() # Number of pixels with magnetization
lam = 1e-6 # Regularisation parameter
# Create empty MagData object for the reconstruction:
- mag_data_rec = MagData(res, (np.zeros(dim), np.zeros(dim), np.zeros(dim)))
+ mag_data_rec = MagData(a, (np.zeros(dim), np.zeros(dim), np.zeros(dim)))
# Function that returns the phase map for a magnetic configuration x:
def F(x):
mag_data_rec.set_vector(mask, x)
- phase = pm.phase_mag_real(res, pj.simple_axis_projection(mag_data_rec), 'slab', b_0)
+ phase = pm.phase_mag_real(a, pj.simple_axis_projection(mag_data_rec), b_0)
return phase.reshape(-1)
# Cost function which should be minimized:
diff --git a/scripts/compare methods/compare_methods.py b/scripts/compare methods/compare_methods.py
index dd32b50dcfe677f84845008b385478cb1fd2dbb4..c949a8ae9ab21e95c452955341baa483cf547887 100644
--- a/scripts/compare methods/compare_methods.py
+++ b/scripts/compare methods/compare_methods.py
@@ -17,6 +17,7 @@ import pyramid.phasemapper as pm
import pyramid.holoimage as hi
import pyramid.analytic as an
from pyramid.magdata import MagData
+from pyramid.kernel import Kernel
from pyramid.phasemap import PhaseMap
@@ -30,7 +31,7 @@ def compare_methods():
'''
# Input parameters:
b_0 = 1.0 # in T
- res = 1.0 # in nm
+ a = 1.0 # in nm
phi = pi/4
padding = 3
density = 10
@@ -41,56 +42,52 @@ def compare_methods():
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)
mag_shape = mc.Shapes.slab(dim, center, width)
- phase_ana = an.phase_mag_slab(dim, res, phi, center, width, b_0)
+ phase_ana = an.phase_mag_slab(dim, a, phi, center, width, b_0)
elif geometry == 'disc':
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
mag_shape = mc.Shapes.disc(dim, center, radius, height)
- phase_ana = an.phase_mag_disc(dim, res, phi, center, radius, height, b_0)
+ phase_ana = an.phase_mag_disc(dim, a, phi, center, radius, height, b_0)
elif geometry == 'sphere':
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)
+ phase_ana = an.phase_mag_sphere(dim, a, phi, center, radius, b_0)
# Project the magnetization data:
- mag_data = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
+ mag_data = MagData(a, 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)
+ phase_map_ana = PhaseMap(a, 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
+ phase_map_fft = PhaseMap(a, pm.phase_mag_fourier(a, projection, padding, b_0))
+ 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
+ phase_map_slab = PhaseMap(a, pm.phase_mag_real(a, projection, b_0, 'slab'))
+ 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
+ phase_map_disc = PhaseMap(a, pm.phase_mag_real(a, projection, b_0, 'disc'))
+ 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 Conv):', time.time() - start_time
+ kernel = Kernel(np.shape(projection[0]), a, b_0, 'slab')
+ phase_map_slab_conv = PhaseMap(a, pm.phase_mag(a, projection, b_0, kernel))
+ 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
+ phase_map_disc_conv = PhaseMap(a, pm.phase_mag(a, projection, b_0))
+ print 'Time for real space approach (Disc Conv): ', 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
+ kernel = Kernel(np.shape(projection[0]), a, b_0, 'slab')
+ 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_disc_fast = PhaseMap(res, phase)
- print 'Time for real space approach (Disc Fast):', time.time() - start_time
+ phase_map_slab_fast = PhaseMap(a, pm.phase_mag(a, projection, b_0, kernel))
+ print 'Time for real space approach (Slab, no kernel):', time.time() - start_time
+ start_time = time.time()
+ kernel = Kernel(np.shape(projection[0]), a, b_0)
+ print 'Time for kernel (Disc): ', time.time() - start_time
+ start_time = time.time()
+ phase_map_disc_fast = PhaseMap(a, pm.phase_mag(a, projection, b_0, kernel))
+ print 'Time for real space approach (Disc, no kernel):', 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')
@@ -101,11 +98,11 @@ def compare_methods():
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)
- phase_map_diff_disc = PhaseMap(res, phase_map_ana.phase-phase_map_disc.phase)
- phase_map_diff_slab_conv = PhaseMap(res, phase_map_ana.phase-phase_map_slab_conv.phase)
- phase_map_diff_disc_conv = PhaseMap(res, phase_map_ana.phase-phase_map_disc_conv.phase)
+ phase_map_diff_fft = PhaseMap(a, phase_map_ana.phase-phase_map_fft.phase)
+ phase_map_diff_slab = PhaseMap(a, phase_map_ana.phase-phase_map_slab.phase)
+ phase_map_diff_disc = PhaseMap(a, phase_map_ana.phase-phase_map_disc.phase)
+ phase_map_diff_slab_conv = PhaseMap(a, phase_map_ana.phase-phase_map_slab_conv.phase)
+ phase_map_diff_disc_conv = PhaseMap(a, phase_map_ana.phase-phase_map_disc_conv.phase)
RMS_fft = np.std(phase_map_diff_fft.phase)
RMS_slab = np.std(phase_map_diff_slab.phase)
RMS_disc = np.std(phase_map_diff_disc.phase)
diff --git a/scripts/compare methods/compare_pixel_fields.py b/scripts/compare methods/compare_pixel_fields.py
index a3f680504c0d7765b74d50f4ac3db6298b48dc26..b24e54a2187016ded9fe506059ba93d2fcf19fdf 100644
--- a/scripts/compare methods/compare_pixel_fields.py
+++ b/scripts/compare methods/compare_pixel_fields.py
@@ -1,6 +1,6 @@
#! python
# -*- coding: utf-8 -*-
-"""Compare the phase map of one pixel for different real space approaches."""
+"""Compare the phase map of one pixel for different approaches."""
import pdb
@@ -31,7 +31,7 @@ def compare_pixel_fields():
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
- res = 1.0 # in nm
+ a = 1.0 # in nm
phi = 0 # in rad
dim = (1, 64, 64)
pixel = (0, int(dim[1]/2), int(dim[2]/2))
@@ -40,46 +40,46 @@ def compare_pixel_fields():
def get_fourier_kernel():
PHI_0 = 2067.83 # magnetic flux in T*nm²
b_0 = 1
- coeff = - 1j * b_0 * res**2 / (2*PHI_0)
- nyq = 1 / res # nyquist frequency
+ coeff = - 1j * b_0 * a**2 / (2*PHI_0)
+ nyq = 1 / a # nyquist frequency
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) #* (8*(res/2)**3)*np.sinc(res/2*f_uu)*np.sinc(res/2*f_vv)
+ phase_fft = coeff * f_vv / (f_uu**2 + f_vv**2 + 1e-30) #* (8*(a/2)**3)*np.sinc(a/2*f_uu)*np.sinc(a/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))
return phase_fft_kernel
# Create magnetic data, project it, get the phase map and display the holography image:
- mag_data = MagData(res, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, pixel), phi))
+ mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, pixel), phi))
mag_data.save_to_llg(directory + '/mag_dist_single_pixel.txt')
projection = pj.simple_axis_projection(mag_data)
# Kernel of a disc in real space:
- phase_map_disc = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc'), 'mrad')
+ phase_map_disc = PhaseMap(a, pm.phase_mag_real(a, projection, 'disc'), 'mrad')
# phase_map_disc.display('Phase of one Pixel (Disc)', limit=limit)
# Kernel of a slab in real space:
- phase_map_slab = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'), 'mrad')
+ phase_map_slab = PhaseMap(a, pm.phase_mag_real(a, projection, 'slab'), 'mrad')
# phase_map_slab.display('Phase of one Pixel (Slab)', limit=limit)
# Kernel of the Fourier method:
- phase_map_fft = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding=0), 'mrad')
+ phase_map_fft = PhaseMap(a, pm.phase_mag_fourier(a, 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 = PhaseMap(a, 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_disc.phase - phase_map_fft.phase, 'mrad')
+ phase_map_diff = PhaseMap(a, phase_map_disc.phase - phase_map_fft.phase, 'mrad')
# phase_map_diff.display('Phase difference of one Pixel (Disc - Fourier)')
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 = PhaseMap(a, 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 = PhaseMap(a, phase_inv_disc)
# phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, Power)')
import matplotlib.pyplot as plt
@@ -87,7 +87,7 @@ def compare_pixel_fields():
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1)
- x = np.linspace(-dim[1]/res/2, dim[1]/res/2-1, dim[1])
+ x = np.linspace(-dim[1]/a/2, dim[1]/a/2-1, 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')
@@ -110,11 +110,11 @@ def compare_pixel_fields():
phase = np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0)
- nyq = 0.5 / res
+ nyq = 0.5 / a
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)
+ kern = 4*(a/2)**2*np.sinc(f_uu*a/2)*np.sinc(f_vv*a/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])
diff --git a/scripts/create distributions/create_alternating_filaments.py b/scripts/create distributions/create_alternating_filaments.py
index 3c4073067ae282186293f2adee2837c452a1ccbb..08e77fe9c647dd08ec6320ac98ad09bb2b3a318f 100644
--- a/scripts/create distributions/create_alternating_filaments.py
+++ b/scripts/create distributions/create_alternating_filaments.py
@@ -28,11 +28,11 @@ def create_alternating_filaments():
# Input parameters:
filename = directory + '/mag_dist_alt_filaments.txt'
dim = (1, 21, 21) # in px (z, y, x)
- res = 10.0 # in nm
+ a = 10.0 # in nm
phi = pi/2
spacing = 5
# Create empty MagData object:
- mag_data = MagData(res)
+ mag_data = MagData(a)
count = int((dim[1]-1) / spacing) + 1
for i in range(count):
pos = i * spacing
diff --git a/scripts/create distributions/create_core_shell_disc.py b/scripts/create distributions/create_core_shell_disc.py
index 41d1231ad5378f1b51339c77abb3037e732b574c..a8720895377448ec086316acf70b43cbc93fad5e 100644
--- a/scripts/create distributions/create_core_shell_disc.py
+++ b/scripts/create distributions/create_core_shell_disc.py
@@ -1,3 +1,4 @@
+#! python
# -*- coding: utf-8 -*-
"""Create a core-shell disc."""
@@ -32,7 +33,7 @@ def create_core_shell_disc():
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_core_shell_disc.txt'
- res = 1.0 # in nm
+ a = 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)
@@ -43,7 +44,7 @@ def create_core_shell_disc():
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 = MagData(a, 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')
@@ -51,8 +52,8 @@ def create_core_shell_disc():
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'))
+ phase_map_z = PhaseMap(a, pm.phase_mag(a, projection_z))
+ phase_map_x = PhaseMap(a, pm.phase_mag(a, projection_x))
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)')
diff --git a/scripts/create distributions/create_core_shell_sphere.py b/scripts/create distributions/create_core_shell_sphere.py
index dd6c0aaba8fc8d3f805d846a26ad18bb38db7dd2..e35aef9c6256a45442ae62deb318c19a90ea6a94 100644
--- a/scripts/create distributions/create_core_shell_sphere.py
+++ b/scripts/create distributions/create_core_shell_sphere.py
@@ -1,3 +1,4 @@
+#! python
# -*- coding: utf-8 -*-
"""Create a core-shell disc."""
@@ -32,7 +33,7 @@ def create_core_shell_sphere():
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_core_shell_sphere.txt'
- res = 1.0 # in nm
+ a = 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)
@@ -43,7 +44,7 @@ def create_core_shell_sphere():
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 = MagData(a, 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')
@@ -51,8 +52,8 @@ def create_core_shell_sphere():
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'))
+ phase_map_z = PhaseMap(a, pm.phase_mag(a, projection_z))
+ phase_map_x = PhaseMap(a, pm.phase_mag(a, projection_x))
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)')
diff --git a/scripts/create distributions/create_logo.py b/scripts/create distributions/create_logo.py
index f800275d4aa7ce88d1f837e51a07a958e53c3270..fe6cde5f940dadc1f24e5ed954c4a6684a23ec64 100644
--- a/scripts/create distributions/create_logo.py
+++ b/scripts/create distributions/create_logo.py
@@ -27,7 +27,7 @@ def create_logo():
'''
# Input parameters:
- res = 10.0 # in nm
+ a = 10.0 # in nm
phi = -pi/2 # in rad
density = 10
dim = (1, 128, 128)
@@ -41,10 +41,10 @@ def create_logo():
right = np.fliplr(left)
mag_shape[0, ...] = np.logical_and(np.logical_and(left, right), bottom)
# Create magnetic data, project it, get the phase map and display the holography image:
- mag_data = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
+ mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
mag_data.quiver_plot()
projection = pj.simple_axis_projection(mag_data)
- phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
+ phase_map = PhaseMap(a, pm.phase_mag(a, projection))
hi.display(hi.holo_image(phase_map, density), 'PYRAMID - LOGO', interpolation='bilinear')
diff --git a/scripts/create distributions/create_multiple_samples.py b/scripts/create distributions/create_multiple_samples.py
index 9590e784a2a1892e9ea9a024b08b146c73cf99a5..1236fb1efc8a9e87b962ae2605b2295bdb8ddd00 100644
--- a/scripts/create distributions/create_multiple_samples.py
+++ b/scripts/create distributions/create_multiple_samples.py
@@ -1,3 +1,4 @@
+#! python
# -*- coding: utf-8 -*-
"""Create random magnetic distributions."""
@@ -30,7 +31,7 @@ def create_multiple_samples():
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_multiple_samples.txt'
- res = 10.0 # nm
+ a = 10.0 # nm
dim = (64, 128, 128)
# Slab:
center = (32, 32, 32) # in px (z, y, x), index starts with 0!
@@ -46,7 +47,7 @@ def create_multiple_samples():
radius = 24 # in px
mag_shape_sphere = mc.Shapes.sphere(dim, center, radius)
# Create empty MagData object and add magnetized objects:
- mag_data = MagData(res)
+ mag_data = MagData(a)
mag_data.add_magnitude(mc.create_mag_dist_homog(mag_shape_slab, pi/4))
mag_data.add_magnitude(mc.create_mag_dist_homog(mag_shape_disc, pi/2))
mag_data.add_magnitude(mc.create_mag_dist_homog(mag_shape_sphere, pi))
@@ -54,7 +55,7 @@ def create_multiple_samples():
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
projection = pj.simple_axis_projection(mag_data)
- phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
+ phase_map = PhaseMap(a, pm.phase_mag(a, projection))
phase_map.display()
hi.display(hi.holo_image(phase_map, 0.5))
diff --git a/scripts/create distributions/create_random_pixels.py b/scripts/create distributions/create_random_pixels.py
index e6eea7ad67dda36ed75cc9b2b9e8fbe103db9ba5..570c427dd493acd4a2ca4261c94e77ce49b61bdc 100644
--- a/scripts/create distributions/create_random_pixels.py
+++ b/scripts/create distributions/create_random_pixels.py
@@ -1,3 +1,4 @@
+#! python
# -*- coding: utf-8 -*-
"""Create random magnetic distributions."""
@@ -7,7 +8,6 @@ import traceback
import sys
import os
-import numpy as np
from numpy import pi
import pyramid.magcreator as mc
@@ -34,21 +34,21 @@ def create_random_pixels():
# Input parameters:
count = 10
dim = (1, 128, 128)
- res = 10 # in nm
+ a = 10 # in nm
rnd.seed(12)
# Create empty MagData object and add pixels:
- mag_data = MagData(res)
+ mag_data = MagData(a)
for i in range(count):
pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
mag_shape = mc.Shapes.pixel(dim, pixel)
phi = 2 * pi * rnd.random()
- magnitude = 1 # TODO: rnd.random()
+ magnitude = 1
mag_data.add_magnitude(mc.create_mag_dist_homog(mag_shape, phi, magnitude=magnitude))
# Plot magnetic distribution, phase map and holographic contour map:
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
projection = pj.simple_axis_projection(mag_data)
- phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
+ phase_map = PhaseMap(a, pm.phase_mag(a, projection))
hi.display(hi.holo_image(phase_map, 10))
diff --git a/scripts/create distributions/create_random_slabs.py b/scripts/create distributions/create_random_slabs.py
index c54b27c9635c8de05d00d8f5fe3bf0ba6a0b1520..513fcb97765336d01d67c9c97838e7fa3b87caac 100644
--- a/scripts/create distributions/create_random_slabs.py
+++ b/scripts/create distributions/create_random_slabs.py
@@ -1,3 +1,4 @@
+#! python
# -*- coding: utf-8 -*-
"""Create random magnetic distributions."""
@@ -7,7 +8,6 @@ import traceback
import sys
import os
-import numpy as np
from numpy import pi
import pyramid.magcreator as mc
@@ -33,11 +33,11 @@ def create_random_slabs():
filename = directory + '/mag_dist_random_pixels.txt'
count = 10
dim = (1, 128, 128)
- res = 10 # in nm
+ a = 10 # in nm
rnd.seed(42)
w_max = 10
# Create empty MagData object and add slabs:
- mag_data = MagData(res)
+ mag_data = MagData(a)
for i in range(count):
width = (1, rnd.randint(1, w_max), rnd.randint(1, w_max))
center = (rnd.randrange(int(width[0]/2), dim[0]-int(width[0]/2)),
@@ -51,7 +51,7 @@ def create_random_slabs():
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
projection = pj.simple_axis_projection(mag_data)
- phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
+ phase_map = PhaseMap(a, pm.phase_mag(a, projection))
hi.display(hi.holo_image(phase_map, 10))
diff --git a/scripts/create distributions/create_sample.py b/scripts/create distributions/create_sample.py
index cfc625a49e2ea419a1a13f1b4d5a222bba81f80f..1478606ee154da4541467c9be332ff80a1e43272 100644
--- a/scripts/create distributions/create_sample.py
+++ b/scripts/create distributions/create_sample.py
@@ -1,3 +1,4 @@
+#! python
# -*- coding: utf-8 -*-
"""Create magnetic distributions with simple geometries."""
@@ -26,10 +27,10 @@ def create_sample():
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
- key = 'pixel'
+ key = 'sphere'
filename = directory + '/mag_dist_' + key + '.txt'
- dim = (1, 1, 5) # in px (z, y, x)
- res = 1.0 # in nm
+ dim = (64, 64, 64) # in px (z, y, x)
+ a = 1.0 # in nm
phi = pi/2
# Geometry parameters:
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!
@@ -51,7 +52,7 @@ def create_sample():
mag_shape = mc.Shapes.pixel(dim, pixel)
# Create magnetic distribution
magnitude = mc.create_mag_dist_homog(mag_shape, phi)
- mag_data = MagData(res, magnitude)
+ mag_data = MagData(a, magnitude)
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
diff --git a/scripts/create distributions/create_vortex.py b/scripts/create distributions/create_vortex.py
index 794a02b920c612c3f34e8dc7d745883880c752ce..9b6879aadd1f1f208006470d10f31aa65627bcf4 100644
--- a/scripts/create distributions/create_vortex.py
+++ b/scripts/create distributions/create_vortex.py
@@ -1,3 +1,4 @@
+#! python
# -*- coding: utf-8 -*-
"""Create the Pyramid-Logo."""
@@ -30,7 +31,7 @@ def create_vortex():
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_vortex.txt'
- res = 10.0 # in nm
+ a = 10.0 # in nm
density = 1
dim = (64, 64, 64)
center = (int(dim[0]/2)-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
@@ -38,12 +39,11 @@ def create_vortex():
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))
+ mag_data = MagData(a, mc.create_mag_dist_vortex(mag_shape))
mag_data.quiver_plot()
- mag_data.quiver_plot3d()
mag_data.save_to_llg(filename)
projection = pj.simple_axis_projection(mag_data)
- phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
+ phase_map = PhaseMap(a, pm.phase_mag(a, projection))
hi.display_combined(phase_map, density, 'Vortex State')
phase_slice = phase_map.phase[dim[1]/2, :]
fig = plt.figure()
diff --git a/scripts/get_jacobi_projection.py b/scripts/get_jacobi_projection.py
new file mode 100644
index 0000000000000000000000000000000000000000..079c39ab40854ab769c18e05fda7dcc4ea0259a4
--- /dev/null
+++ b/scripts/get_jacobi_projection.py
@@ -0,0 +1,20 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Nov 29 16:45:13 2013
+
+@author: Jan
+"""
+
+import pyramid.magcreator as mc
+from pyramid.magdata import MagData
+import pyramid.projector as pj
+from pyramid.projector import Projection
+
+
+a = 1.0
+dim = (2, 2, 2)
+px = (0, 0, 0)
+mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, px), 0))
+
+proj = pj.simple_axis_projection(mag_data)
+
diff --git a/scripts/gui/gui_create_logo.py b/scripts/gui/gui_create_logo.py
index c7b7e52a843d974a4adff3165f93946c3aa236a4..071d866f1f22fae05630fe88847f4dff04a87afc 100644
--- a/scripts/gui/gui_create_logo.py
+++ b/scripts/gui/gui_create_logo.py
@@ -107,7 +107,7 @@ def create_logo(dim, axis):
'''
# Input parameters:
- res = 10.0 # in nm
+ a = 10.0 # in nm
phi = -pi/2 # in rad
density = 10
# Create magnetic shape:
@@ -120,9 +120,9 @@ def create_logo(dim, axis):
right = np.fliplr(left)
mag_shape[0, ...] = np.logical_and(np.logical_and(left, right), bottom)
# Create magnetic data, project it, get the phase map and display the holography image:
- mag_data = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
+ mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
projection = pj.simple_axis_projection(mag_data)
- phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
+ phase_map = PhaseMap(a, pm.phase_mag_real(a, projection, 'slab'))
hi.display(hi.holo_image(phase_map, density), 'PYRAMID - LOGO', axis)
diff --git a/scripts/kernel.py b/scripts/kernel.py
index 6fc4bd05ed4b3a8246bfb363a88e878778b673e4..4e0291bfa886d364f342df8b4c2f48cb91e5bd12 100644
--- a/scripts/kernel.py
+++ b/scripts/kernel.py
@@ -5,12 +5,14 @@ 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
@@ -28,7 +30,6 @@ 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')
diff --git a/scripts/paper 1/ch5-0-evaluation_and_comparison.py b/scripts/paper 1/ch5-0-evaluation_and_comparison.py
index 1f53da28d2bb45814db3c1444469530f93ae05e2..bdae2361b751c18660361391eae14d978f1a4d0d 100644
--- a/scripts/paper 1/ch5-0-evaluation_and_comparison.py
+++ b/scripts/paper 1/ch5-0-evaluation_and_comparison.py
@@ -21,7 +21,7 @@ from pyramid.magdata import MagData
import matplotlib.pyplot as plt
-from matplotlib.ticker import FixedFormatter, IndexLocator, LinearLocator
+from matplotlib.ticker import FixedFormatter, IndexLocator
def run():
@@ -60,7 +60,7 @@ def run():
else:
print '--CREATE MAGNETIC DISTRIBUTIONS'
# Input parameters:
- res = 1.0 # in nm
+ a = 1.0 # in nm
phi = pi/2
dim = (16, 128, 128) # in px (z, y, x)
# Create magnetic shape:
@@ -69,10 +69,10 @@ def run():
height = dim[0]/2 # in px
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 = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
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 = MagData(a, mc.create_mag_dist_vortex(mag_shape, center))
mag_data_vort.scale_down(2)
# Mayavi-Plots:
mag_data_disc.quiver_plot3d()
diff --git a/scripts/paper 1/ch5-1-evaluation_real_space.py b/scripts/paper 1/ch5-1-evaluation_real_space.py
index 2bf7b23d18e2a3ee5e2f63f6a27d910e62e03795..babbdcbac83173d0d34f09b82c35f37221461444 100644
--- a/scripts/paper 1/ch5-1-evaluation_real_space.py
+++ b/scripts/paper 1/ch5-1-evaluation_real_space.py
@@ -47,7 +47,7 @@ def run():
print 'CH5-1 ANALYTIC SOLUTIONS'
# Input parameters:
- res = 0.125 # in nm
+ a = 0.125 # in nm
phi = pi/2
dim = (128, 1024, 1024) # in px (z, y, x)
density = 100
@@ -57,10 +57,10 @@ def run():
height = dim[0]/2 # in px
print '--CALCULATE ANALYTIC SOLUTIONS'
# Get analytic solution:
- phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
- phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
- phase_map_ana_disc = PhaseMap(res, phase_ana_disc)
- phase_map_ana_vort = PhaseMap(res, phase_ana_vort)
+ phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
+ phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
+ phase_map_ana_disc = PhaseMap(a, phase_ana_disc)
+ phase_map_ana_vort = PhaseMap(a, phase_ana_vort)
print '--PLOT/SAVE ANALYTIC SOLUTIONS'
hi.display_combined(phase_map_ana_disc, density,
'Analytic solution: hom. magn. disc', 'bilinear')
@@ -85,7 +85,7 @@ def run():
print 'CH5-1 PHASE SLICES REAL SPACE'
# Input parameters:
- res = 0.25 # in nm
+ a = 0.25 # in nm
phi = pi/2
density = 100
dim = (64, 512, 512) # in px (z, y, x)
@@ -108,8 +108,8 @@ def run():
y_d = []
dy_d = []
# Analytic solution:
- L = dim[1] * res # in px/nm
- Lz = 0.5 * dim[0] * res # in px/nm
+ L = dim[1] * a # in px/nm
+ Lz = 0.5 * dim[0] * a # in px/nm
R = 0.25 * L # in px/nm
x0 = L / 2 # in px/nm
@@ -123,16 +123,15 @@ def run():
y_d.append(F_disc(x_d[0]))
dy_d.append(np.zeros_like(x_d[0]))
# Create MagData (Disc):
- mag_data_disc = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
+ mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
for i in range(5):
mag_data_disc.scale_down()
- print '----res =', mag_data_disc.res, 'nm', 'dim =', mag_data_disc.dim
+ print '----a =', mag_data_disc.a, 'nm', 'dim =', mag_data_disc.dim
projection = pj.simple_axis_projection(mag_data_disc)
- phase_map = PhaseMap(mag_data_disc.res,
- pm.phase_mag_real(mag_data_disc.res, projection, 'slab'))
- hi.display_combined(phase_map, density, 'Disc, res = {} nm'.format(mag_data_disc.res))
- x_d.append(np.linspace(mag_data_disc.res * 0.5,
- mag_data_disc.res * (mag_data_disc.dim[1]-0.5),
+ phase_map = PhaseMap(mag_data_disc.a, pm.phase_mag(mag_data_disc.a, projection))
+ hi.display_combined(phase_map, density, 'Disc, a = {} nm'.format(mag_data_disc.a))
+ x_d.append(np.linspace(mag_data_disc.a * 0.5,
+ mag_data_disc.a * (mag_data_disc.dim[1]-0.5),
mag_data_disc.dim[1]))
slice_pos = int(mag_data_disc.dim[1]/2)
y_d.append(phase_map.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
@@ -143,8 +142,8 @@ def run():
y_v = []
dy_v = []
# Analytic solution:
- L = dim[1] * res # in px/nm
- Lz = 0.5 * dim[0] * res # in px/nm
+ L = dim[1] * a # in px/nm
+ Lz = 0.5 * dim[0] * a # in px/nm
R = 0.25 * L # in px/nm
x0 = L / 2 # in px/nm
@@ -157,16 +156,15 @@ def run():
y_v.append(F_vort(x_v[0]))
dy_v.append(np.zeros_like(x_v[0]))
# Create MagData (Vortex):
- mag_data_vort = MagData(res, mc.create_mag_dist_vortex(mag_shape))
+ mag_data_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape))
for i in range(5):
mag_data_vort.scale_down()
- print '----res =', mag_data_vort.res, 'nm', 'dim =', mag_data_vort.dim
+ print '----a =', mag_data_vort.a, 'nm', 'dim =', mag_data_vort.dim
projection = pj.simple_axis_projection(mag_data_vort)
- phase_map = PhaseMap(mag_data_vort.res,
- pm.phase_mag_real(mag_data_vort.res, projection, 'slab'))
- hi.display_combined(phase_map, density, 'Disc, res = {} nm'.format(mag_data_vort.res))
- x_v.append(np.linspace(mag_data_vort.res * 0.5,
- mag_data_vort.res * (mag_data_vort.dim[1]-0.5),
+ phase_map = PhaseMap(mag_data_vort.a, pm.phase_mag(mag_data_vort.a, projection))
+ hi.display_combined(phase_map, density, 'Disc, a = {} nm'.format(mag_data_vort.a))
+ x_v.append(np.linspace(mag_data_vort.a * 0.5,
+ mag_data_vort.a * (mag_data_vort.dim[1]-0.5),
mag_data_vort.dim[1]))
slice_pos = int(mag_data_vort.dim[1]/2)
y_v.append(phase_map.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
@@ -296,7 +294,7 @@ def run():
print 'CH5-1 PHASE DIFFERENCES REAL SPACE'
# Input parameters:
- res = 1.0 # in nm
+ a = 1.0 # in nm
phi = pi/2
dim = (16, 128, 128) # 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 with 0!
@@ -310,16 +308,16 @@ def run():
else:
print '--CREATE MAGNETIC DISTRIBUTIONS'
# Create magnetic shape (4 times the size):
- res_big = res / 2
+ a_big = a / 2
dim_big = (dim[0]*2, dim[1]*2, dim[2]*2)
center_big = (dim_big[0]/2-0.5, dim_big[1]/2.-0.5, dim_big[2]/2.-0.5)
radius_big = dim_big[1]/4 # in px
height_big = dim_big[0]/2 # in px
mag_shape = mc.Shapes.disc(dim_big, center_big, radius_big, height_big)
# Create MagData (4 times the size):
- mag_data_disc = MagData(res_big, mc.create_mag_dist_homog(mag_shape, phi))
- mag_data_vort = MagData(res_big, mc.create_mag_dist_vortex(mag_shape, center_big))
- # Scale mag_data, resolution and dimensions:
+ mag_data_disc = MagData(a_big, mc.create_mag_dist_homog(mag_shape, phi))
+ mag_data_vort = MagData(a_big, mc.create_mag_dist_vortex(mag_shape, center_big))
+ # Scale mag_data, grid spacing and dimensions:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
print '--SAVE MAGNETIC DISTRIBUTIONS'
@@ -331,18 +329,18 @@ def run():
projection_disc = pj.simple_axis_projection(mag_data_disc)
projection_vort = pj.simple_axis_projection(mag_data_vort)
# Get analytic solutions:
- phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
- phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
+ phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
+ phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# Create norm for the plots:
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Calculations (Disc):
- phase_num_disc = pm.phase_mag_real(res, projection_disc, 'slab')
- phase_diff_disc = PhaseMap(res, (phase_num_disc-phase_ana_disc), 'mrad')
+ phase_num_disc = pm.phase_mag(a, projection_disc)
+ phase_diff_disc = PhaseMap(a, (phase_num_disc-phase_ana_disc), 'mrad')
RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Calculations (Vortex):
- phase_num_vort = pm.phase_mag_real(res, projection_vort, 'slab')
- phase_diff_vort = PhaseMap(res, (phase_num_vort-phase_ana_vort), 'mrad')
+ phase_num_vort = pm.phase_mag(a, projection_vort)
+ phase_diff_vort = PhaseMap(a, (phase_num_vort-phase_ana_vort), 'mrad')
RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
print '--PLOT/SAVE PHASE DIFFERENCES'
diff --git a/scripts/paper 1/ch5-2-evaluation_fourier_space.py b/scripts/paper 1/ch5-2-evaluation_fourier_space.py
index 13f4e2f03146556e545ec8c0231a18ca83845309..4a1ef4596f874aaf1cb39fb8d9be64a597f09651 100644
--- a/scripts/paper 1/ch5-2-evaluation_fourier_space.py
+++ b/scripts/paper 1/ch5-2-evaluation_fourier_space.py
@@ -47,7 +47,7 @@ def run():
print 'CH5-1 PHASE SLICES FOURIER SPACE'
# Input parameters:
- res = 0.25 # in nm
+ a = 0.25 # in nm
phi = pi/2
density = 100
dim = (64, 512, 512) # in px (z, y, x)
@@ -72,8 +72,8 @@ def run():
dy_d0 = []
dy_d10 = []
# Analytic solution:
- L = dim[1] * res # in px/nm
- Lz = 0.5 * dim[0] * res # in px/nm
+ L = dim[1] * a # in px/nm
+ Lz = 0.5 * dim[0] * a # in px/nm
R = 0.25 * L # in px/nm
x0 = L / 2 # in px/nm
@@ -89,23 +89,23 @@ def run():
dy_d0.append(np.zeros_like(x_d[0]))
dy_d10.append(np.zeros_like(x_d[0]))
# Create MagData (Disc):
- mag_data_disc = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
+ mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
for i in range(5):
mag_data_disc.scale_down()
- print '----res =', mag_data_disc.res, 'nm', 'dim =', mag_data_disc.dim
+ print '----a =', mag_data_disc.a, 'nm', 'dim =', mag_data_disc.dim
projection = pj.simple_axis_projection(mag_data_disc)
- phase_map0 = PhaseMap(mag_data_disc.res,
- pm.phase_mag_fourier(mag_data_disc.res, projection,
+ phase_map0 = PhaseMap(mag_data_disc.a,
+ pm.phase_mag_fourier(mag_data_disc.a, projection,
padding=0), 'mrad')
- phase_map10 = PhaseMap(mag_data_disc.res,
- pm.phase_mag_fourier(mag_data_disc.res, projection,
+ phase_map10 = PhaseMap(mag_data_disc.a,
+ pm.phase_mag_fourier(mag_data_disc.a, projection,
padding=10), 'mrad')
hi.display_combined(phase_map0, density,
- 'Disc, res = {} nm'.format(mag_data_disc.res))
+ 'Disc, a = {} nm'.format(mag_data_disc.a))
hi.display_combined(phase_map10, density,
- 'Disc, res = {} nm'.format(mag_data_disc.res))
- x_d.append(np.linspace(mag_data_disc.res * 0.5,
- mag_data_disc.res * (mag_data_disc.dim[1]-0.5),
+ 'Disc, a = {} nm'.format(mag_data_disc.a))
+ x_d.append(np.linspace(mag_data_disc.a * 0.5,
+ mag_data_disc.a * (mag_data_disc.dim[1]-0.5),
mag_data_disc.dim[1]))
slice_pos = int(mag_data_disc.dim[1]/2)
y_d0.append(phase_map0.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
@@ -120,8 +120,8 @@ def run():
dy_v0 = []
dy_v10 = []
# Analytic solution:
- L = dim[1] * res # in px/nm
- Lz = 0.5 * dim[0] * res # in px/nm
+ L = dim[1] * a # in px/nm
+ Lz = 0.5 * dim[0] * a # in px/nm
R = 0.25 * L # in px/nm
x0 = L / 2 # in px/nm
@@ -136,23 +136,23 @@ def run():
dy_v0.append(np.zeros_like(x_v[0]))
dy_v10.append(np.zeros_like(x_v[0]))
# Create MagData (Vortex):
- mag_data_vort = MagData(res, mc.create_mag_dist_vortex(mag_shape))
+ mag_data_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape))
for i in range(5):
mag_data_vort.scale_down()
- print '----res =', mag_data_vort.res, 'nm', 'dim =', mag_data_vort.dim
+ print '----a =', mag_data_vort.a, 'nm', 'dim =', mag_data_vort.dim
projection = pj.simple_axis_projection(mag_data_vort)
- phase_map0 = PhaseMap(mag_data_vort.res,
- pm.phase_mag_fourier(mag_data_vort.res, projection,
+ phase_map0 = PhaseMap(mag_data_vort.a,
+ pm.phase_mag_fourier(mag_data_vort.a, projection,
padding=0), 'mrad')
- phase_map10 = PhaseMap(mag_data_vort.res,
- pm.phase_mag_fourier(mag_data_vort.res, projection,
+ phase_map10 = PhaseMap(mag_data_vort.a,
+ pm.phase_mag_fourier(mag_data_vort.a, projection,
padding=10), 'mrad')
hi.display_combined(phase_map0, density,
- 'Disc, res = {} nm'.format(mag_data_vort.res))
+ 'Disc, a = {} nm'.format(mag_data_vort.a))
hi.display_combined(phase_map10, density,
- 'Disc, res = {} nm'.format(mag_data_vort.res))
- x_v.append(np.linspace(mag_data_vort.res * 0.5,
- mag_data_vort.res * (mag_data_vort.dim[1]-0.5),
+ 'Disc, a = {} nm'.format(mag_data_vort.a))
+ x_v.append(np.linspace(mag_data_vort.a * 0.5,
+ mag_data_vort.a * (mag_data_vort.dim[1]-0.5),
mag_data_vort.dim[1]))
slice_pos = int(mag_data_vort.dim[1]/2)
y_v0.append(phase_map0.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
@@ -400,7 +400,7 @@ def run():
print 'CH5-2 PHASE DIFFERENCES FOURIER SPACE'
# Input parameters:
- res = 1.0 # in nm
+ a = 1.0 # in nm
phi = pi/2
dim = (16, 128, 128) # 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 with 0!
@@ -414,16 +414,16 @@ def run():
else:
print '--CREATE MAGNETIC DISTRIBUTIONS'
# Create magnetic shape (4 times the size):
- res_big = res / 2
+ a_big = a / 2
dim_big = (dim[0]*2, dim[1]*2, dim[2]*2)
center_big = (dim_big[0]/2-0.5, dim_big[1]/2.-0.5, dim_big[2]/2.-0.5)
radius_big = dim_big[1]/4 # in px
height_big = dim_big[0]/2 # in px
mag_shape = mc.Shapes.disc(dim_big, center_big, radius_big, height_big)
# Create MagData (4 times the size):
- mag_data_disc = MagData(res_big, mc.create_mag_dist_homog(mag_shape, phi))
- mag_data_vort = MagData(res_big, mc.create_mag_dist_vortex(mag_shape, center_big))
- # Scale mag_data, resolution and dimensions:
+ mag_data_disc = MagData(a_big, mc.create_mag_dist_homog(mag_shape, phi))
+ mag_data_vort = MagData(a_big, mc.create_mag_dist_vortex(mag_shape, center_big))
+ # Scale mag_data, grid spacing and dimensions:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
print '--SAVE MAGNETIC DISTRIBUTIONS'
@@ -435,8 +435,8 @@ def run():
projection_disc = pj.simple_axis_projection(mag_data_disc)
projection_vort = pj.simple_axis_projection(mag_data_vort)
# Get analytic solutions:
- phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
- phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
+ phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
+ phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
@@ -445,15 +445,15 @@ def run():
bounds = np.array([-100, -50, -25, -5, 0, 5, 25, 50, 100])
norm = BoundaryNorm(bounds, RdBu.N)
# Disc:
- phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=0)
- phase_diff_disc = PhaseMap(res, (phase_num_disc-phase_ana_disc), 'mrad')
+ phase_num_disc = pm.phase_mag_fourier(a, projection_disc, padding=0)
+ phase_diff_disc = PhaseMap(a, (phase_num_disc-phase_ana_disc), 'mrad')
RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
phase_diff_disc.display('Homog. magn. disc, RMS = {:3.2f} mrad'.format(RMS_disc),
axis=axes[0], limit=np.max(bounds), norm=norm)
axes[0].set_aspect('equal')
# Vortex:
- phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=0)
- phase_diff_vort = PhaseMap(res, (phase_num_vort-phase_ana_vort), 'mrad')
+ phase_num_vort = pm.phase_mag_fourier(a, projection_vort, padding=0)
+ phase_diff_vort = PhaseMap(a, (phase_num_vort-phase_ana_vort), 'mrad')
RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
phase_diff_vort.display('Vortex state disc, RMS = {:3.2f} mrad'.format(RMS_vort),
axis=axes[1], limit=np.max(bounds), norm=norm)
@@ -471,15 +471,15 @@ def run():
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Disc:
- phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=10)
- phase_diff_disc = PhaseMap(res, (phase_num_disc-phase_ana_disc), 'mrad')
+ phase_num_disc = pm.phase_mag_fourier(a, projection_disc, padding=10)
+ phase_diff_disc = PhaseMap(a, (phase_num_disc-phase_ana_disc), 'mrad')
RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
phase_diff_disc.display('Homog. magn. disc, RMS = {:3.2f} mrad'.format(RMS_disc),
axis=axes[0], limit=np.max(bounds), norm=norm)
axes[0].set_aspect('equal')
# Vortex:
- phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=10)
- phase_diff_vort = PhaseMap(res, (phase_num_vort-phase_ana_vort), 'mrad')
+ phase_num_vort = pm.phase_mag_fourier(a, projection_vort, padding=10)
+ phase_diff_vort = PhaseMap(a, (phase_num_vort-phase_ana_vort), 'mrad')
RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
phase_diff_vort.display('Vortex state disc, RMS = {:3.2f} mrad'.format(RMS_vort),
axis=axes[1], limit=np.max(bounds), norm=norm)
@@ -494,7 +494,7 @@ def run():
print 'CH5-2 FOURIER PADDING VARIATION'
# Input parameters:
- res = 1.0 # in nm
+ a = 1.0 # in nm
phi = pi/2
dim = (16, 128, 128) # 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 with 0!
@@ -508,16 +508,16 @@ def run():
else:
print '--CREATE MAGNETIC DISTRIBUTIONS'
# Create magnetic shape (4 times the size):
- res_big = res / 2
+ a_big = a / 2
dim_big = (dim[0]*2, dim[1]*2, dim[2]*2)
center_big = (dim_big[0]/2-0.5, dim_big[1]/2.-0.5, dim_big[2]/2.-0.5)
radius_big = dim_big[1]/4 # in px
height_big = dim_big[0]/2 # in px
mag_shape = mc.Shapes.disc(dim_big, center_big, radius_big, height_big)
# Create MagData (4 times the size):
- mag_data_disc = MagData(res_big, mc.create_mag_dist_homog(mag_shape, phi))
- mag_data_vort = MagData(res_big, mc.create_mag_dist_vortex(mag_shape, center_big))
- # Scale mag_data, resolution and dimensions:
+ mag_data_disc = MagData(a_big, mc.create_mag_dist_homog(mag_shape, phi))
+ mag_data_vort = MagData(a_big, mc.create_mag_dist_vortex(mag_shape, center_big))
+ # Scale mag_data, grid spacing and dimensions:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
print '--SAVE MAGNETIC DISTRIBUTIONS'
@@ -528,8 +528,8 @@ def run():
projection_disc = pj.simple_axis_projection(mag_data_disc)
projection_vort = pj.simple_axis_projection(mag_data_vort)
# Get analytic solutions:
- phase_ana_disc = an.phase_mag_disc(dim, res, phi, center, radius, height)
- phase_ana_vort = an.phase_mag_vortex(dim, res, center, radius, height)
+ phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
+ phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# List of applied paddings:
padding_list = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
@@ -539,16 +539,16 @@ def run():
data_disc[0, :] = padding_list
for (i, padding) in enumerate(padding_list):
key = ', '.join(['Padding->RMS|duration', 'Fourier', 'padding={}'.format(padding_list[i]),
- 'res={}'.format(res), 'dim={}'.format(dim), 'phi={}'.format(phi), 'disc'])
+ 'a={}'.format(a), 'dim={}'.format(dim), 'phi={}'.format(phi), 'disc'])
if key in data_shelve:
data_disc[:, i] = data_shelve[key]
else:
print '----calculate and save padding =', padding_list[i]
start_time = time.time()
- phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding_list[i])
+ phase_num_disc = pm.phase_mag_fourier(a, projection_disc, padding=padding_list[i])
data_disc[2, i] = time.time() - start_time
phase_diff = (phase_num_disc-phase_ana_disc)
- phase_map_diff = PhaseMap(res, phase_diff, 'mrad')
+ phase_map_diff = PhaseMap(a, phase_diff, 'mrad')
phase_map_diff.display()
data_disc[1, i] = np.sqrt(np.mean(phase_map_diff.phase**2))*1E3 # *1E3: rad to mraddd
data_shelve[key] = data_disc[:, i]
@@ -580,7 +580,7 @@ def run():
axes[1].set_xlabel('padding', fontsize=15)
axes[1].set_ylabel('duration [s]', fontsize=15)
axes[1].set_xlim(-0.5, 16.5)
- axes[1].set_ylim(-0.05, 1.5)
+ axes[1].set_ylim(-0.05, 3)
axes[1].xaxis.set_major_locator(MaxNLocator(nbins=10, integer=True))
axes[1].yaxis.set_major_locator(MaxNLocator(nbins=10))
axes[1].tick_params(axis='both', which='major', labelsize=14)
@@ -595,16 +595,16 @@ def run():
data_vort[0, :] = padding_list
for (i, padding) in enumerate(padding_list):
key = ', '.join(['Padding->RMS|duration', 'Fourier', 'padding={}'.format(padding_list[i]),
- 'res={}'.format(res), 'dim={}'.format(dim), 'phi={}'.format(phi), 'vort'])
+ 'a={}'.format(a), 'dim={}'.format(dim), 'phi={}'.format(phi), 'vort'])
if key in data_shelve:
data_vort[:, i] = data_shelve[key]
else:
print '----calculate and save padding =', padding_list[i]
start_time = time.time()
- phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding_list[i])
+ phase_num_vort = pm.phase_mag_fourier(a, projection_vort, padding=padding_list[i])
data_vort[2, i] = time.time() - start_time
phase_diff = (phase_num_vort-phase_ana_vort)
- phase_map_diff = PhaseMap(res, phase_diff, 'mrad')
+ phase_map_diff = PhaseMap(a, phase_diff, 'mrad')
phase_map_diff.display()
data_vort[1, i] = np.sqrt(np.mean(phase_map_diff.phase**2))*1E3 # *1E3: rad to mrad
data_shelve[key] = data_vort[:, i]
@@ -636,7 +636,7 @@ def run():
axes[1].set_xlabel('padding', fontsize=15)
axes[1].set_ylabel('duration [s]', fontsize=15)
axes[1].set_xlim(-0.5, 16.5)
- axes[1].set_ylim(-0.05, 1.5)
+ axes[1].set_ylim(-0.05, 3)
axes[1].xaxis.set_major_locator(MaxNLocator(nbins=10, integer=True))
axes[1].yaxis.set_major_locator(MaxNLocator(nbins=10))
axes[1].tick_params(axis='both', which='major', labelsize=14)
diff --git a/scripts/paper 1/ch5-3-comparison_of_methods.py b/scripts/paper 1/ch5-3-comparison_of_methods.py
index 7ba66098201465433c58a20d551eee4eda1b2a6c..6d230f049699fe388d30d5c117f6512aa225f03e 100644
--- a/scripts/paper 1/ch5-3-comparison_of_methods.py
+++ b/scripts/paper 1/ch5-3-comparison_of_methods.py
@@ -50,7 +50,7 @@ def run():
else:
# Input parameters:
steps = 6
- res = 0.25 # in nm
+ a = 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!
@@ -61,44 +61,44 @@ def run():
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))
+ mag_data_disc = MagData(a, 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))
+ mag_data_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape, center))
# Create Data Arrays:
- res_list = [res*2**i for i in np.linspace(1, steps, steps)]
- data_disc_fourier0 = np.vstack((res_list, np.zeros((2, steps))))
- data_vort_fourier0 = np.vstack((res_list, np.zeros((2, steps))))
- data_disc_fourier1 = np.vstack((res_list, np.zeros((2, steps))))
- data_vort_fourier1 = np.vstack((res_list, np.zeros((2, steps))))
- data_disc_fourier10 = np.vstack((res_list, np.zeros((2, steps))))
- data_vort_fourier10 = np.vstack((res_list, np.zeros((2, steps))))
- data_disc_real_s = np.vstack((res_list, np.zeros((2, steps))))
- data_vort_real_s = np.vstack((res_list, np.zeros((2, steps))))
- data_disc_real_d = np.vstack((res_list, np.zeros((2, steps))))
- data_vort_real_d = np.vstack((res_list, np.zeros((2, steps))))
+ a_list = [a*2**i for i in np.linspace(1, steps, steps)]
+ data_disc_fourier0 = np.vstack((a_list, np.zeros((2, steps))))
+ data_vort_fourier0 = np.vstack((a_list, np.zeros((2, steps))))
+ data_disc_fourier1 = np.vstack((a_list, np.zeros((2, steps))))
+ data_vort_fourier1 = np.vstack((a_list, np.zeros((2, steps))))
+ data_disc_fourier10 = np.vstack((a_list, np.zeros((2, steps))))
+ data_vort_fourier10 = np.vstack((a_list, np.zeros((2, steps))))
+ data_disc_real_s = np.vstack((a_list, np.zeros((2, steps))))
+ data_vort_real_s = np.vstack((a_list, np.zeros((2, steps))))
+ data_disc_real_d = np.vstack((a_list, np.zeros((2, steps))))
+ data_vort_real_d = np.vstack((a_list, np.zeros((2, steps))))
for i in range(steps):
- # Scale mag_data, resolution and dimensions:
+ # Scale mag_data, grid spacing and dimensions:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
dim = mag_data_disc.dim
- res = mag_data_disc.res
+ a = mag_data_disc.a
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 '--a =', a, '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)
+ phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
# Fourier unpadded:
padding = 0
start_time = time.clock()
- phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+ phase_num_disc = pm.phase_mag_fourier(a, 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
@@ -106,7 +106,7 @@ def run():
# Fourier padding 1:
padding = 1
start_time = time.clock()
- phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+ phase_num_disc = pm.phase_mag_fourier(a, projection_disc, padding=padding)
data_disc_fourier1[2, i] = time.clock() - start_time
print '------time (disc, fourier1) =', data_disc_fourier1[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
@@ -114,35 +114,35 @@ def run():
# Fourier padding 10:
padding = 10
start_time = time.clock()
- phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+ phase_num_disc = pm.phase_mag_fourier(a, 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:
start_time = time.clock()
- phase_num_disc = pm.phase_mag_real(res, projection_disc, 'slab')
+ phase_num_disc = pm.phase_mag_real(a, projection_disc, geometry='slab')
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:
start_time = time.clock()
- phase_num_disc = pm.phase_mag_real(res, projection_disc, 'disc')
+ phase_num_disc = pm.phase_mag_real(a, projection_disc, geometry='disc')
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'
+ print '----CALCULATE RMS/DURATION VORTEX STATE 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)
+ phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# Fourier unpadded:
padding = 0
start_time = time.clock()
- phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+ phase_num_vort = pm.phase_mag_fourier(a, 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
@@ -150,7 +150,7 @@ def run():
# Fourier padding 1:
padding = 1
start_time = time.clock()
- phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+ phase_num_vort = pm.phase_mag_fourier(a, projection_vort, padding=padding)
data_vort_fourier1[2, i] = time.clock() - start_time
print '------time (vortex, fourier1) =', data_vort_fourier1[2, i]
phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
@@ -158,21 +158,21 @@ def run():
# Fourier padding 10:
padding = 10
start_time = time.clock()
- phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+ phase_num_vort = pm.phase_mag_fourier(a, 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:
start_time = time.clock()
- phase_num_vort = pm.phase_mag_real(res, projection_vort, 'slab')
+ phase_num_vort = pm.phase_mag_real(a, projection_vort, geometry='slab')
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:
start_time = time.clock()
- phase_num_vort = pm.phase_mag_real(res, projection_vort, 'disc')
+ phase_num_vort = pm.phase_mag_real(a, projection_vort, geometry='disc')
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
@@ -192,20 +192,20 @@ def run():
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]:
+ # Plot duration against a (disc) [top/left]:
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_fourier1[0], data_disc_fourier1[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('res [nm]', fontsize=15)
+ axes[1, 0].set_xlabel('a [nm]', fontsize=15)
axes[1, 0].set_ylabel('RMS [mrad]', fontsize=15)
axes[1, 0].set_xlim(-0.5, 16.5)
axes[1, 0].tick_params(axis='both', which='major', labelsize=14)
axes[1, 0].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
- # Plot RMS against res (disc) [bottom/left]:
+ # Plot RMS against a (disc) [bottom/left]:
plt.tick_params(axis='both', which='major', labelsize=14)
axes[0, 0].set_yscale('log')
axes[0, 0].plot(data_disc_fourier0[0], data_disc_fourier0[2], ':bs')
@@ -219,19 +219,19 @@ def run():
axes[0, 0].tick_params(axis='both', which='major', labelsize=14)
axes[0, 0].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
- # Plot duration against res (vortex) [top/right]:
+ # Plot duration against a (vortex) [top/right]:
axes[1, 1].set_yscale('log')
axes[1, 1].plot(data_vort_fourier0[0], data_vort_fourier0[1], ':bs')
axes[1, 1].plot(data_vort_fourier1[0], data_vort_fourier1[1], ':bo')
axes[1, 1].plot(data_vort_fourier10[0], data_vort_fourier10[1], ':b^')
axes[1, 1].plot(data_vort_real_s[0], data_vort_real_s[1], '--rs')
axes[1, 1].plot(data_vort_real_d[0], data_vort_real_d[1], '--ro')
- axes[1, 1].set_xlabel('res [nm]', fontsize=15)
+ axes[1, 1].set_xlabel('a [nm]', fontsize=15)
axes[1, 1].set_xlim(-0.5, 16.5)
axes[1, 1].tick_params(axis='both', which='major', labelsize=14)
axes[1, 1].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
- # Plot RMS against res (vortex) [bottom/right]:
+ # Plot RMS against a (vortex) [bottom/right]:
axes[0, 1].set_yscale('log')
axes[0, 1].plot(data_vort_fourier0[0], data_vort_fourier0[2], ':bs',
label='Fourier padding=0')
diff --git a/scripts/paper 1/ch5-4-comparison_of_methods_new.py b/scripts/paper 1/ch5-4-comparison_of_methods_new.py
index b99f7646387ce0a346d9671bf1ab0b9b4e14f07a..6be62fcf7f172aaf5c53b626f1bdc66ec6a69775 100644
--- a/scripts/paper 1/ch5-4-comparison_of_methods_new.py
+++ b/scripts/paper 1/ch5-4-comparison_of_methods_new.py
@@ -24,6 +24,7 @@ import pyramid.projector as pj
import pyramid.phasemapper as pm
import pyramid.analytic as an
from pyramid.magdata import MagData
+from pyramid.kernel import Kernel
import shelve
@@ -50,7 +51,7 @@ def run():
else:
# Input parameters:
steps = 5
- res = 0.25 # in nm
+ a = 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!
@@ -61,9 +62,9 @@ def run():
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))
+ mag_data_disc = MagData(a, 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))
+ mag_data_vort = MagData(a, 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)]
@@ -79,26 +80,26 @@ def run():
data_vort_real_d = np.vstack((dim_list, np.zeros((2, steps))))
for i in range(steps):
- # Scale mag_data, resolution and dimensions:
+ # Scale mag_data, grid spacing and dimensions:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
dim = mag_data_disc.dim
- res = mag_data_disc.res
+ a = mag_data_disc.a
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 '--a =', a, '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)
+ phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
# Fourier unpadded:
padding = 0
start_time = time.clock()
- phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+ phase_num_disc = pm.phase_mag_fourier(a, 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
@@ -106,7 +107,7 @@ def run():
# Fourier padding 3:
padding = 3
start_time = time.clock()
- phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+ phase_num_disc = pm.phase_mag_fourier(a, 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
@@ -114,25 +115,23 @@ def run():
# Fourier padding 10:
padding = 10
start_time = time.clock()
- phase_num_disc = pm.phase_mag_fourier(res, projection_disc, padding=padding)
+ phase_num_disc = pm.phase_mag_fourier(a, 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)
+ kernel = Kernel((dim[1], dim[2]), a, geometry='slab')
start_time = time.clock()
- phase_num_disc = pm.phase_mag_real_fast(res, projection_disc, (v_kern_f, u_kern_f))
+ phase_num_disc = pm.phase_mag(a, projection_disc, kernel=kernel)
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)
+ kernel = Kernel((dim[1], dim[2]), a, geometry='slab')
start_time = time.clock()
- phase_num_disc = pm.phase_mag_real_fast(res, projection_disc, (v_kern_f, u_kern_f))
+ phase_num_disc = pm.phase_mag(a, projection_disc, kernel=kernel)
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
@@ -142,11 +141,11 @@ def run():
# 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)
+ phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# Fourier unpadded:
padding = 0
start_time = time.clock()
- phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+ phase_num_vort = pm.phase_mag_fourier(a, 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
@@ -154,7 +153,7 @@ def run():
# Fourier padding 3:
padding = 3
start_time = time.clock()
- phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+ phase_num_vort = pm.phase_mag_fourier(a, 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
@@ -162,25 +161,23 @@ def run():
# Fourier padding 10:
padding = 10
start_time = time.clock()
- phase_num_vort = pm.phase_mag_fourier(res, projection_vort, padding=padding)
+ phase_num_vort = pm.phase_mag_fourier(a, 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)
+ kernel = Kernel((dim[1], dim[2]), a, geometry='slab')
start_time = time.clock()
- phase_num_vort = pm.phase_mag_real_fast(res, projection_vort, (v_kern_f, u_kern_f))
+ phase_num_vort = pm.phase_mag(a, projection_vort, kernel=kernel)
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)
+ kernel = Kernel((dim[1], dim[2]), a, geometry='disc')
start_time = time.clock()
- phase_num_vort = pm.phase_mag_real_fast(res, projection_vort, (v_kern_f, u_kern_f))
+ phase_num_vort = pm.phase_mag(a, projection_vort, kernel=kernel)
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
@@ -200,7 +197,7 @@ def run():
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]:
+ # Plot duration against a (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')
@@ -217,7 +214,7 @@ def run():
axes[1, 0].xaxis.set_minor_locator(NullLocator())
axes[1, 0].grid()
- # Plot RMS against res (disc) [bottom/left]:
+ # Plot RMS against a (disc) [bottom/left]:
plt.tick_params(axis='both', which='major', labelsize=14)
axes[0, 0].set_xscale('log')
axes[0, 0].set_yscale('log')
@@ -235,7 +232,7 @@ def run():
axes[0, 0].xaxis.set_minor_locator(NullLocator())
axes[0, 0].grid()
- # Plot duration against res (vortex) [top/right]:
+ # Plot duration against a (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',
@@ -256,7 +253,7 @@ def run():
axes[1, 1].xaxis.set_minor_locator(NullLocator())
axes[1, 1].grid()
- # Plot RMS against res (vortex) [bottom/right]:
+ # Plot RMS against a (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',
diff --git a/scripts/reconstruction/get_jacobi.py b/scripts/reconstruction/get_jacobi.py
index 87a9b471ec3f0bceea2fbd7f03d4edfe345d6f15..c4919005a6a50178a8926915e18b4038ec879ed9 100644
--- a/scripts/reconstruction/get_jacobi.py
+++ b/scripts/reconstruction/get_jacobi.py
@@ -6,9 +6,6 @@ Created on Wed Apr 03 11:15:38 2013
"""
import time
-import pdb
-import traceback
-import sys
import os
import numpy as np
@@ -21,85 +18,65 @@ from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
-def get_jacobi():
- '''Calculate and display the phase map from a given magnetization.
- Arguments:
- None
- Returns:
- None
-
- '''
- directory = '../../output/reconstruction'
- if not os.path.exists(directory):
- os.makedirs(directory)
- filename = directory + '/jacobi.npy'
- # TODO: Input via GUI
- b_0 = 1.0 # in T
- dim = (1, 3, 3) # in px (y,x)
- res = 10.0 # in nm
- phi = pi/4
-
- center = (0, int(dim[1]/2), int(dim[2]/2)) # in px (y,x) index starts with 0!
- width = (0, 1, 1) # in px (y,x)
-
- mag_data = MagData(res, mc.create_mag_dist_homog(mc.Shapes.slab(dim, center, width), phi))
- projection = pj.simple_axis_projection(mag_data)
- print 'Projection calculated!'
-
- '''NUMERICAL SOLUTION'''
- # numerical solution Real Space:
- dim_proj = np.shape(projection[0])
- size = np.prod(dim_proj)
- kernel = pm.Kernel(dim_proj, res, 'disc')
-
- tic = time.clock()
- kernel.multiply_jacobi(np.ones(2*size)) # column-wise
- toc = time.clock()
- print 'Time for one multiplication with the Jacobi-Matrix: ', toc - tic
-
- jacobi = np.zeros((size, 2*size))
- tic = time.clock()
- phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc', b_0, jacobi=jacobi))
- toc = time.clock()
- phase_map.display()
- np.savetxt(filename, jacobi)
- print 'Time for Jacobi-Matrix during phase calculation: ', toc - tic
-
- tic = time.clock()
- jacobi_test = kernel.get_jacobi()
- toc = time.clock()
- print 'Time for Jacobi-Matrix from the Kernel: ', toc - tic
-
- unity = np.eye(2*size)
- jacobi_test2 = np.zeros((size, 2*size))
- tic = time.clock()
- for i in range(unity.shape[1]):
- jacobi_test2[:, i] = kernel.multiply_jacobi(unity[:, i]) # column-wise
- toc = time.clock()
- print 'Time for getting the Jacobi-Matrix (vector-wise): ', toc - tic
-
- unity_transp = np.eye(size)
- jacobi_transp = np.zeros((2*size, size))
- tic = time.clock()
- for i in range(unity_transp.shape[1]):
- jacobi_transp[:, i] = kernel.multiply_jacobi_T(unity_transp[:, i]) # column-wise
- toc = time.clock()
- print 'Time for getting the transp. Jacobi-Matrix (vector-wise):', toc - tic
-
- print 'Methods (during vs kernel) in accordance? ', \
- np.logical_not(np.all(jacobi-jacobi_test))
- print 'Methods (during vs vector-wise) in accordance? ', \
- np.logical_not(np.all(jacobi-jacobi_test2))
- print 'Methods (transponed Jacobi) in accordance? ', \
- np.logical_not(np.all(jacobi.T-jacobi_transp))
-
- return (jacobi, jacobi_test, jacobi_test2, jacobi_transp)
-
-
-if __name__ == "__main__":
- try:
- jacobi = get_jacobi()
- except:
- type, value, tb = sys.exc_info()
- traceback.print_exc()
- pdb.post_mortem(tb)
+directory = '../../output/reconstruction'
+if not os.path.exists(directory):
+ os.makedirs(directory)
+filename = directory + '/jacobi.npy'
+b_0 = 1.0 # in T
+dim = (1, 8, 8) # in px (y,x)
+res = 10.0 # in nm
+phi = pi/4
+
+center = (0, int(dim[1]/2), int(dim[2]/2)) # in px (y,x) index starts with 0!
+width = (0, 1, 1) # in px (y,x)
+
+mag_data = MagData(res, mc.create_mag_dist_homog(mc.Shapes.slab(dim, center, width), phi))
+projection = pj.simple_axis_projection(mag_data)
+print 'Projection calculated!'
+
+'''NUMERICAL SOLUTION'''
+# numerical solution Real Space:
+dim_proj = np.shape(projection[0])
+size = np.prod(dim_proj)
+kernel = pm.Kernel(dim_proj, res, 'disc')
+
+tic = time.clock()
+kernel.multiply_jacobi(np.ones(2*size)) # column-wise
+toc = time.clock()
+print 'Time for one multiplication with the Jacobi-Matrix: ', toc - tic
+
+jacobi = np.zeros((size, 2*size))
+tic = time.clock()
+phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, b_0, 'disc', jacobi=jacobi))
+toc = time.clock()
+phase_map.display()
+np.savetxt(filename, jacobi)
+print 'Time for Jacobi-Matrix during phase calculation: ', toc - tic
+
+tic = time.clock()
+jacobi_test = kernel.get_jacobi()
+toc = time.clock()
+print 'Time for Jacobi-Matrix from the Kernel: ', toc - tic
+
+unity = np.eye(2*size)
+jacobi_test2 = np.zeros((size, 2*size))
+tic = time.clock()
+for i in range(unity.shape[1]):
+ jacobi_test2[:, i] = kernel.multiply_jacobi(unity[:, i]) # column-wise
+toc = time.clock()
+print 'Time for getting the Jacobi-Matrix (vector-wise): ', toc - tic
+
+unity_transp = np.eye(size)
+jacobi_transp = np.zeros((2*size, size))
+tic = time.clock()
+for i in range(unity_transp.shape[1]):
+ jacobi_transp[:, i] = kernel.multiply_jacobi_T(unity_transp[:, i]) # column-wise
+toc = time.clock()
+print 'Time for getting the transp. Jacobi-Matrix (vector-wise):', toc - tic
+
+print 'Methods (during vs kernel) in accordance? ', \
+ np.logical_not(np.all(jacobi-jacobi_test))
+print 'Methods (during vs vector-wise) in accordance? ', \
+ np.logical_not(np.all(jacobi-jacobi_test2))
+print 'Methods (transponed Jacobi) in accordance? ', \
+ np.logical_not(np.all(jacobi.T-jacobi_transp))
diff --git a/scripts/reconstruction/reconstruct_core_shell.py b/scripts/reconstruction/reconstruct_core_shell.py
index 56d2c033dcf75c55b02766f547bc32424e5ccf09..91e09676bf4373784a1d3e13ebc042cdf8e69c5b 100644
--- a/scripts/reconstruction/reconstruct_core_shell.py
+++ b/scripts/reconstruction/reconstruct_core_shell.py
@@ -5,9 +5,6 @@ Created on Wed Nov 13 10:25:53 2013
@author: Jan
"""
-import sys
-import traceback
-import pdb
import os
import numpy as np
@@ -19,13 +16,12 @@ import pyramid.phasemapper as pm
import pyramid.projector as pj
import pyramid.holoimage as hi
from pyramid.magdata import MagData
+from pyramid.kernel import Kernel
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:
@@ -41,26 +37,27 @@ else:
mag_data = MagData.load_from_llg(directory +
'/magnetic distributions/mag_dist_core_shell_disc.txt')
mag_data.quiver_plot()
- res = mag_data.res
+ a = mag_data.a
dim = mag_data.dim
slices = 30
angles = np.linspace(0, 2*pi, slices, endpoint=True)
projections = []
phase_maps = []
+ kernel = Kernel((dim[1], dim[2]), a)
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_map = PhaseMap(a, pm.phase_mag(a, projection, kernel=kernel), 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')
+ 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)
@@ -87,7 +84,7 @@ 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)
+a = mag_data.a
+mag_data_reconstruct = MagData(a, (np.zeros(dim),)*3)
mag_vector = mag_data_reconstruct.get_vector(mask)
diff --git a/scripts/reconstruction/reconstruct_random_pixels.py b/scripts/reconstruction/reconstruct_random_pixels.py
index 710293e4da4457d88efb966e08ea1be12044f1d9..a3a1aedf8346794af8cb1f8feecbf1d159df1149 100644
--- a/scripts/reconstruction/reconstruct_random_pixels.py
+++ b/scripts/reconstruction/reconstruct_random_pixels.py
@@ -3,10 +3,6 @@
import random as rnd
-import pdb
-import traceback
-import sys
-import numpy as np
from numpy import pi
import pyramid.magcreator as mc
@@ -18,48 +14,31 @@ from pyramid.phasemap import PhaseMap
import pyramid.reconstructor as rc
-def reconstruct_random_distribution():
- '''Calculate and reconstruct a random magnetic distribution.
- Arguments:
- None
- Returns:
- None
-
- '''
- # Input parameters:
- n_pixel = 5
- dim = (1, 16, 16)
- b_0 = 1 # in T
- res = 10.0 # in nm
- rnd.seed(18)
-
- # Create empty MagData object and add random pixels:
- mag_data = MagData(res)
- for i in range(n_pixel):
- pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
- mag_shape = mc.Shapes.pixel(dim, pixel)
- phi = 2 * pi * rnd.random()
- magnitude = rnd.random()
- mag_data.add_magnitude(mc.create_mag_dist_homog(mag_shape, phi, magnitude))
- # Plot magnetic distribution, phase map and holographic contour map:
- mag_data.quiver_plot()
- projection = pj.simple_axis_projection(mag_data)
- phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab', b_0))
- hi.display_combined(phase_map, 10, 'Generated Distribution')
-
- # Reconstruct the magnetic distribution:
- mag_data_rec = rc.reconstruct_simple_leastsq(phase_map, mag_data.get_mask(), b_0)
-
- # Display the reconstructed phase map and holography image:
- projection_rec = pj.simple_axis_projection(mag_data_rec)
- phase_map_rec = PhaseMap(res, pm.phase_mag_real(res, projection_rec, 'slab', b_0))
- hi.display_combined(phase_map_rec, 10, 'Reconstructed Distribution')
-
-
-if __name__ == "__main__":
- try:
- reconstruct_random_distribution()
- except:
- type, value, tb = sys.exc_info()
- traceback.print_exc()
- pdb.post_mortem(tb)
+# Input parameters:
+n_pixel = 5
+dim = (1, 16, 16)
+b_0 = 1 # in T
+a = 10.0 # in nm
+rnd.seed(18)
+
+# Create empty MagData object and add random pixels:
+mag_data = MagData(a)
+for i in range(n_pixel):
+ pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
+ mag_shape = mc.Shapes.pixel(dim, pixel)
+ phi = 2 * pi * rnd.random()
+ magnitude = rnd.random()
+ mag_data.add_magnitude(mc.create_mag_dist_homog(mag_shape, phi, magnitude))
+# Plot magnetic distribution, phase map and holographic contour map:
+mag_data.quiver_plot()
+projection = pj.simple_axis_projection(mag_data)
+phase_map = PhaseMap(a, pm.phase_mag(a, projection, b_0))
+hi.display_combined(phase_map, 10, 'Generated Distribution')
+
+# Reconstruct the magnetic distribution:
+mag_data_rec = rc.reconstruct_simple_leastsq(phase_map, mag_data.get_mask(), b_0)
+
+# Display the reconstructed phase map and holography image:
+projection_rec = pj.simple_axis_projection(mag_data_rec)
+phase_map_rec = PhaseMap(a, pm.phase_mag(a, projection_rec, b_0))
+hi.display_combined(phase_map_rec, 10, 'Reconstructed Distribution')
diff --git a/scripts/test_arb_projection.py b/scripts/test_arb_projection.py
index fa71062af5674e0491d12f725fa110e1872c4bdd..c652331cc2d894f43f44c6a9ce95b6547050b994 100644
--- a/scripts/test_arb_projection.py
+++ b/scripts/test_arb_projection.py
@@ -5,7 +5,6 @@ Created on Tue Sep 03 12:55:40 2013
@author: Jan
"""
-'''2D Case Bresenham's line algorithm'''
import numpy as np
from numpy import pi
diff --git a/scripts/test_h5py.py b/scripts/test_h5py.py
index 40271b226546467672ddde6727f53aaf23ecad0e..7132dba75c8f4fa279a98a5aa9b6fdd4ed9f4a79 100644
--- a/scripts/test_h5py.py
+++ b/scripts/test_h5py.py
@@ -5,33 +5,23 @@ Created on Thu Nov 21 17:35:37 2013
@author: Jan
"""
+
import h5py
import numpy as np
-with h5py.File('testfile.hdf5', 'w') as f:
+with h5py.File('testfile.hdf5', 'w') as f:
dset = f.create_dataset('testdataset', (100,), dtype='i')
-
print 'dset.shape:', dset.shape
-
print 'dset.dtype:', dset.dtype
-
dset[...] = np.arange(100)
-
print 'dset[0]:', dset[0]
-
print 'dset[10]:', dset[10]
-
print 'dset[0:100:10]:', dset[0:100:10]
-
print 'dset.name:', dset.name
-
print 'f.name:', f.name
-
grp = f.create_group('subgroup')
-
dset_big = grp.create_dataset('another_set', (1000, 1000, 1000), dtype='f')
-
for i in range(dset_big.shape[0]):
print 'i:', i
dset_big[i, ...] = np.ones(dset_big.shape[1:])
diff --git a/scripts/test_methods.py b/scripts/test_methods.py
index f52fcfeeb770fc2b2bce131fcfd5828bf23a8c58..c7af2e6872061092cf1ba9e3daded5aad459916d 100644
--- a/scripts/test_methods.py
+++ b/scripts/test_methods.py
@@ -2,10 +2,6 @@
# -*- coding: utf-8 -*-
"""Compare the different methods to create phase maps."""
-import pdb
-import traceback
-import sys
-
from numpy import pi
import pyramid.magcreator as mc
@@ -16,72 +12,52 @@ from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
-def compare_methods():
- '''Calculate and display the phase map from a given magnetization.
- Arguments:
- None
- Returns:
- None
-
- '''
- # Input parameters:
- res = 10.0 # in nm
- phi = 0
- theta = pi/2
- tilt = pi/4
- density = 0.25
- dim = (128, 128, 128) # in px (z, y, x)
- # Create magnetic shape:
- 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]/3., dim[1]/4., dim[2]/2.) # in px (z, y, x)
- mag_shape = mc.Shapes.slab(dim, center, width)
- elif geometry == 'disc':
- 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
- mag_shape = mc.Shapes.disc(dim, center, radius, height)
- elif geometry == 'sphere':
- center = (dim[0]/2, dim[1]/2, dim[2]/2) # in px (z, y, x) index starts with 0!
- 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))
-
-# density = 0.3
-# mag_data = MagData.load_from_llg('../output/magnetic distributions/mag_dist_sphere.txt')
-
- res = mag_data.res
- mag_data.quiver_plot(ax_slice=dim[1]/2)
+# Input parameters:
+a = 10.0 # in nm
+phi = 0
+theta = pi/2
+tilt = pi/4
+density = 0.25
+dim = (128, 128, 128) # in px (z, y, x)
+# Create magnetic shape:
+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]/3., dim[1]/4., dim[2]/2.) # in px (z, y, x)
+ mag_shape = mc.Shapes.slab(dim, center, width)
+elif geometry == 'disc':
+ 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
+ mag_shape = mc.Shapes.disc(dim, center, radius, height)
+elif geometry == 'sphere':
+ center = (dim[0]/2, dim[1]/2, dim[2]/2) # in px (z, y, x) index starts with 0!
+ radius = dim[0]/4 # in px
+ mag_shape = mc.Shapes.sphere(dim, center, radius)
+# Project the magnetization data:
+mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, phi, theta))
+
+a = mag_data.a
+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)
- print 'Total projection time:', time.time() - start
- # Construct phase maps:
- phase_map_mag = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding=1))
- phase_map_elec = PhaseMap(res, pm.phase_elec(res, projection, v_0=3))
- # Display the combinated plots with phasemap and holography image:
- hi.display_combined(phase_map_mag, density, title='Magnetic Phase')
- hi.display_combined(phase_map_elec, density, title='Electric Phase')
-
- 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:
- compare_methods()
- except:
- type, value, tb = sys.exc_info()
- traceback.print_exc()
- pdb.post_mortem(tb)
+import time
+start = time.time()
+projection = pj.single_tilt_projection(mag_data, tilt)
+print 'Total projection time:', time.time() - start
+# Construct phase maps:
+phase_map_mag = PhaseMap(a, pm.phase_mag_fourier(a, projection, padding=1))
+phase_map_elec = PhaseMap(a, pm.phase_elec(a, projection, v_0=3))
+# Display the combinated plots with phasemap and holography image:
+hi.display_combined(phase_map_mag, density, title='Magnetic Phase')
+hi.display_combined(phase_map_elec, density, title='Electric Phase')
+
+phase_map = PhaseMap(a, 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()
diff --git a/scripts/test_vadim.py b/scripts/test_vadim.py
index ce1f3d5d9621ff73aedd566582c11653979ced15..229313c2f88afad9979393a1a3814ec526969a5a 100644
--- a/scripts/test_vadim.py
+++ b/scripts/test_vadim.py
@@ -19,6 +19,7 @@ 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
@@ -33,7 +34,7 @@ def calculate_charge_batch(phase_map):
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)
+ grad_y, grad_x = np.gradient(phase_map.phase, phase_map.a, phase_map.a)
xi, yi = np.zeros(t-b), np.zeros(t-b)
for i, ti in enumerate(np.arange(b, t)):
xi[i] = ti
@@ -44,10 +45,10 @@ directory = '../output/vadim/'
if not os.path.exists(directory):
os.makedirs(directory)
-res = 1.0 # in nm
+a = 1.0 # in nm
phi = pi/4
density = 30
-dim = (256, 256, 64) # in px (z, y, x)
+dim = (128, 128, 128) # in px (z, y, x)
l = dim[2]/4.
r = dim[2]/4. + dim[2]/2.
@@ -59,104 +60,104 @@ center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index
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))
+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))
mag_shape_ellipse = mc.Shapes.ellipsoid(dim, (center[0], 0*center[1], center[2]), (radius*5, radius*10, 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))
-mag_data_ellipse = MagData(res, mc.create_mag_dist_homog(mag_shape_ellipse, phi))
+mag_data_sphere = MagData(a, mc.create_mag_dist_homog(mag_shape_sphere, phi))
+mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape_disc, phi))
+mag_data_vortex = MagData(a, mc.create_mag_dist_vortex(mag_shape_disc))
+mag_data_slab = MagData(a, mc.create_mag_dist_homog(mag_shape_slab, phi))
+mag_data_ellipse = MagData(a, mc.create_mag_dist_homog(mag_shape_ellipse, 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)
+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)
projection_ellipse = pj.simple_axis_projection(mag_data_ellipse)
-## 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')
+# Magnetic phase map of a homogeneously magnetized disc:
+phase_map_mag_disc = PhaseMap(a, pm.phase_mag(a, 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(a, pm.phase_mag(a, 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(a, pm.phase_elec(a, 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(a, pm.phase_elec(a, 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(a, pm.phase_elec(a, 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')
# MIP phase of an ellipsoid:
-phase_map_mip_ellipsoid = PhaseMap(res, pm.phase_elec(res, projection_ellipse, v_0=1, v_acc=300000))
+phase_map_mip_ellipsoid = PhaseMap(a, pm.phase_elec(a, projection_ellipse, v_0=1, v_acc=300000))
phase_map_mip_ellipsoid.save_to_txt(directory+'phase_map_mip_ellipsoid.txt')
axis, _ = hi.display_combined(phase_map_mip_ellipsoid, density)
axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
@@ -171,35 +172,3 @@ np.savetxt(directory+'charge_integral_mip_ellipsoid.txt', np.array([x, y]).T)
plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mip_ellipsoid.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]))