Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • empyre/empyre
  • weber/empyre
  • wessels/empyre
  • bryan/empyre
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 2710 deletions
pyramid/__init__.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Package for the creation and reconstruction of magnetic distributions and resulting phase maps.
Modules
-------
magcreator
Create simple magnetic distributions.
magdata
Class for the storage of magnetizatin data.
projector
Create projections of a given magnetization distribution.
phasemapper
Create magnetic and electric phase maps from magnetization data.
phasemap
Class for the storage of phase data.
analytic
Create phase maps for magnetic distributions with analytic solutions.
holoimage
Create holographic contour maps from a given phase map.
reconstructor
Reconstruct magnetic distributions from given phasemaps.
Subpackages
-----------
numcore
Provides fast numerical functions for core routines.
"""
pyramid/analytic.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create phase maps for magnetic distributions with analytic solutions.
This module provides methods for the calculation of the magnetic phase for simple geometries for
which the analytic solutions are known. These can be used for comparison with the phase
calculated by the functions from the :mod:`~pyramid.phasemapper` module.
"""
import numpy as np
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):
'''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.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
The center of the slab in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the slab in pixel coordinates `(z, y, x)`.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phase_map : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
'''
# Function for the phase:
def phi_mag(x, y):
def F_0(x, y):
a = np.log(x**2 + y**2 + 1E-30)
b = np.arctan(x / (y+1E-30))
return x*a - 2*x + 2*y*b
return coeff * Lz * (- np.cos(phi) * (F_0(x-x0-Lx/2, y-y0-Ly/2)
- F_0(x-x0+Lx/2, y-y0-Ly/2)
- F_0(x-x0-Lx/2, y-y0+Ly/2)
+ F_0(x-x0+Lx/2, y-y0+Ly/2))
+ np.sin(phi) * (F_0(y-y0-Ly/2, x-x0-Lx/2)
- F_0(y-y0+Ly/2, x-x0-Lx/2)
- F_0(y-y0-Ly/2, x-x0+Lx/2)
+ 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]
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)
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):
'''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.
phi : float
The azimuthal angle, describing the direction of the magnetization.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`.
radius : float
The radius of the disc in pixel coordinates.
height : float
The height of the disc in pixel coordinates.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phase_map : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
'''
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
r[center[1], center[2]] = 1E-30
result = coeff * Lz * ((y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
result *= np.where(r <= R, 1, (R / r) ** 2)
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
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)
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):
'''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.
phi : float
The azimuthal angle, describing the direction of the magnetization.
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.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phase_map : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
'''
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
r[center[1], center[2]] = 1E-30
result = coeff * R ** 3 / r ** 2 * ((y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
result *= np.where(r > R, 1, (1 - (1 - (r / R) ** 2) ** (3. / 2.)))
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
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)
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):
'''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.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`, which is also the vortex center.
radius : float
The radius of the disc in pixel coordinates.
height : float
The height of the disc in pixel coordinates.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phase_map : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
'''
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
result = coeff * np.where(r <= R, r - R, 0)
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
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)
xx, yy = np.meshgrid(x, y)
# Return phase:
return phi_mag(xx, yy)
pyramid/holoimage.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create holographic contour maps from a given phase map.
This module converts phase maps into holographic contour map. This basically means taking the
cosine of the (optionally amplified) phase and encoding the direction of the 2-dimensional
gradient via color. The directional encoding can be seen by using the :func:`~.make_color_wheel`
function. Use the :func:`~.holoimage` function to create these holographic contour maps. It is
possible to use these as input for the :func:`~.display` to plot them, or just pass the
:class:`~pyramid.phasemap.PhaseMap` object to the :func:`~.display_combined` function to plot
the phase map and the holographic contour map next to each other.
"""
import numpy as np
from numpy import pi
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import NullLocator
from PIL import Image
from pyramid.phasemap import PhaseMap
CDICT = {'red': [(0.00, 1.0, 0.0),
(0.25, 1.0, 1.0),
(0.50, 1.0, 1.0),
(0.75, 0.0, 0.0),
(1.00, 0.0, 1.0)],
'green': [(0.00, 0.0, 0.0),
(0.25, 0.0, 0.0),
(0.50, 1.0, 1.0),
(0.75, 1.0, 1.0),
(1.00, 0.0, 1.0)],
'blue': [(0.00, 1.0, 1.0),
(0.25, 0.0, 0.0),
(0.50, 0.0, 0.0),
(0.75, 0.0, 0.0),
(1.00, 1.0, 1.0)]}
HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
def holo_image(phase_map, density=1):
'''Create a holographic contour map from a :class:`~pyramid.phasemap.PhaseMap` object.
Parameters
----------
phase_map : :class:`~pyramid.phasemap.PhaseMap`
A :class:`~pyramid.phasemap.PhaseMap` object storing the phase information.
density : float, optional
The gain factor for determining the number of contour lines. The default is 1.
Returns
-------
holo_image : :class:`~PIL.image`
The resulting holographic contour map with color encoded gradient.
'''
assert isinstance(phase_map, PhaseMap), 'phase_map has to be a PhaseMap object!'
# Extract the phase (considering the unit):
phase = phase_map.phase
# Calculate the holography image intensity:
img_holo = (1 + np.cos(density * phase)) / 2
# Calculate the phase gradients, expressed by magnitude and angle:
phase_grad_y, phase_grad_x = np.gradient(phase, phase_map.res, phase_map.res)
phase_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:
phase_magnitude = np.sin(phase_magnitude/phase_magnitude.max() * pi / 2)
# Color code the angle and create the holography image:
rgba = HOLO_CMAP(phase_angle)
rgb = (255.999 * img_holo.T * phase_magnitude.T * rgba[:, :, :3].T).T.astype(np.uint8)
holo_image = Image.fromarray(rgb)
return holo_image
def make_color_wheel():
'''Display a color wheel to illustrate the color coding of the gradient direction.
Parameters
----------
None
Returns
-------
None
'''
x = np.linspace(-256, 256, num=512)
y = np.linspace(-256, 256, num=512)
xx, yy = np.meshgrid(x, y)
r = np.sqrt(xx ** 2 + yy ** 2)
# Create the wheel:
color_wheel_magnitude = (1 - np.cos(r * pi/360)) / 2
color_wheel_magnitude *= 0 * (r > 256) + 1 * (r <= 256)
color_wheel_angle = (1 - np.arctan2(xx, -yy)/pi) / 2
# Color code the angle and create the holography image:
rgba = HOLO_CMAP(color_wheel_angle)
rgb = (255.999 * color_wheel_magnitude.T * rgba[:, :, :3].T).T.astype(np.uint8)
color_wheel = Image.fromarray(rgb)
# Plot the color wheel:
fig = plt.figure(figsize=(4, 4))
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.imshow(color_wheel, origin='lower')
axis.xaxis.set_major_locator(NullLocator())
axis.yaxis.set_major_locator(NullLocator())
def display(holo_image, title='Holographic Contour Map', axis=None, interpolation='none'):
'''Display the color coded holography image.
Parameters
----------
holo_image : :class:`~PIL.image`
The resulting holographic contour map with color encoded gradient.
title : string, optional
The title of the plot. The default is 'Holographic Contour Map'.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method. No interpolation is used in the default case.
Returns
-------
None
'''
# If no axis is specified, a new figure is created:
if axis is None:
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1, aspect='equal')
# Plot the image and set axes:
axis.imshow(holo_image, origin='lower', interpolation=interpolation)
# Set the title and the axes labels:
axis.set_title(title)
plt.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlabel('x-axis [px]', fontsize=15)
axis.set_ylabel('y-axis [px]', fontsize=15)
def display_combined(phase_map, density=1, title='Combined Plot', interpolation='none'):
'''Display a given phase map and the resulting color coded holography image in one plot.
Parameters
----------
phase_map : :class:`~pyramid.phasemap.PhaseMap`
A :class:`~pyramid.phasemap.PhaseMap` object storing the phase information.
density : float, optional
The gain factor for determining the number of contour lines. The default is 1.
title : string, optional
The title of the plot. The default is 'Combined Plot'.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
Returns
-------
None
'''
# Create combined plot and set title:
fig = plt.figure(figsize=(16, 7))
fig.suptitle(title, fontsize=20)
# Plot holography image:
holo_axis = fig.add_subplot(1, 2, 1, aspect='equal')
display(holo_image(phase_map, density), axis=holo_axis, interpolation=interpolation)
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect='equal')
fig.subplots_adjust(right=0.85)
phase_map.display(axis=phase_axis)
pyramid/magcreator.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create simple magnetic distributions.
The :mod:`~.magcreator` module is responsible for the creation of simple distributions of
magnetic moments. In the :class:`~.Shapes` class, you can find several general shapes for the
3-dimensional volume that should be magnetized (e.g. slabs, spheres, discs or single pixels).
These magnetic shapes are then used as input for the creating functions (or you could specify the
volume yourself as a 3-dimensional boolean matrix or a matrix with values in the range from 0 to 1,
which modifies the magnetization amplitude). The specified volume can either be magnetized
homogeneously with the :func:`~.create_mag_dist_homog` function by specifying the magnetization
direction, or in a vortex state with the :func:`~.create_mag_dist_vortex` by specifying the vortex
center.
"""
import numpy as np
from numpy import pi
class Shapes:
'''Class containing functions for generating magnetic shapes.
The :class:`~.Shapes` class is a collection of some mehtods 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.
'''
@classmethod
def slab(cls, dim, center, width):
'''Create the shape of a slab.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the slab in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the slab in pixel coordinates `(z, y, x)`.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a tuple of length 3!'
mag_shape = np.array([[[abs(x - center[2]) <= width[2] / 2
and abs(y - center[1]) <= width[1] / 2
and abs(z - center[0]) <= width[0] / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
return mag_shape
@classmethod
def disc(cls, dim, center, radius, height, axis='z'):
'''Create the shape of a zylindrical disc in x-, y-, or z-direction.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
The center of the disc in pixel coordinates `(z, y, x)`.
radius : float
The radius of the disc in pixel coordinates.
height : float
The height of the disc in pixel coordinates.
axis : {'z', 'y', 'x'}, optional
The orientation of the disc axis. The default is 'z'.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
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!'
assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as a string)!'
if axis == 'z':
mag_shape = np.array([[[np.hypot(x - center[2], y - center[1]) <= radius
and abs(z - center[0]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
elif axis == 'y':
mag_shape = np.array([[[np.hypot(x - center[2], z - center[0]) <= radius
and abs(y - center[1]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
elif axis == 'x':
mag_shape = np.array([[[np.hypot(y - center[1], z - center[0]) <= radius
and abs(x - center[2]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
return mag_shape
@classmethod
def sphere(cls, dim, center, radius):
'''Create the shape of a sphere.
Parameters
----------
dim : tuple (N=3)
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.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
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!'
mag_shape = np.array([[[np.sqrt((x-center[2])**2
+ (y-center[1])**2
+ (z-center[0])**2) <= radius
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
return mag_shape
@classmethod
def filament(cls, dim, pos, axis='y'):
'''Create the shape of a filament.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
pos : tuple (N=2)
The position of the filament in pixel coordinates `(coord1, coord2)`.
`coord1` and `coord2` stand for the two axes, which are perpendicular to `axis`. For
the default case (`axis = y`), it is `(coord1, coord2) = (z, x)`.
axis : {'y', 'x', 'z'}, optional
The orientation of the filament axis. The default is 'y'.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pos) == (2,), 'Parameter pos has to be a tuple of length 2!'
assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as a string)!'
mag_shape = np.zeros(dim)
if axis == 'z':
mag_shape[:, pos[0], pos[1]] = 1
elif axis == 'y':
mag_shape[pos[0], :, pos[1]] = 1
elif axis == 'x':
mag_shape[pos[0], pos[1], :] = 1
return mag_shape
@classmethod
def pixel(cls, dim, pixel):
'''Create the shape of a single pixel.
Parameters
----------
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
pixel : tuple (N=3)
The coordinates of the pixel `(z, y, x)`.
Returns
-------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shape as a 3D-array with values between 1 and 0.
'''
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pixel) == (3,), 'Parameter pixel has to be a tuple of length 3!'
mag_shape = np.zeros(dim)
mag_shape[pixel] = 1
return mag_shape
def create_mag_dist_homog(mag_shape, phi, theta=pi/2, magnitude=1):
'''Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
Parameters
----------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shapes (see Shapes.* for examples).
phi : float
The azimuthal angle, describing the direction of the magnetized object.
theta : float, optional
The polar angle, describing the direction of the magnetized object.
The default is pi/2, which means, the z-component is zero.
magnitude : float, optional
The relative magnitude for the magnetic shape. The default is 1.
Returns
-------
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The magnetic distribution as a tuple of the 3 components in
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
'''
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
z_mag = np.ones(dim) * np.cos(theta) * mag_shape * magnitude
y_mag = np.ones(dim) * np.sin(theta) * np.sin(phi) * mag_shape * magnitude
x_mag = np.ones(dim) * np.sin(theta) * np.cos(phi) * mag_shape * magnitude
return z_mag, y_mag, x_mag
def create_mag_dist_vortex(mag_shape, center=None, magnitude=1):
'''Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
Parameters
----------
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shapes (see :class:`.~Shapes` for examples).
center : tuple (N=2 or N=3), optional
The vortex center, given in 2D `(x, y)` or 3D `(x, y, z)`, where `z` is discarded.
Is set to the center of the field of view if not specified.
Vortex center has to be between two pixels.
magnitude : float, optional
The relative magnitude for the magnetic shape. The default is 1.
Returns
-------
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The magnetic distribution as a tuple of the 3 components in
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
'''
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
if center is None:
center = (int(dim[1]/2)-0.5, int(dim[2]/2)-0.5) # center has to be between (!) two pixels
if len(center) == 3: # if a 3D-center is given, just take the x and y components
center = (center[1], center[2])
assert len(center) == 2, 'Vortex center has to be defined in 3D or 2D or not at all!'
x = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
y = np.linspace(-center[0], dim[1]-1-center[0], dim[1])
xx, yy = np.meshgrid(x, y)
phi = np.arctan2(yy, xx) - pi/2
z_mag = np.zeros(dim)
y_mag = -np.ones(dim) * np.sin(phi) * mag_shape * magnitude
x_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
return z_mag, y_mag, x_mag
pyramid/magdata.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""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.
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from mayavi import mlab
import netCDF4
class MagData:
'''Class for storing magnetization data.
Represents 3-dimensional magnetic distributions with 3 components which are stored as a
tuple of size 3 in `magnitude`. The object can be created empty by omitting the `magnitude`
paramter. The `magnitude` can be added later by using the :func:`~.add_magnitude` function.
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.
Attributes
----------
res : float
The resolution of the grid (grid spacing in nm)
dim : tuple (N=3)
Dimensions of the grid.
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The `z`-, `y`- and `x`-component of the magnetization vector for every 3D-gridpoint
as a tuple.
'''
def __init__(self, res, magnitude=None):
'''Constructor for a :class:`~.MagData` object for storing magnetization data.
Parameters
----------
res : float
The resolution of the grid (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])
assert len(dim) == 3, 'Magnitude has to be defined for a 3-dimensional grid!'
assert np.shape(magnitude[1]) == np.shape(magnitude[2]) == dim, \
'Dimensions of the magnitude components do not match!'
self.magnitude = magnitude
self.dim = dim
else:
self.magnitude = None
self.dim = None
self.res = res
def add_magnitude(self, magnitude):
'''Add a given magnitude to the magnitude of the :class:`~.MagData`.
Parameters
----------
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
The `z`-, `y`- and `x`-component of the magnetization vector for every 3D-gridpoint
as a tuple. If the :class:`~.MagData` object already has a `magnitude`, the added
one has to match the dimensions `dim`, otherwise the added magnitude sets `dim`.
Returns
-------
None
'''
if self.magnitude is not None: # Add magnitude to existing one
assert np.shape(magnitude) == (3,) + self.dim, \
'Added magnitude has to have the same dimensions!'
z_mag, y_mag, x_mag = self.magnitude
z_new, y_new, x_new = magnitude
z_mag += z_new
y_mag += y_new
x_mag += x_new
self.magnitude = (z_mag, y_mag, x_mag)
else: # Magnitude is set for the first time and the dimensions are calculated
dim = np.shape(magnitude[0])
assert len(dim) == 3, 'Magnitude has to be defined for a 3-dimensional grid!'
assert np.shape(magnitude[1]) == np.shape(magnitude[2]) == dim, \
'Dimensions of the magnitude components do not match!'
self.magnitude = magnitude
self.dim = dim
def get_mask(self, threshold=0):
'''Mask all pixels where the amplitude of the magnetization lies above `threshold`.
Parameters
----------
threshold : float, optional
A pixel only gets masked, if it lies above this threshold . The default is 0.
Returns
-------
mask : :class:`~numpy.ndarray` (N=3, boolean)
Mask of the pixels where the amplitude of the magnetization lies above `threshold`.
'''
return np.sqrt(np.sum(np.array(self.magnitude)**2, axis=0)) > threshold
def get_vector(self, mask):
'''Returns the magnetic components arranged in a vector, specified by a mask.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean)
Masks the pixels from which the components should be taken.
Returns
-------
vector : :class:`~numpy.ndarray` (N=1)
The vector containing magnetization components of the specified pixels.
Order is: first all `x`-, then all `y`-, then all `z`-components.
'''
return np.concatenate([self.magnitude[2][mask],
self.magnitude[1][mask],
self.magnitude[0][mask]])
def set_vector(self, mask, vector):
'''Set the magnetic components of the masked pixels to the values specified by `vector`.
Parameters
----------
mask : :class:`~numpy.ndarray` (N=3, boolean)
Masks the pixels from which the components should be taken.
vector : :class:`~numpy.ndarray` (N=1)
The vector containing magnetization components of the specified pixels.
Order is: first all `x`-, then all `y-, then all `z`-components.
Returns
-------
None
'''
assert np.size(vector) % 3 == 0, 'Vector has to contain all 3 components for every pixel!'
count = np.size(vector)/3
self.magnitude[2][mask] = vector[:count] # x-component
self.magnitude[1][mask] = vector[count:2*count] # y-component
self.magnitude[0][mask] = vector[2*count:] # z-component
def scale_down(self, n=1):
'''Scale down the magnetic distribution by averaging over two pixels along each axis.
Parameters
----------
n : int, optional
Number of times the magnetic distribution is scaled down. The default is 1.
Returns
-------
None
Notes
-----
Acts in place and changes dimensions and resolution accordingly.
Only possible, if each axis length is a power of 2!
'''
assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
assert self.dim[0] % 2**n == 0 and self.dim[1] % 2**n == 0 and self.dim[2] % 2**n == 0, \
'For downscaling, every dimension must be a multiple of 2!'
for t in range(n):
# Create coarser grid for the magnetization:
z_mag = self.magnitude[0].reshape(self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2)
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.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'):
'''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'.
Returns
-------
None
'''
dim = self.dim
res = self.res * 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]
xx = xx.reshape(-1)
yy = yy.reshape(-1)
zz = zz.reshape(-1)
x_mag = np.reshape(self.magnitude[2], (-1))
y_mag = np.reshape(self.magnitude[1], (-1))
z_mag = np.reshape(self.magnitude[0], (-1))
# Save data to file:
data = np.array([xx, yy, zz, x_mag, y_mag, z_mag]).T
with open(filename, 'w') as mag_file:
mag_file.write('LLGFileCreator: %s\n' % filename.replace('.txt', ''))
mag_file.write(' %d %d %d\n' % (dim[2], dim[1], dim[0]))
mag_file.writelines('\n'.join(' '.join('{:7.6e}'.format(cell)
for cell in row) for row in data))
@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')
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))
def save_to_netcdf4(self, filename='..\output\magdata_output.nc'):
'''Save magnetization data in a file with NetCDF4-format.
Parameters
----------
filename : string, optional
The name of the NetCDF4-file in which to store the magnetization data.
The default is 'magdata_output.nc'.
Returns
-------
None
'''
mag_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
mag_file.res = self.res
mag_file.createDimension('z_dim', self.dim[0])
mag_file.createDimension('y_dim', self.dim[1])
mag_file.createDimension('x_dim', self.dim[2])
z_mag = mag_file.createVariable('z_mag', 'f', ('z_dim', 'y_dim', 'x_dim'))
y_mag = mag_file.createVariable('y_mag', 'f', ('z_dim', 'y_dim', 'x_dim'))
x_mag = mag_file.createVariable('x_mag', 'f', ('z_dim', 'y_dim', 'x_dim'))
z_mag[:] = self.magnitude[0]
y_mag[:] = self.magnitude[1]
x_mag[:] = self.magnitude[2]
print mag_file
mag_file.close()
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.
Parameters
----------
title : string, optional
The title for the plot.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
filename : string, optional
The filename, specifying the location where the image is saved. If not specified,
the image is shown instead.
proj_axis : {'z', 'y', 'x'}, optional
The axis, from which a slice is plotted. The default is 'z'.
ax_slice : int, optional
The slice-index of the axis specified in `proj_axis`. Is set to the center of
`proj_axis` if not specified.
Returns
-------
None
'''
assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
'Axis has to be x, y or z (as string).'
if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
if ax_slice is None:
ax_slice = int(self.dim[0]/2)
mag_slice_u = self.magnitude[2][ax_slice, ...]
mag_slice_v = self.magnitude[1][ax_slice, ...]
u_label = 'x [px]'
v_label = 'y [px]'
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
if ax_slice is None:
ax_slice = int(self.dim[1]/2)
mag_slice_u = self.magnitude[2][:, ax_slice, :]
mag_slice_v = self.magnitude[0][:, ax_slice, :]
u_label = 'x [px]'
v_label = 'z [px]'
elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
if ax_slice is None:
ax_slice = int(self.dim[2]/2)
mag_slice_u = self.magnitude[1][..., ax_slice]
mag_slice_v = self.magnitude[0][..., ax_slice]
u_label = 'y [px]'
v_label = 'z [px]'
# If no axis is specified, a new figure is created:
if axis is None:
fig = plt.figure(figsize=(8.5, 7))
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.quiver(mag_slice_u, mag_slice_v, pivot='middle', angles='xy', scale_units='xy',
scale=1, headwidth=6, headlength=7)
axis.set_xlim(-1, np.shape(mag_slice_u)[1])
axis.set_ylim(-1, np.shape(mag_slice_u)[0])
axis.set_title(title, fontsize=18)
axis.set_xlabel(u_label, fontsize=15)
axis.set_ylabel(v_label, fontsize=15)
axis.tick_params(axis='both', which='major', labelsize=14)
axis.xaxis.set_major_locator(MaxNLocator(nbins=8, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=8, integer=True))
plt.show()
def quiver_plot3d(self):
'''Plot the magnetization as 3D-vectors in a quiverplot.
Parameters
----------
None
Returns
-------
None
'''
res = self.res
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]
xx = xx.reshape(-1)
yy = yy.reshape(-1)
zz = zz.reshape(-1)
x_mag = np.reshape(self.magnitude[2], (-1))
y_mag = np.reshape(self.magnitude[1], (-1))
z_mag = np.reshape(self.magnitude[0], (-1))
# Plot them as vectors:
mlab.figure()
plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow', scale_factor=10.0)
mlab.outline(plot)
mlab.axes(plot)
mlab.colorbar()
pyramid/numcore/__init__.py deleted 100644 → 0
View file @ 6e9bc36e
"""Package for numerical core routines which enable fast computations.
Modules
-------
phase_mag_real
Provides numerical core routines for the real space approach in
:func:`~pyramid.phasemapper.phase_mag_real` of module :mod:`~pyramid.phasemapper`.
Notes
-----
Packages are written in `pyx`-format for the use with :mod:`~cython`.
"""
from phase_mag_real import *
pyramid/numcore/phase_mag_real.pyx deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Numerical core routines for the phase calculation using the real space approach.
Provides a helper function to speed up :func:`~pyramid.phasemapper.phase_mag_real` of module
:mod:`~pyramid.phasemapper`, by using C-speed for the for-loops and by omitting boundary and
wraparound checks.
"""
import numpy as np
import math
cimport cython
cimport numpy as np
@cython.boundscheck(False)
@cython.wraparound(False)
def phase_mag_real_core(
unsigned int v_dim, unsigned int u_dim,
double[:, :] v_phi, double[:, :] u_phi,
double[:, :] v_mag, double[:, :] u_mag,
double[:, :] phase, float threshold):
'''Numerical core routine for the phase calculation using the real space approach.
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.
v_mag, u_mag : :class:`~numpy.ndarray` (N=2)
Magnetization components in `u`- and `v`-direction.
phase : :class:`~numpy.ndarray` (N=2)
Matrix in which the resulting magnetic phase map should be stored.
threshold : float
The `threshold` determines which pixels contribute to the magnetic phase.
Returns
-------
None
'''
cdef unsigned int i, j, p, q, p_c, q_c
cdef double u_m, v_m
for j in range(v_dim):
for i in range(u_dim):
u_m = u_mag[j, i]
v_m = v_mag[j, i]
p_c = u_dim - 1 - i
q_c = v_dim - 1 - j
if abs(u_m) > threshold:
for q in range(v_dim):
for p in range(u_dim):
phase[q, p] += u_m * u_phi[q_c + q, p_c + p]
if abs(v_m) > threshold:
for q in range(v_dim):
for p in range(u_dim):
phase[q, p] -= v_m * v_phi[q_c + q, p_c + p]
pyramid/phasemap.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-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`
for further information.
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import netCDF4
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.
Attributes
----------
res : float
The resolution of the grid (grid spacing in nm)
dim : tuple (N=2)
Dimensions of the grid.
phase : :class:`~numpy.ndarray` (N=2)
Matrix containing the phase shift.
unit : {'rad', 'mrad', 'µrad'}, optional
Set the unit of the phase map. This is important for the :func:`~.display` function,
because the phase is scaled accordingly. Does not change the phase itself, which is
always in `rad`.
'''
UNITDICT = {'rad': 1E0,
'mrad': 1E3,
'µrad': 1E6}
def __init__(self, res, 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.
phase : :class:`~numpy.ndarray` (N=2)
Matrix containing the phase shift.
unit : {'rad', 'mrad', 'µrad'}, optional
Set the unit of the phase map. This is important for the :func:`~.display` function,
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.dim = dim
self.unit = unit
self.phase = phase
def set_unit(self, unit):
'''Set the unit for the phase map.
Parameters
----------
unit : {'rad', 'mrad'}, optional
Set the unit of the phase map. This is important for the :func:`~.display` function,
because the phase is scaled accordingly. Does not change the phase itself, which is
always in `rad`.
Returns
-------
None
'''
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.
Parameters
----------
filename : string
The name of the file in which to store the phase map data.
The default is 'phasemap_output.txt'.
Returns
-------
None
'''
with open(filename, 'w') as phase_file:
phase_file.write('{}\n'.format(filename.replace('.txt', '')))
phase_file.write('resolution = {} nm\n'.format(self.res))
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.
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')
res = phase_file.res
phase = phase_file.variables['phase'][:]
phase_file.close()
return PhaseMap(res, phase)
def save_to_netcdf4(self, filename='..\output\phasemap_output.nc'):
'''Save :class:`~.PhaseMap` data in a file with NetCDF4-format.
Parameters
----------
filename : string, optional
The name of the NetCDF4-file in which to store the phase data.
The default is 'phasemap_output.nc'.
Returns
-------
None
'''
phase_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
phase_file.res = self.res
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()
def display(self, title='Phase Map', cmap='RdBu', limit=None, norm=None, axis=None):
'''Display the phasemap as a colormesh.
Parameters
----------
title : string, optional
The title of the plot. The default is 'Phase Map'.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
The default is 'RdBu'.
limit : float, optional
Plotlimit for the phase in both negative and positive direction (symmetric around 0).
If not specified, the maximum amplitude the phase is used.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
Returns
-------
None
'''
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
if limit is None:
limit = np.max(np.abs(phase))
# If no axis is specified, a new figure is created:
if axis is None:
fig = plt.figure(figsize=(8.5, 7))
axis = fig.add_subplot(1, 1, 1, aspect='equal')
# Plot the phasemap:
im = axis.pcolormesh(phase, cmap=cmap, vmin=-limit, vmax=limit, norm=norm)
# Set the axes ticks and labels:
axis.set_xlim(0, np.shape(phase)[1])
axis.set_ylim(0, np.shape(phase)[0])
ticks = (axis.get_xticks()*self.res).astype(int)
axis.set_xticklabels(ticks)
ticks = (axis.get_yticks()*self.res).astype(int)
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_yticklabels(ticks)
axis.set_title(title, fontsize=18)
axis.set_xlabel('x [nm]', fontsize=15)
axis.set_ylabel('y [nm]', fontsize=15)
# Add colorbar:
fig = plt.gcf()
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.82, 0.15, 0.02, 0.7])
cbar = fig.colorbar(im, cax=cbar_ax)
cbar.ax.tick_params(labelsize=14)
cbar.set_label('phase shift [{}]'.format(self.unit), fontsize=15)
# Show plot:
plt.show()
def display3d(self, title='Phase Map', cmap='RdBu'):
'''Display the phasemap as a 3-D surface with contourplots.
Parameters
----------
title : string, optional
The title of the plot. The default is 'Phase Map'.
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
The default is 'RdBu'.
Returns
-------
None
'''
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
# Create figure and axis:
fig = plt.figure()
axis = Axes3D(fig)
# Plot surface and contours:
u = range(self.dim[1])
v = range(self.dim[0])
uu, vv = np.meshgrid(u, v)
axis.plot_surface(uu, vv, phase, rstride=4, cstride=4, alpha=0.7, cmap=cmap,
linewidth=0, antialiased=False)
axis.contourf(uu, vv, phase, 15, zdir='z', offset=np.min(phase), cmap=cmap)
axis.view_init(45, -135)
axis.set_xlabel('x-axis [px]')
axis.set_ylabel('y-axis [px]')
axis.set_zlabel('phase shift [{}]'.format(self.unit))
# Show Plot:
plt.show()
pyramid/phasemapper.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create magnetic and electric phase maps from magnetization data.
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).
"""
import numpy as np
from numpy import pi
import pyramid.numcore as nc
from scipy.signal import fftconvolve
PHI_0 = -2067.83 # magnetic flux in T*nm²
H_BAR = 6.626E-34 # Planck constant in J*s
M_E = 9.109E-31 # electron mass in kg
Q_E = 1.602E-19 # electron charge in C
C = 2.998E8 # speed of light in m/s
def phase_mag_fourier(res, projection, 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.
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.
padding : int, optional
Factor for the zero padding. The default is 0 (no padding). For a factor of n the number
of pixels is increase by ``(1+n)**2``. Padded zeros are cropped at the end.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
phase : :class:`~numpy.ndarray` (N=2)
The phase as a 2-dimensional array.
'''
v_dim, u_dim = np.shape(projection[0])
v_mag, u_mag = projection[:-1]
# Create zero padded matrices:
u_pad = u_dim/2 * padding
v_pad = v_dim/2 * padding
u_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
v_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
u_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = u_mag
v_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = v_mag
# Fourier transform of the two components:
u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_big), axes=0)
v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_big), axes=0)
# Calculate the Fourier transform of the phase:
nyq = 1 / res # nyquist frequency
f_u = np.linspace(0, nyq/2, u_mag_fft.shape[1])
f_v = np.linspace(-nyq/2, nyq/2, u_mag_fft.shape[0], endpoint=False)
f_uu, f_vv = np.meshgrid(f_u, f_v)
coeff = (1j*b_0) / (2*PHI_0)
phase_fft = coeff * res * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30)
# Transform to real space and revert padding:
phase_big = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
phase = phase_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
return phase
def phase_mag_real(res, projection, method='discs', b_0=1, jacobi=None):
'''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.
'''
# Function for creating the lookup-tables:
def phi_lookup(method, n, m, res, b_0):
if method == 'slab':
def F_h(n, m):
a = np.log(res**2 * (n**2 + m**2))
b = np.arctan(n / m)
return n*a - 2*n + 2*m*b
return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
- F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
elif method == 'disc':
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
# Process input parameters:
v_dim, u_dim = np.shape(projection[0])
v_mag, u_mag = projection[:-1]
coeff = -b_0 * res**2 / (2*PHI_0)
# Create lookup-tables for the phase of one pixel:
u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
uu, vv = np.meshgrid(u, v)
u_phi = phi_lookup(method, uu, vv, res, b_0)
v_phi = phi_lookup(method, vv, uu, res, b_0)
# Calculation of the phase:
phase = np.zeros((v_dim, u_dim))
threshold = 0
if jacobi is not None: # With Jacobian matrix (slower)
jacobi[:] = 0 # Jacobi matrix --> zeros
for j in range(v_dim):
for i in range(u_dim):
u_phase = u_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
jacobi[:, i+u_dim*j] = u_phase.reshape(-1)
if abs(u_mag[j, i]) > threshold:
phase += u_mag[j, i] * u_phase
v_phase = v_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
jacobi[:, u_dim*v_dim+i+u_dim*j] = -v_phase.reshape(-1)
if abs(v_mag[j, i]) > threshold:
phase -= v_mag[j, i] * v_phase
else: # Without Jacobi matrix (faster)
nc.phase_mag_real_core(v_dim, u_dim, v_phi, u_phi, v_mag, u_mag, phase, threshold)
# Return the phase:
return phase
def phase_mag_real_conv(res, projection, method='disc', b_0=1, jacobi=None):
'''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!!!
# Function for creating the lookup-tables:
def phi_lookup(method, n, m, res, b_0):
if method == 'slab':
def F_h(n, m):
a = np.log(res**2 * (n**2 + m**2))
b = np.arctan(n / m)
return n*a - 2*n + 2*m*b
return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
- F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
elif method == 'disc':
in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
# Process input parameters:
v_dim, u_dim = np.shape(projection[0])
v_mag, u_mag = projection[:-1]
coeff = -b_0 * res**2 / (2*PHI_0)
# Create lookup-tables for the phase of one pixel:
u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
uu, vv = np.meshgrid(u, v)
u_phi = phi_lookup(method, uu, vv, res, b_0)
v_phi = phi_lookup(method, vv, uu, res, b_0)
# Return the phase:
return fftconvolve(u_mag, u_phi, 'same') - fftconvolve(v_mag, v_phi, 'same')
# def fftconvolve(in1, in2, mode="full"):
# """Convolve two N-dimensional arrays using FFT.
#
# Convolve `in1` and `in2` using the fast Fourier transform method, with
# the output size determined by the `mode` argument.
#
# This is generally much faster than `convolve` for large arrays (n > ~500),
# but can be slower when only a few output values are needed, and can only
# output float arrays (int or object array inputs will be cast to float).
#
# Parameters
# ----------
# in1 : array_like
# First input.
# in2 : array_like
# Second input. Should have the same number of dimensions as `in1`;
# if sizes of `in1` and `in2` are not equal then `in1` has to be the
# larger array.
# mode : str {'full', 'valid', 'same'}, optional
# A string indicating the size of the output:
#
# ``full``
# The output is the full discrete linear convolution
# of the inputs. (Default)
# ``valid``
# The output consists only of those elements that do not
# rely on the zero-padding.
# ``same``
# The output is the same size as `in1`, centered
# with respect to the 'full' output.
#
# Returns
# -------
# out : array
# An N-dimensional array containing a subset of the discrete linear
# convolution of `in1` with `in2`.
#
# """
# in1 = asarray(in1)
# in2 = asarray(in2)
#
# if rank(in1) == rank(in2) == 0: # scalar inputs
# return in1 * in2
# elif not in1.ndim == in2.ndim:
# raise ValueError("in1 and in2 should have the same rank")
# elif in1.size == 0 or in2.size == 0: # empty arrays
# return array([])
#
# s1 = array(in1.shape)
# s2 = array(in2.shape)
# complex_result = (np.issubdtype(in1.dtype, np.complex) or
# np.issubdtype(in2.dtype, np.complex))
# size = s1 + s2 - 1
#
# if mode == "valid":
# for d1, d2 in zip(s1, s2):
# if not d1 >= d2:
# warnings.warn("in1 should have at least as many items as in2 in "
# "every dimension for 'valid' mode. In scipy "
# "0.13.0 this will raise an error",
# DeprecationWarning)
#
# # Always use 2**n-sized FFT
# fsize = 2 ** np.ceil(np.log2(size)).astype(int)
# print('fsize =', fsize)
# print('s1 =', s1)
# print('s2 =', s2)
# fslice = tuple([slice(0, int(sz)) for sz in size])
# if not complex_result:
# ret = irfftn(rfftn(in1, fsize) *
# rfftn(in2, fsize), fsize)[fslice].copy()
# ret = ret.real
# else:
# ret = ifftn(fftn(in1, fsize) * fftn(in2, fsize))[fslice].copy()
#
# if mode == "full":
# return ret
# elif mode == "same":
# return _centered(ret, s1)
# elif mode == "valid":
# return _centered(ret, abs(s1 - s2) + 1)
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
#def phase_mag_real2(res, projection, method, b_0=1, jacobi=None):
# '''Calculate phasemap from magnetization data (real space approach).
# Arguments:
# res - the resolution of the grid (grid spacing) in nm
# projection - projection of a magnetic distribution (created with pyramid.projector)
# method - string, describing the method to use for the pixel field ('slab' or 'disc')
# b_0 - magnetic induction corresponding to a magnetization Mo in T (default: 1)
# jacobi - matrix in which to save the jacobi matrix (default: None, faster computation)
# Returns:
# the phasemap as a 2 dimensional array
#
# '''
# # Function for creating the lookup-tables:
# def phi_lookup(method, n, m, res, b_0):
# if method == 'slab':
# def F_h(n, m):
# a = np.log(res**2 * (n**2 + m**2))
# b = np.arctan(n / m)
# return n*a - 2*n + 2*m*b
# return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
# - F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
# elif method == 'disc':
# in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
# return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
# # Process input parameters:
# v_dim, u_dim = np.shape(projection[0])
# v_mag, u_mag = projection
# coeff = -b_0 * res**2 / (2*PHI_0)
# # Create lookup-tables for the phase of one pixel:
# u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
# v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
# uu, vv = np.meshgrid(u, v)
# phi_u = phi_lookup(method, uu, vv, res, b_0)
# phi_v = phi_lookup(method, vv, uu, res, b_0)
# # Calculation of the phase:
# phase = np.zeros((v_dim, u_dim))
# threshold = 0
# if jacobi is not None: # With Jacobian matrix (slower)
# jacobi[:] = 0 # Jacobi matrix --> zeros
# ############################### TODO: NUMERICAL CORE ####################################
# for j in range(v_dim):
# for i in range(u_dim):
# phase_u = phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
# jacobi[:, i+u_dim*j] = phase_u.reshape(-1)
# if abs(u_mag[j, i]) > threshold:
# phase += u_mag[j, i] * phase_u
# phase_v = phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
# jacobi[:, u_dim*v_dim+i+u_dim*j] = -phase_v.reshape(-1)
# if abs(v_mag[j, i]) > threshold:
# phase -= v_mag[j, i] * phase_v
# ############################### TODO: NUMERICAL CORE ####################################
# else: # Without Jacobi matrix (faster)
## phasecopy = phase.copy()
## start_time = time.time()
## numcore.phase_mag_real_helper_1(v_dim, u_dim, phi_u, phi_v,
## u_mag, v_mag, phasecopy, threshold)
## print 'with numcore : ', time.time() - start_time
## start_time = time.time()
# for j in range(v_dim):
# for i in range(u_dim):
# if abs(u_mag[j, i]) > threshold:
# phase += u_mag[j, i] * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
# if abs(v_mag[j, i]) > threshold:
# phase -= v_mag[j, i] * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
## print 'without numcore: ', time.time() - start_time
## print 'num. difference: ', ((phase - phasecopy) ** 2).sum()
# # Return the phase:
# return phase
#
#
#def phase_mag_real_alt(res, projection, method, b_0=1, jacobi=None): # TODO: Modify
# '''Calculate phasemap from magnetization data (real space approach).
# Arguments:
# res - the resolution of the grid (grid spacing) in nm
# projection - projection of a magnetic distribution (created with pyramid.projector)
# method - string, describing the method to use for the pixel field ('slab' or 'disc')
# b_0 - magnetic induction corresponding to a magnetization Mo in T (default: 1)
# jacobi - matrix in which to save the jacobi matrix (default: None, faster computation)
# Returns:
# the phasemap as a 2 dimensional array
#
# '''
# # Function for creating the lookup-tables:
# def phi_lookup(method, n, m, res, b_0):
# if method == 'slab':
# def F_h(n, m):
# a = np.log(res**2 * (n**2 + m**2))
# b = np.arctan(n / m)
# return n*a - 2*n + 2*m*b
# return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
# - F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
# elif method == 'disc':
# in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
# return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
#
# # Function for the phase contribution of one pixel:
# def phi_mag(i, j):
# return (np.cos(beta[j, i]) * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
# - np.sin(beta[j, i]) * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
#
# # Function for the derivative of the phase contribution of one pixel:
# def phi_mag_deriv(i, j):
# return -(np.sin(beta[j, i]) * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
# + np.cos(beta[j, i]) * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
#
# # Process input parameters:
# v_dim, u_dim = np.shape(projection[0])
# v_mag, u_mag = projection
# beta = np.arctan2(v_mag, u_mag)
# mag = np.hypot(u_mag, v_mag)
# coeff = -b_0 * res**2 / (2*PHI_0)
# # Create lookup-tables for the phase of one pixel:
# u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
# v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
# uu, vv = np.meshgrid(u, v)
# phi_u = phi_lookup(method, uu, vv, res, b_0)
# phi_v = phi_lookup(method, vv, uu, res, b_0)
# # Calculation of the phase:
# phase = np.zeros((v_dim, u_dim))
# threshold = 0
# if jacobi is not None: # With Jacobian matrix (slower)
# jacobi[:] = 0 # Jacobi matrix --> zeros
# ############################### TODO: NUMERICAL CORE ####################################
# for j in range(v_dim):
# for i in range(u_dim):
# phase_cache = phi_mag(i, j)
# jacobi[:, i+u_dim*j] = phase_cache.reshape(-1)
# if mag[j, i] > threshold:
# phase += mag[j, i]*phase_cache
# jacobi[:, u_dim*v_dim+i+u_dim*j] = (mag[j, i]*phi_mag_deriv(i, j)).reshape(-1)
# ############################### TODO: NUMERICAL CORE ####################################
# else: # Without Jacobi matrix (faster)
# for j in range(v_dim):
# for i in range(u_dim):
# if abs(mag[j, i]) > threshold:
# phase += mag[j, i] * phi_mag(i, j)
# # Return the phase:
# return phase
pyramid/projector.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create projections of a given magnetization distribution.
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).
"""
import time
import numpy as np
from numpy import pi
import itertools
from pyramid.magdata import MagData
def simple_axis_projection(mag_data, axis='z', threshold=0):
'''
Project a magnetization distribution along one of the main axes of the 3D-grid.
Parameters
----------
mag_data : :class:`~pyramid.magdata.MagData`
A :class:`~pyramid.magdata.MagData` object storing the magnetization distribution,
which should be projected.
axis : {'z', 'y', 'x'}, optional
The projection direction as a string.
threshold : float, optional
A pixel only gets masked, if it lies above this threshold. The default is 0.
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 resolution 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)!'
if axis == 'z':
projection = (mag_data.magnitude[1].sum(0), # y_mag -> v_mag
mag_data.magnitude[2].sum(0), # x_mag -> u_mag
mag_data.get_mask(threshold).sum(0)) # thickness profile
elif axis == 'y':
projection = (mag_data.magnitude[0].sum(1), # z_mag -> v_mag
mag_data.magnitude[2].sum(1), # x_mag -> u_mag
mag_data.get_mask(threshold).sum(1)) # thickness profile
elif axis == 'x':
projection = (mag_data.magnitude[0].sum(2), # z_mag -> v_mag
mag_data.magnitude[1].sum(2), # y_mag -> u_mag
mag_data.get_mask(threshold).sum(2)) # thickness profile
return projection
def single_tilt_projection(mag_data, tilt=0, threshold=0):
# TODO: Docstring!!!
assert isinstance(mag_data, MagData), 'Parameter mag_data has to be a MagData object!'
# assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as string)!'
# Set starting variables:
dim = (mag_data.dim[0], mag_data.dim[2])
z_mag, y_mag, x_mag = mag_data.magnitude
mask = mag_data.get_mask()
projection = (np.zeros((mag_data.dim[1], mag_data.dim[2])),
np.zeros((mag_data.dim[1], mag_data.dim[2])),
np.zeros((mag_data.dim[1], mag_data.dim[2])))
def get_position(p, m, b, size):
x, y = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
def get_impact(pos, r, size):
return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
if 0 <= x < size]
def get_weight(delta, rho):
a, b = delta-rho, delta+rho
if a >= 1 or b <= -1: # TODO: Should not be necessary!
print 'Outside of bounds:', delta
return 0
# Upper boundary:
if b >= 1:
w_b = 0.5
else:
w_b = (b*np.sqrt(1-b**2) + np.arctan(b/np.sqrt(1-b**2))) / pi
# Lower boundary:
if a <= -1:
w_a = -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
# Creating coordinate list of all voxels:
xi = range(dim[0])
yj = range(dim[1])
ii, jj = np.meshgrid(xi, yj)
voxels = list(itertools.product(yj, xi))
# Calculate positions along the projected pixel coordinate system:
direction = (-np.cos(tilt), np.sin(tilt))
center = (dim[0]/2., dim[1]/2.)
m = direction[0]/(direction[1]+1E-30)
b = center[0] - m * center[1]
positions = get_position(voxels, m, b, dim[0])
# Calculate weights:
r = 1/np.sqrt(np.pi) # radius of the voxel circle
rho = 0.5 / r # TODO: ratio of radii
weights = []
for i, voxel in enumerate(voxels):
voxel_weights = []
impacts = get_impact(positions[i], r, dim[0])
for impact in impacts:
distance = np.abs(impact+0.5 - positions[i])
delta = distance / r
voxel_weights.append((impact, get_weight(delta, rho)))
weights.append((voxel, voxel_weights))
# Calculate projection with the calculated weights for the voxels:
for i, weight in enumerate(weights):
voxel = weights[i][0]
voxel_weights = weights[i][1]
for voxel_weight in voxel_weights:
pixel, weight = voxel_weight
# 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]]
return projection
pyramid/reconstructor.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Reconstruct magnetic distributions from given phasemaps.
This module reconstructs 3-dimensional magnetic distributions (as :class:`~pyramid.magdata.MagData`
objects) from a given set of phase maps (represented by :class:`~pyramid.phasemap.PhaseMap`
objects) by using several model based reconstruction algorithms which use the forward model
provided by :mod:`~pyramid.projector` and :mod:`~pyramid.phasemapper` and a priori knowledge of
the distribution.
So far, only a simple least square algorithm for known pixel locations for 2-dimensional problems
is implemented (:func:`~.reconstruct_simple_leastsq`), but more complex solutions are planned.
"""
import numpy as np
from scipy.optimize import leastsq
import pyramid.projector as pj
import pyramid.phasemapper as pm
from pyramid.magdata import MagData
def reconstruct_simple_leastsq(phase_map, mask, b_0=1):
'''Reconstruct a magnetic distribution for a 2-D problem with known pixel locations.
Parameters
----------
phase_map : :class:`~pyramid.phasemap.PhaseMap`
A :class:`~pyramid.phasemap.PhaseMap` object, representing the phase from which to
reconstruct the magnetic distribution.
mask : :class:`~numpy.ndarray` (N=3)
A boolean matrix (or a matrix consisting of ones and zeros), representing the
positions of the magnetized voxels in 3 dimensions.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
Returns
-------
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a
:class:`~pyramid.magdata.MagData` object.
Notes
-----
Only works for a single phase_map, if the positions of the magnetized voxels are known and
for slice thickness of 1 (constraint for the `z`-dimension).
'''
# Read in parameters:
y_m = phase_map.phase.reshape(-1) # Measured phase map as a vector
res = phase_map.res # Resolution
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)))
# 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)
return phase.reshape(-1)
# Cost function which should be minimized:
def J(x_i):
y_i = F(x_i)
term1 = (y_i - y_m)
term2 = lam * x_i
return np.concatenate([term1, term2])
# Reconstruct the magnetization components:
x_rec, _ = leastsq(J, np.zeros(3*count))
mag_data_rec.set_vector(mask, x_rec)
return mag_data_rec
scripts/compare methods/compare_methods.py deleted 100644 → 0
View file @ 6e9bc36e
#! python
# -*- coding: utf-8 -*-
"""Compare the different methods to create phase maps."""
import time
import pdb
import traceback
import sys
import numpy as np
from numpy import pi
import pyramid.magcreator as mc
import pyramid.projector as pj
import pyramid.phasemapper as pm
import pyramid.holoimage as hi
import pyramid.analytic as an
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:
b_0 = 1.1 # in T
res = 10.0 # in nm
phi = pi/4
padding = 3
density = 1
dim = (8, 512, 512) # in px (z, y, x)
# Create magnetic shape:
geometry = 'slab'
if geometry == 'slab':
center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z,y,x) index starts at 0!
width = (dim[0]/2, dim[1]/2, dim[2]/2) # in px (z, y, x)
mag_shape = mc.Shapes.slab(dim, center, width)
phase_ana = an.phase_mag_slab(dim, res, 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)
elif geometry == 'sphere':
center = (50, 50, 50) # in px (z, y, x) index starts with 0!
radius = 25 # in px
mag_shape = mc.Shapes.sphere(dim, center, radius)
phase_ana = an.phase_mag_sphere(dim, res, phi, center, radius, b_0)
# Project the magnetization data:
mag_data = MagData(res, mc.create_mag_dist_homog(mag_shape, phi))
mag_data.quiver_plot(ax_slice=int(center[0]))
projection = pj.simple_axis_projection(mag_data)
# Construct phase maps:
phase_map_ana = PhaseMap(res, phase_ana)
start_time = time.time()
phase_map_fft = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding, b_0))
print 'Time for Fourier space approach: ', time.time() - start_time
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
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
start_time = time.time()
phase_map_slab_conv = PhaseMap(res, pm.phase_mag_real_conv(res, projection, 'slab', b_0))
print 'Time for real space approach (Slab Convolve):', time.time() - start_time
start_time = time.time()
phase_map_disc_conv = PhaseMap(res, pm.phase_mag_real_conv(res, projection, 'disc', b_0))
print 'Time for real space approach (Disc Convolve):', time.time() - start_time
# 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')
hi.display_combined(phase_map_slab, density, 'Real Space (Slab)')
hi.display_combined(phase_map_disc, density, 'Real Space (Disc)')
hi.display_combined(phase_map_slab_conv, density, 'Real Space (Slab Convolve)')
hi.display_combined(phase_map_disc_conv, density, 'Real Space (Disc Convolve)')
# 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)
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)
RMS_slab_conv = np.std(phase_map_diff_slab_conv.phase)
RMS_disc_conv = np.std(phase_map_diff_disc_conv.phase)
phase_map_diff_fft.display('Fourier Space (RMS = {:3.2e})'.format(RMS_fft))
phase_map_diff_slab.display('Real Space (Slab) (RMS = {:3.2e})'.format(RMS_slab))
phase_map_diff_disc.display('Real Space (Disc) (RMS = {:3.2e})'.format(RMS_disc))
phase_map_diff_slab_conv.display('Real Space (Disc) (RMS = {:3.2e})'.format(RMS_slab_conv))
phase_map_diff_disc_conv.display('Real Space (Disc) (RMS = {:3.2e})'.format(RMS_disc_conv))
if __name__ == "__main__":
try:
compare_methods()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
scripts/compare methods/compare_pixel_fields.py deleted 100644 → 0
View file @ 6e9bc36e
#! python
# -*- coding: utf-8 -*-
"""Compare the phase map of one pixel for different real space approaches."""
import pdb
import traceback
import sys
import os
import numpy as np
import pyramid.magcreator as mc
import pyramid.projector as pj
import pyramid.phasemapper as pm
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
def compare_pixel_fields():
'''Calculate and display the phase map for different real space approaches.
Arguments:
None
Returns:
None
'''
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
res = 1.0 # in nm
phi = 0 # in rad
dim = (1, 32, 32)
pixel = (0, int(dim[1]/2), int(dim[2]/2))
limit = 0.35
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
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)
# 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.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.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.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.display('Phase of one Pixel (Fourier)', limit=limit)
# Kernel of the Fourier method, calculated directly:
phase_map_fft_kernel = PhaseMap(res, get_fourier_kernel(), 'mrad')
phase_map_fft_kernel.display('Phase of one Pixel (Fourier Kernel)', limit=limit)
# Kernel differences:
print 'Fourier Kernel, direct and indirect method are identical:', \
np.all(phase_map_fft_kernel.phase - phase_map_fft.phase) == 0
phase_map_diff = PhaseMap(res, phase_map_fft.phase - phase_map_disc.phase, 'mrad')
phase_map_diff.display('Phase difference of one Pixel (Disc - Fourier)')
phase_inv_fft_r = np.real(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))
phase_map_inv_fft = PhaseMap(res, phase_inv_fft_r)
phase_map_inv_fft.display('FT of the Phase of one Pixel (Fourier, real)')
phase_inv_fft_i = np.imag(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))
phase_map_inv_fft = PhaseMap(res, phase_inv_fft_i)
phase_map_inv_fft.display('FT of the Phase of one Pixel (Fourier, imag)')
phase_inv_disc_r = np.real(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))
phase_map_inv_disc = PhaseMap(res, phase_inv_disc_r)
phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, real)')
phase_inv_disc_i = np.imag(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))
phase_map_inv_disc = PhaseMap(res, phase_inv_disc_i)
phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, imag)')
if __name__ == "__main__":
try:
compare_pixel_fields()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
scripts/create distributions/create_alternating_filaments.py deleted 100644 → 0
View file @ 6e9bc36e
#! python
# -*- coding: utf-8 -*-
"""Create magnetic distribution of alternating filaments"""
import pdb
import traceback
import sys
import os
from numpy import pi
import pyramid.magcreator as mc
from pyramid.magdata import MagData
def create_alternating_filaments():
'''Calculate, display and save the magnetic distribution of alternating filaments to file.
Arguments:
None
Returns:
None
'''
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_alt_filaments.txt'
dim = (1, 21, 21) # in px (z, y, x)
res = 10.0 # in nm
phi = pi/2
spacing = 5
# Create empty MagData object:
mag_data = MagData(res)
count = int((dim[1]-1) / spacing) + 1
for i in range(count):
pos = i * spacing
mag_shape = mc.Shapes.filament(dim, (0, pos))
mag_data.add_magnitude(mc.create_mag_dist_homog(mag_shape, phi))
phi *= -1 # Switch the angle
# Plot magnetic distribution
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
# print np.max(mag_data.magnitude)
if __name__ == "__main__":
try:
create_alternating_filaments()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
scripts/create distributions/create_logo.py deleted 100644 → 0
View file @ 6e9bc36e
#! python
# -*- coding: utf-8 -*-
"""Create the Pyramid-Logo."""
import pdb
import traceback
import sys
import numpy as np
from numpy import pi
import pyramid.magcreator as mc
import pyramid.projector as pj
import pyramid.phasemapper as pm
import pyramid.holoimage as hi
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
def create_logo():
'''Calculate and display the Pyramid-Logo.
Arguments:
None
Returns:
None
'''
# Input parameters:
res = 10.0 # in nm
phi = -pi/2 # in rad
density = 10
dim = (1, 128, 128)
# Create magnetic shape:
mag_shape = np.zeros(dim)
x = range(dim[2])
y = range(dim[1])
xx, yy = np.meshgrid(x, y)
bottom = (yy >= 0.25*dim[1])
left = (yy <= 0.75/0.5 * dim[1]/dim[2] * xx)
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.quiver_plot()
projection = pj.simple_axis_projection(mag_data)
phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
hi.display(hi.holo_image(phase_map, density), 'PYRAMID - LOGO', interpolation='bilinear')
if __name__ == "__main__":
try:
create_logo()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
scripts/create distributions/create_multiple_samples.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create random magnetic distributions."""
import pdb
import traceback
import sys
import os
from numpy import pi
import pyramid.magcreator as mc
from pyramid.magdata import MagData
import pyramid.phasemapper as pm
import pyramid.projector as pj
import pyramid.holoimage as hi
from pyramid.phasemap import PhaseMap
def create_multiple_samples():
'''Calculate, display and save a random magnetic distribution to file.
Arguments:
None
Returns:
None
'''
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_multiple_samples.txt'
res = 10.0 # nm
dim = (64, 128, 128)
# Slab:
center = (32, 32, 32) # in px (z, y, x), index starts with 0!
width = (48, 48, 48) # in px (z, y, x)
mag_shape_slab = mc.Shapes.slab(dim, center, width)
# Disc:
center = (32, 32, 96) # in px (z, y, x), index starts with 0!
radius = 24 # in px
height = 24 # in px
mag_shape_disc = mc.Shapes.disc(dim, center, radius, height)
# Sphere:
center = (32, 96, 64) # in px (z, y, x), index starts with 0!
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.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))
# Plot the 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.display()
hi.display(hi.holo_image(phase_map, 0.5))
if __name__ == "__main__":
try:
create_multiple_samples()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
scripts/create distributions/create_random_pixels.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create random magnetic distributions."""
import random as rnd
import pdb
import traceback
import sys
import os
import numpy as np
from numpy import pi
import pyramid.magcreator as mc
from pyramid.magdata import MagData
import pyramid.phasemapper as pm
import pyramid.projector as pj
import pyramid.holoimage as hi
from pyramid.phasemap import PhaseMap
def create_random_pixels():
'''Calculate, display and save a random magnetic distribution to file.
Arguments:
None
Returns:
None
'''
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_random_pixels.txt'
# Input parameters:
count = 10
dim = (1, 128, 128)
res = 10 # in nm
rnd.seed(12)
# Create empty MagData object and add pixels:
mag_data = MagData(res)
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()
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'))
hi.display(hi.holo_image(phase_map, 10))
if __name__ == "__main__":
try:
create_random_pixels()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
scripts/create distributions/create_random_slabs.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create random magnetic distributions."""
import random as rnd
import pdb
import traceback
import sys
import os
import numpy as np
from numpy import pi
import pyramid.magcreator as mc
from pyramid.magdata import MagData
import pyramid.phasemapper as pm
import pyramid.projector as pj
import pyramid.holoimage as hi
from pyramid.phasemap import PhaseMap
def create_random_slabs():
'''Calculate, display and save a random magnetic distribution to file.
Arguments:
None
Returns:
None
'''
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_random_pixels.txt'
count = 10
dim = (1, 128, 128)
res = 10 # in nm
rnd.seed(42)
w_max = 10
# Create empty MagData object and add slabs:
mag_data = MagData(res)
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)),
rnd.randrange(int(width[1]/2), dim[1]-int(width[1]/2)),
rnd.randrange(int(width[2]/2), dim[2]-int(width[2]/2)))
mag_shape = mc.Shapes.slab(dim, center, width)
phi = 2 * pi * rnd.random()
magnitude = 1 # TODO: rnd.random()
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'))
hi.display(hi.holo_image(phase_map, 10))
if __name__ == "__main__":
try:
create_random_slabs()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
scripts/create distributions/create_sample.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create magnetic distributions with simple geometries."""
import pdb
import traceback
import sys
import os
from numpy import pi
import pyramid.magcreator as mc
from pyramid.magdata import MagData
def create_sample():
'''Calculate, display and save simple magnetic distributions to file.
Arguments:
None
Returns:
None
'''
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
key = 'sphere'
filename = directory + '/mag_dist_' + key + '.txt'
dim = (128, 128, 128) # in px (z, y, x)
res = 10.0 # in nm
phi = pi/4
# Geometry parameters:
center = (64, 64, 64) # in px (z, y, x), index starts with 0!
width = (1, 50, 50) # in px (z, y, x)
radius = 25 # in px
height = 1 # in px
pos = (0, 63) # in px (tuple of length 2)
pixel = (0, 63, 63) # in px (z, y, x), index starts with 0!
# Determine the magnetic shape:
if key == 'slab':
mag_shape = mc.Shapes.slab(dim, center, width)
elif key == 'disc':
mag_shape = mc.Shapes.disc(dim, center, radius, height)
elif key == 'sphere':
mag_shape = mc.Shapes.sphere(dim, center, radius)
elif key == 'filament':
mag_shape = mc.Shapes.filament(dim, pos)
elif key == 'pixel':
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.quiver_plot()
mag_data.save_to_llg(filename)
if __name__ == "__main__":
try:
create_sample()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)
scripts/create distributions/create_vortex.py deleted 100644 → 0
View file @ 6e9bc36e
# -*- coding: utf-8 -*-
"""Create the Pyramid-Logo."""
import pdb
import traceback
import sys
import os
import matplotlib.pyplot as plt
import pyramid.magcreator as mc
import pyramid.projector as pj
import pyramid.phasemapper as pm
import pyramid.holoimage as hi
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
def create_vortex():
'''Calculate and display the Pyramid-Logo.
Arguments:
None
Returns:
None
'''
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_vortex.txt'
res = 10.0 # in nm
density = 1
dim = (1, 128, 128)
center = (0, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
radius = 0.25 * dim[1]
height = 1
# 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.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'))
hi.display_combined(phase_map, density, 'Vortex State')
phase_slice = phase_map.phase[dim[1]/2, :]
fig = plt.figure()
fig.add_subplot(111)
plt.plot(range(dim[1]), phase_slice)
if __name__ == "__main__":
try:
create_vortex()
except:
type, value, tb = sys.exc_info()
traceback.print_exc()
pdb.post_mortem(tb)