pyramid/magcreator.py

-# -*- coding: utf-8 -*-
-"""Create simple LLG Files which describe magnetization in 2D (z-Dim=1)."""
-
-
-import numpy as np
-from numpy import pi
-
-
-class Shapes:
-
-    '''Class containing functions for generating shapes.'''
-
-    @classmethod
-    def soft_slab(cls, dim, center, width):
-        '''Get the magnetic shape of a slab.
-        Arguments:
-            dim    - the dimensions of the grid, shape(z, y, x)
-            center - the center of the slab in pixel coordinates, shape(z, y, x)
-            width  - the width of the slab in pixel coordinates, shape(z, y, x)
-        Returns:
-            the magnetic shape as a 3D-array with ones and zeros
-
-        '''
-        assert np.shape(dim) == (3,), 'Parameter dim has to be a shape of 3 dimensions!'
-        assert np.shape(center) == (3,), 'Parameter center has to be a shape of 3 dimensions!'
-        assert np.shape(width) == (3,), 'Parameter width has to be a shape of 3 dimensions!'
-        z0, y0, x0 = center
-        Lz, Ly, Lx = width
-        mag_shape = np.zeros(dim)
-        for z in range(dim[0]):
-            for y in range(dim[1]):
-                for x in range(dim[2]):
-                    mag_shape[z, y, x] = (max(min(x+0.5, x0+Lx/2.)-max(x-0.5, x0-Lx/2.), 0)
-                                        * max(min(y+0.5, y0+Ly/2.)-max(y-0.5, y0-Ly/2.), 0)
-                                        * max(min(z+0.5, z0+Lz/2.)-max(z-0.5, z0-Lz/2.), 0))
-        return mag_shape
-
-    @classmethod
-    def slab(cls, dim, center, width):
-        '''Get the magnetic shape of a slab.
-        Arguments:
-            dim    - the dimensions of the grid, shape(z, y, x)
-            center - the center of the slab in pixel coordinates, shape(z, y, x)
-            width  - the width of the slab in pixel coordinates, shape(z, y, x)
-        Returns:
-            the magnetic shape as a 3D-array with ones and zeros
-
-        '''
-        assert np.shape(dim) == (3,), 'Parameter dim has to be a shape of 3 dimensions!'
-        assert np.shape(center) == (3,), 'Parameter center has to be a shape of 3 dimensions!'
-        assert np.shape(width) == (3,), 'Parameter width has to be a shape of 3 dimensions!'
-        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'):
-        '''Get the magnetic shape of a zylindrical disc in x-, y-, or z-direction.
-        Arguments:
-            dim    - the dimensions of the grid, shape(z, y, x)
-            center - the center of the disc in pixel coordinates, shape(z, y, x)
-            radius - the radius of the disc in pixel coordinates (scalar value)
-            height - the height of the disc in pixel coordinates (scalar value)
-            axis   - the orientation of the disc axis, (String: 'x', 'y', 'z'), default = 'z'
-        Returns:
-            the magnetic shape as a 3D-array with ones and zeros
-
-        '''
-        assert np.shape(dim) == (3,), 'Parameter dim has to be a shape of 3 dimensions!'
-        assert np.shape(center) == (3,), 'Parameter center has to be a shape of 3 dimensions!'
-        assert radius > 0 and np.shape(radius) == (), 'Radius has to be positive scalar value!'
-        assert height > 0 and np.shape(height) == (), 'Height has to be positive scalar value!'
-        assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as 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):
-        '''Get the magnetic shape of a sphere.
-        Arguments:
-            dim    - the dimensions of the grid, shape(z, y, x)
-            center - the center of the disc in pixel coordinates, shape(z, y, x)
-            radius - the radius of the disc in pixel coordinates (scalar value)
-            height - the height of the disc in pixel coordinates (scalar value)
-            axis   - the orientation of the disc axis, (String: 'x', 'y', 'z'), default = 'z'
-        Returns:
-            the magnetic shape as a 3D-array with ones and zeros
-
-        '''
-        assert np.shape(dim) == (3,), 'Parameter dim has to be a shape of 3 dimensions!'
-        assert np.shape(center) == (3,), 'Parameter center has to be a shape of 3 dimensions!'
-        assert radius > 0 and np.shape(radius) == (), 'Radius has to be 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'):
-        '''Get the magnetic shape of a single filament in x-, y-, or z-direction.
-        Arguments:
-            dim  - the dimensions of the grid, shape(z, y, x)
-            pos  - the position of the filament in pixel coordinates, shape(coord1, coord2)
-            axis - the orientation of the filament axis, (String: 'x', 'y', 'z'), default = 'y'
-        Returns:
-            the magnetic shape as a 3D-array with ones and zeros
-
-        '''
-        assert np.shape(dim) == (3,), 'Parameter dim has to be a shape of 3 dimensions!'
-        assert np.shape(pos) == (2,), 'Parameter pos has to be a shape of 2 dimensions!'
-        assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as 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):
-        '''Get the magnetic shape of a single magnetized pixel.
-        Arguments:
-            dim   - the dimensions of the grid, shape(z, y, x)
-            pixel - the coordinates of the magnetized pixel, shape(z, y, x)
-        Returns:
-            the magnetic shape as a 3D-array with ones and zeros
-
-        '''
-        assert np.shape(dim) == (3,), 'Parameter dim has to be a shape of 3 dimensions!'
-        assert np.shape(pixel) == (3,), 'Parameter pixel has to be a shape of 3 dimensions!'
-        mag_shape = np.zeros(dim)
-        mag_shape[pixel] = 1
-        return mag_shape
-
-
-def create_mag_dist(mag_shape, phi, theta=pi/2, magnitude=1):
-    '''Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
-    Arguments:
-        mag_shape - the magnetic shapes (numpy arrays, see Shapes.* for examples)
-        phi       - the azimuthal angle, describing the direction of the magnetized object
-        theta     - the polar angle, describing the direction of the magnetized object
-                    (optional, is set to pi/2 if not specified -> z-component is zero)
-        magnitude - the relative magnitudes for the magnetic shape (optional, one if not specified)
-    Returns:
-        the 3D magnetic distribution as a MagData object (see pyramid.magdata for reference)
-
-    '''
-    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_comb(mag_shape_list, phi_list, theta_list=None, magnitude_list=None):
-    '''Create a 3-dimensional magnetic distribution from a list of homogeneous magnetized objects.
-    Arguments:
-        mag_shape_list - a list of magnetic shapes (numpy arrays, see Shapes.* for examples)
-        phi_list       - a list of azimuthal angles, describing the direction of the
-                         magnetized object
-        theta_list     - a list of polar angles, describing the direction of the magnetized object
-                         (optional, is set to pi/2 if not specified -> z-component is zero)
-        magnitude_list - a list of relative magnitudes for the magnetic shapes
-                         (optional, if not specified, every relative magnitude is set to one)
-    Returns:
-        the 3D magnetic distribution as a MagData object (see pyramid.magdata for reference)
-
-    '''
-    # If no relative magnitude is specified, 1 is assumed for every homog. object:
-    if magnitude_list is None:
-        magnitude_list = np.ones(len(phi_list))
-    # If no relative magnitude is specified, 1 is assumed for every homog. object:
-    if theta_list is None:
-        theta_list = np.ones(len(phi_list)) * pi/2
-    # For every shape of a homogeneous object a relative magnitude and angle have to be set:
-    assert np.shape(mag_shape_list)[0] == len(phi_list) == len(theta_list) == len(magnitude_list),\
-        'Lists have not the same length!'
-    dim = np.shape(mag_shape_list[0])  # Has to be the shape of ALL mag_shapes!
-    assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
-    assert np.array([mag_shape_list[i].shape == dim for i in range(len(mag_shape_list))]).all(),\
-        'Magnetic shapes must describe distributions with the same size!'
-    # Start with a zero distribution:
-    x_mag = np.zeros(dim)
-    y_mag = np.zeros(dim)
-    z_mag = np.zeros(dim)
-    # Add every specified homogeneous object:
-    for i in range(np.size(phi_list)):
-        mag_object = create_mag_dist(mag_shape_list[i], phi_list[i], magnitude_list[i])
-        z_mag += mag_object[0]
-        y_mag += mag_object[1]
-        x_mag += mag_object[2]
-    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.
-    Arguments:
-        mag_shape - the magnetic shapes (numpy arrays, see Shapes.* for examples)
-        center    - the vortex center, given in 2D (x,y) or 3D (x,y,z), where z is discarded
-                    (optional, is set to the center of the field of view if not specified, always
-                    between two pixels!)
-        magnitude - the relative magnitudes for the magnetic shape (optional, one if not specified)
-    Returns:
-        the 3D magnetic distribution as a MagData object (see pyramid.magdata for reference)
-
-    '''
-    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!'
-    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

-# -*- coding: utf-8 -*-
-"""Class for creating objects to store magnetizatin data."""
-
-
-import numpy as np
-import tables.netcdf3 as nc
-import matplotlib.pyplot as plt
-from matplotlib.ticker import MaxNLocator
-
-
-class MagData:
-
-    '''An object storing magnetization data.'''
-
-    def __init__(self, res, magnitude):  # TODO: electrostatic component!
-        '''Constructor for a MagData object for storing magnetization data.
-        Arguments:
-            res       - the resolution of the grid (grid spacing) in nm
-            magnitude - the z-, y- and x-component of the magnetization vector for every
-                        3D-gridpoint as a tuple
-        Returns:
-            MagData object
-
-        '''
-        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.res = res
-        self.dim = dim
-        self.magnitude = magnitude
-
-    def get_vector(self, mask):
-        # TODO: DOCSTRING!
-        return np.concatenate([self.magnitude[2][mask],
-                               self.magnitude[1][mask],
-                               self.magnitude[0][mask]])
-
-    def set_vector(self, mask, vector):
-        # TODO: DOCSTRING!
-        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 get_mask(self, threshold=0):
-        # TODO: DOCSTRING!
-        z_mask = abs(self.magnitude[0]) > threshold
-        x_mask = abs(self.magnitude[1]) > threshold
-        y_mask = abs(self.magnitude[2]) > threshold
-        return np.logical_or(np.logical_or(x_mask, y_mask), z_mask)
-
-    def scale_down(self, n=1):
-        # TODO: DOCSTRING!
-        # Starting magnetic distribution:
-        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 DataMag object from LLG-file (classmethod).
-        Arguments:
-            filename - the name of the LLG-file from which to load the data
-        Returns.
-            MagData object
-
-        '''
-        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.
-        Arguments:
-            filename - the name of the LLG-file in which to store the magnetization data
-                       (default: '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_netcdf(cls, filename):
-        '''Construct MagData object from a NetCDF-file (classmethod).
-        Arguments:
-            filename - name of the file from which to load the data
-        Returns:
-            PhaseMap object
-
-        '''
-        f = nc.NetCDFFile(filename, 'r')
-        res = getattr(f, 'res')
-        z_mag = f.variables['z_mag'].getValue()
-        y_mag = f.variables['y_mag'].getValue()
-        x_mag = f.variables['x_mag'].getValue()
-        f.close()
-        return MagData(res, (z_mag, y_mag, x_mag))
-
-    def save_to_netcdf(self, filename='..\output\magdata_output.nc'):
-        '''Save magnetization data in a file with NetCDF-format.
-        Arguments:
-            filename - the name of the file in which to store the phase map data
-                       (default: 'phasemap_output.txt')
-        Returns:
-            None
-
-        '''
-        f = nc.NetCDFFile(filename, 'w')
-        setattr(f, 'res', self.res)
-        f.createDimension('z_dim', self.dim[0])
-        f.createDimension('y_dim', self.dim[1])
-        f.createDimension('x_dim', self.dim[2])
-        z_mag = f.createVariable('z_mag', 'f', ('z_dim', 'y_dim', 'x_dim'))
-        y_mag = f.createVariable('y_mag', 'f', ('z_dim', 'y_dim', 'x_dim'))
-        x_mag = f.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 f
-        f.close()
-
-    def quiver_plot(self, title='Magnetic Distribution)', proj_axis='z', ax_slice=None,
-                    filename=None, axis=None): # TODO!!
-        '''Plot a slice of the magnetization as a quiver plot.
-        Arguments:
-            axis     - the axis from which a slice is plotted ('x', 'y' or 'z'), default = 'z'
-            ax_slice - the slice-index of the specified axis (optional, if not specified, is
-                       set to the center of the specified axis)
-            filename - filename, specifying the location where the image is saved (optional, if
-                       not specified, image is shown instead)
-        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(0, np.shape(mag_slice_u)[1])
-        axis.set_ylim(0, 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):
-        '''3D-Quiver-Plot of the magnetization as vectors.
-        Arguments:
-            None
-        Returns:
-            None
-
-        '''
-        from mayavi import mlab
-
-        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

-from phase_mag_real import *

pyramid/numcore/phase_mag_real.pyx

-import numpy as np
-import math
-cimport cython
-cimport numpy as np
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def phase_mag_real_helper(
-        unsigned int v_dim, unsigned int u_dim,
-        double[:, :] phi_u, double[:, :] phi_v,
-        double[:, :] u_mag, double[:, :] v_mag,
-        double[:, :] phase, float threshold):
-    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 * phi_u[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 * phi_v[q_c + q, p_c + p]

pyramid/phasemap.py

-# -*- coding: utf-8 -*-
-"""Class for creating objects to store phase maps."""
-
-
-import numpy as np
-import matplotlib.pyplot as plt
-import tables.netcdf3 as nc
-
-
-class PhaseMap:
-
-    '''An object storing magnetization data.'''
-
-    def __init__(self, res, phase):
-        '''Constructor for a MagData object for storing magnetization data.
-        Arguments:
-            res       - the resolution of the grid (grid spacing) in nm
-            magnitude - the z-, y- and x-component of the magnetization vector for every
-                        3D-gridpoint as a tuple
-        Returns:
-            MagData object
-
-        '''
-        dim = np.shape(phase)
-        assert len(dim) == 2, 'Phasemap has to be 2-dimensional!'
-        self.res = res
-        self.dim = dim
-        self.phase = phase
-
-    @classmethod
-    def load_from_txt(cls, filename):
-        '''Construct PhaseMap object from a human readable txt-file (classmethod).
-        Arguments:
-            filename - name of the file from which to load the data
-        Returns.
-            PhaseMap object
-
-        '''
-        with open(filename, 'r') as f:
-            f.readline()  # Headerline is not used
-            res = int(f.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 PhaseMap data in a file with txt-format.
-        Arguments:
-            filename - the name of the file in which to store the phase map data
-                       (default: 'phasemap_output.txt')
-        Returns:
-            None
-
-        '''
-        with open(filename, 'w') as f:
-            f.write('{}\n'.format(filename.replace('.txt', '')))
-            f.write('resolution = {} nm\n'.format(self.res))
-            np.savetxt(f, self.phase, fmt='%7.6e', delimiter='\t')
-
-    @classmethod
-    def load_from_netcdf(cls, filename):
-        '''Construct PhaseMap object from a NetCDF-file (classmethod).
-        Arguments:
-            filename - name of the file from which to load the data
-        Returns:
-            PhaseMap object
-
-        '''
-        f = nc.NetCDFFile(filename, 'r')
-        res = getattr(f, 'res')
-        phase = f.variables['phase'].getValue()
-        f.close()
-        return PhaseMap(res, phase)
-
-    def save_to_netcdf(self, filename='..\output\phasemap_output.nc'):
-        '''Save PhaseMap data in a file with NetCDF-format.
-        Arguments:
-            filename - the name of the file in which to store the phase map data
-                       (default: 'phasemap_output.txt')
-        Returns:
-            None
-
-        '''
-        f = nc.NetCDFFile(filename, 'w')
-        setattr(f, 'res', self.res)
-        f.createDimension('v_dim', self.dim[0])
-        f.createDimension('u_dim', self.dim[1])
-        phase = f.createVariable('phase', 'f', ('v_dim', 'u_dim'))
-        phase[:] = self.phase
-        f.close()
-
-    def display(self, title='Phase Map', labels=('x-axis [nm]', 'y-axis [nm]', 'phase [rad]'),
-                cmap='RdBu', limit=None, norm=None, axis=None):
-        '''Display the phasemap as a colormesh.
-        Arguments:
-            title - the title of the plot (default: 'Phase Map')
-            axis  - the axis on which to display the plot (default: None, a new figure is created)
-            cmap  - the colormap which is used for the plot (default: 'gray')
-        Returns:
-            None
-
-        ''' # TODO: Docstring!
-        
-        # TODO: ALWAYS CENTERED around 0
-        if limit is None:
-            limit = np.max(np.abs(self.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(self.phase, cmap=cmap, vmin=-limit, vmax=limit, norm=norm)
-        # Set the axes ticks and labels:
-        axis.set_xlim(0, np.shape(self.phase)[1])
-        axis.set_ylim(0, np.shape(self.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(labels[0], fontsize=15)
-        axis.set_ylabel(labels[1], 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(labels[2], fontsize=15)
-        
-        plt.show()
-
-    def display3d(self, title='Phase Map', labels=('x-axis [nm]', 'y-axis [nm]', 'phase [rad]'),
-                cmap='RdBu', limit=None, norm=None, axis=None):
-        '''Display the phasemap as a colormesh.
-        Arguments:
-            title - the title of the plot (default: 'Phase Map')
-            axis  - the axis on which to display the plot (default: None, a new figure is created)
-            cmap  - the colormap which is used for the plot (default: 'gray')
-        Returns:
-            None
-
-        ''' # TODO: Docstring!
-        
-        from mpl_toolkits.mplot3d import Axes3D
-
-        fig = plt.figure()
-        ax = Axes3D(fig)#.gca(projection='3d')
-
-        u = range(self.dim[1])
-        v = range(self.dim[0])
-        uu, vv = np.meshgrid(u, v)
-        ax.plot_surface(uu, vv, self.phase, rstride=4, cstride=4, alpha=0.7, cmap='RdBu',
-                        linewidth=0, antialiased=False)
-        ax.contourf(uu, vv, self.phase, 15, zdir='z', offset=np.min(self.phase), cmap='RdBu')
-        ax.view_init(45, -135)
-        ax.set_xlabel('x-axis [px]')
-        ax.set_ylabel('y-axis [px]')
-        ax.set_zlabel('phase [mrad]')
-
-        plt.show()

pyramid/phasemapper.py

-# -*- coding: utf-8 -*-
-"""Create and display a phase map from magnetization data."""
-
-
-import numpy as np
-import numcore
-
-from numpy import pi
-
-# Physical constants
-PHI_0 = -2067.83    # magnetic flux in T*nm²
-H_BAR = 6.626E-34  # Planck constant in J*s
-M_E = 9.109E-31    # electron mass in kg
-Q_E = 1.602E-19    # electron charge in C
-C = 2.998E8        # speed of light in m/s
-
-
-def phase_mag_fourier(res, projection, b_0=1, padding=0):
-    '''Calculate phasemap from magnetization data (Fourier space approach).
-    Arguments:
-        res        - the resolution of the grid (grid spacing) in nm
-        projection - projection of a magnetic distribution (created with pyramid.projector)
-        b_0        - magnetic induction corresponding to a magnetization Mo in T (default: 1)
-        padding    - factor for zero padding, the default is 0 (no padding), for a factor of n the
-                     number of pixels is increase by (1+n)**2, padded zeros are cropped at the end
-    Returns:
-        the phasemap as a 2 dimensional array
-
-    '''
-    v_dim, u_dim = np.shape(projection[0])
-    v_mag, u_mag = projection
-    # Create zero padded matrices:
-    u_pad = u_dim/2 * padding
-    v_pad = v_dim/2 * padding
-    u_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
-    v_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
-    u_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = u_mag
-    v_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = v_mag
-    # Fourier transform of the two components:
-    u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_big), axes=0)
-    v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_big), axes=0)
-    # Calculate the Fourier transform of the phase:
-    nyq = 1 / res  # nyquist frequency
-    f_u = np.linspace(0, nyq/2, u_mag_fft.shape[1])
-    f_v = np.linspace(-nyq/2, nyq/2, u_mag_fft.shape[0], endpoint=False)
-    f_uu, f_vv = np.meshgrid(f_u, f_v)
-    coeff = (1j*b_0) / (2*PHI_0)
-    phase_fft = coeff * res * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30)
-    # Transform to real space and revert padding:
-    phase_big = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
-    phase = phase_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
-    return phase
-
-
-def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
-    '''Calculate phasemap from magnetization data (real space approach).
-    Arguments:
-        res        - the resolution of the grid (grid spacing) in nm
-        projection - projection of a magnetic distribution (created with pyramid.projector)
-        method     - String, describing the method to use for the pixel field ('slab' or 'disc')
-        b_0        - magnetic induction corresponding to a magnetization Mo in T (default: 1)
-        jacobi     - matrix in which to save the jacobi matrix (default: None, faster computation)
-    Returns:
-        the phasemap as a 2 dimensional array
-
-    '''
-    # Function for creating the lookup-tables:
-    def phi_lookup(method, n, m, res, b_0):
-        if method == 'slab':
-            def F_h(n, m):
-                a = np.log(res**2 * (n**2 + m**2))
-                b = np.arctan(n / m)
-                return n*a - 2*n + 2*m*b
-            return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
-                                - F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
-        elif method == 'disc':
-            in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
-            return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
-    # Process input parameters:
-    v_dim, u_dim = np.shape(projection[0])
-    v_mag, u_mag = projection
-    coeff = -b_0 * res**2 / (2*PHI_0)
-    # Create lookup-tables for the phase of one pixel:
-    u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
-    v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
-    uu, vv = np.meshgrid(u, v)
-    phi_u = phi_lookup(method, uu, vv, res, b_0)
-    phi_v = phi_lookup(method, vv, uu, res, b_0)
-    # Calculation of the phase:
-    phase = np.zeros((v_dim, u_dim))
-    threshold = 0
-    if jacobi is not None:  # With Jacobian matrix (slower)
-        jacobi[:] = 0  # Jacobi matrix --> zeros
-        ############################### TODO: NUMERICAL CORE  #####################################
-        for j in range(v_dim):
-            for i in range(u_dim):
-                phase_u = phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-                jacobi[:, i+u_dim*j] = phase_u.reshape(-1)
-                if abs(u_mag[j, i]) > threshold:
-                    phase += u_mag[j, i] * phase_u
-                phase_v = phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-                jacobi[:, u_dim*v_dim+i+u_dim*j] = -phase_v.reshape(-1)
-                if abs(v_mag[j, i]) > threshold:
-                    phase -= v_mag[j, i] * phase_v
-        ############################### TODO: NUMERICAL CORE  #####################################
-    else:  # Without Jacobi matrix (faster)
-        numcore.phase_mag_real_helper(v_dim, u_dim, phi_u, phi_v, u_mag, v_mag, phase, threshold)
-    # Return the phase:
-    return phase
-
-
-def phase_mag_real2(res, projection, method, b_0=1, jacobi=None):
-    '''Calculate phasemap from magnetization data (real space approach).
-    Arguments:
-        res        - the resolution of the grid (grid spacing) in nm
-        projection - projection of a magnetic distribution (created with pyramid.projector)
-        method     - String, describing the method to use for the pixel field ('slab' or 'disc')
-        b_0        - magnetic induction corresponding to a magnetization Mo in T (default: 1)
-        jacobi     - matrix in which to save the jacobi matrix (default: None, faster computation)
-    Returns:
-        the phasemap as a 2 dimensional array
-
-    '''
-    # Function for creating the lookup-tables:
-    def phi_lookup(method, n, m, res, b_0):
-        if method == 'slab':
-            def F_h(n, m):
-                a = np.log(res**2 * (n**2 + m**2))
-                b = np.arctan(n / m)
-                return n*a - 2*n + 2*m*b
-            return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
-                                - F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
-        elif method == 'disc':
-            in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
-            return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
-    # Process input parameters:
-    v_dim, u_dim = np.shape(projection[0])
-    v_mag, u_mag = projection
-    coeff = -b_0 * res**2 / (2*PHI_0)
-    # Create lookup-tables for the phase of one pixel:
-    u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
-    v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
-    uu, vv = np.meshgrid(u, v)
-    phi_u = phi_lookup(method, uu, vv, res, b_0)
-    phi_v = phi_lookup(method, vv, uu, res, b_0)
-    # Calculation of the phase:
-    phase = np.zeros((v_dim, u_dim))
-    threshold = 0
-    if jacobi is not None:  # With Jacobian matrix (slower)
-        jacobi[:] = 0  # Jacobi matrix --> zeros
-        ############################### TODO: NUMERICAL CORE  #####################################
-        for j in range(v_dim):
-            for i in range(u_dim):
-                phase_u = phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-                jacobi[:, i+u_dim*j] = phase_u.reshape(-1)
-                if abs(u_mag[j, i]) > threshold:
-                    phase += u_mag[j, i] * phase_u
-                phase_v = phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-                jacobi[:, u_dim*v_dim+i+u_dim*j] = -phase_v.reshape(-1)
-                if abs(v_mag[j, i]) > threshold:
-                    phase -= v_mag[j, i] * phase_v
-        ############################### TODO: NUMERICAL CORE  #####################################
-    else:  # Without Jacobi matrix (faster)
-#        phasecopy = phase.copy()
-#        start_time = time.time()
-#        numcore.phase_mag_real_helper_1(v_dim, u_dim, phi_u, phi_v,
-#                                        u_mag, v_mag, phasecopy, threshold)
-#        print 'with numcore   : ', time.time() - start_time
-#        start_time = time.time()
-        for j in range(v_dim):
-            for i in range(u_dim):
-                if abs(u_mag[j, i]) > threshold:
-                    phase += u_mag[j, i] * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-                if abs(v_mag[j, i]) > threshold:
-                    phase -= v_mag[j, i] * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-#        print 'without numcore: ', time.time() - start_time
-#        print 'num. difference: ', ((phase - phasecopy) ** 2).sum()
-    # Return the phase:
-    return phase
-
-
-def phase_mag_real_alt(res, projection, method, b_0=1, jacobi=None):  # TODO: Modify
-    '''Calculate phasemap from magnetization data (real space approach).
-    Arguments:
-        res        - the resolution of the grid (grid spacing) in nm
-        projection - projection of a magnetic distribution (created with pyramid.projector)
-        method     - String, describing the method to use for the pixel field ('slab' or 'disc')
-        b_0        - magnetic induction corresponding to a magnetization Mo in T (default: 1)
-        jacobi     - matrix in which to save the jacobi matrix (default: None, faster computation)
-    Returns:
-        the phasemap as a 2 dimensional array
-
-    '''
-    # Function for creating the lookup-tables:
-    def phi_lookup(method, n, m, res, b_0):
-        if method == 'slab':
-            def F_h(n, m):
-                a = np.log(res**2 * (n**2 + m**2))
-                b = np.arctan(n / m)
-                return n*a - 2*n + 2*m*b
-            return coeff * 0.5 * (F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5)
-                                - F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5))
-        elif method == 'disc':
-            in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
-            return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
-
-    # Function for the phase contribution of one pixel:
-    def phi_mag(i, j):
-        return (np.cos(beta[j, i]) * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-              - np.sin(beta[j, i]) * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
-
-    # Function for the derivative of the phase contribution of one pixel:
-    def phi_mag_deriv(i, j):
-        return -(np.sin(beta[j, i]) * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-               + np.cos(beta[j, i]) * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
-
-    # Process input parameters:
-    v_dim, u_dim = np.shape(projection[0])
-    v_mag, u_mag = projection
-    beta = np.arctan2(v_mag, u_mag)
-    mag = np.hypot(u_mag, v_mag)
-    coeff = -b_0 * res**2 / (2*PHI_0)
-    # Create lookup-tables for the phase of one pixel:
-    u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
-    v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
-    uu, vv = np.meshgrid(u, v)
-    phi_u = phi_lookup(method, uu, vv, res, b_0)
-    phi_v = phi_lookup(method, vv, uu, res, b_0)
-    # Calculation of the phase:
-    phase = np.zeros((v_dim, u_dim))
-    threshold = 0
-    if jacobi is not None:  # With Jacobian matrix (slower)
-        jacobi[:] = 0  # Jacobi matrix --> zeros
-        ############################### TODO: NUMERICAL CORE  #####################################
-        for j in range(v_dim):
-            for i in range(u_dim):
-                phase_cache = phi_mag(i, j)
-                jacobi[:, i+u_dim*j] = phase_cache.reshape(-1)
-                if mag[j, i] > threshold:
-                    phase += mag[j, i]*phase_cache
-                    jacobi[:, u_dim*v_dim+i+u_dim*j] = (mag[j, i]*phi_mag_deriv(i, j)).reshape(-1)
-        ############################### TODO: NUMERICAL CORE  #####################################
-    else:  # Without Jacobi matrix (faster)
-        for j in range(v_dim):
-            for i in range(u_dim):
-                if abs(mag[j, i]) > threshold:
-                    phase += mag[j, i] * phi_mag(i, j)
-    # Return the phase:
-    return phase
-
-
-def phase_elec(mag_data, v_0=1, v_acc=30000):
-    # TODO: Docstring
-
-    res  = mag_data.res
-    z_dim, y_dim, x_dim = mag_data.dim
-    z_mag, y_mag, x_mag = mag_data.magnitude
-
-    phase = np.logical_or(x_mag, y_mag, z_mag)
-
-    lam = H_BAR / np.sqrt(2 * M_E * v_acc * (1 + Q_E*v_acc / (2*M_E*C**2)))
-
-    Ce = (2*pi*Q_E/lam * (Q_E*v_acc +   M_E*C**2) /
-            (Q_E*v_acc * (Q_E*v_acc + 2*M_E*C**2)))
-
-    phase *= res * v_0 * Ce
-
-    return phase

pyramid/projector.py

-# -*- coding: utf-8 -*-
-"""Planar projection of the magnetization distribution of a MagData object."""
-
-
-from pyramid.magdata import MagData
-
-
-def simple_axis_projection(mag_data, axis='z'):
-    '''Project a magnetization distribution along one of the main axes of the 3D-grid.
-    Arguments:
-        mag_data - a MagData object storing the magnetization distribution
-        axis     - the projection direction (String: 'x', 'y', 'z'), default = 'z'
-    Returns:
-        the in-plane projection of the magnetization as a tuple: (x_mag, y_mag)
-        ()
-
-    '''
-    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
-    elif axis == 'y':
-        projection = (mag_data.magnitude[0].sum(1),  # z_mag -> v_mag
-                      mag_data.magnitude[2].sum(1))  # x_mag -> u_mag
-    elif axis == 'x':
-        projection = (mag_data.magnitude[0].sum(2),  # z_mag -> v_mag
-                      mag_data.magnitude[1].sum(2))  # y_mag -> u_mag
-    return projection
-
-# TODO: proper projection algorithm with two angles and such!
-# CAUTION: the res for the projection does not have to be the res of the 3D-magnetization!
-# Just for a first approach

pyramid/reconstructor.py

-# -*- coding: utf-8 -*-
-"""Reconstruct magnetic distributions with given phasemaps"""
-
-
-import numpy as np
-import pyramid.projector as pj
-import pyramid.phasemapper as pm
-from pyramid.magdata import MagData
-from scipy.optimize import leastsq
-
-
-def reconstruct_simple_leastsq(phase_map, mask, b_0):
-    '''Reconstruct a magnetic distribution where the positions of the magnetized voxels are known
-    from a single phase_map using the least square method (only works for slice thickness = 1)
-    Arguments:
-        phase_map - a PhaseMap object, from which to reconstruct the magnetic distribution
-        mask      - a boolean matrix representing the positions of the magnetized voxels (3D)
-        b_0       - magnetic induction corresponding to a magnetization Mo in T (default: 1)
-    Returns:
-        the reconstructed magnetic distribution (as a MagData object)
-
-    '''
-    # 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_discs.py

-#! python
-# -*- coding: utf-8 -*-
-"""Create the Pyramid-Logo."""
-
-
-import pdb
-import traceback
-import sys
-import numpy as np
-from numpy import pi
-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
-
-
-PHI_0 = -2067.83  # magnetic flux in T*nm²
-
-
-def compare_vortices():
-    ''' DOCSTRING'''  # TODO: Docstring!
-    # Input parameters:
-    dim_list = [(16, 256, 256), (8, 128, 128), (4, 64, 64), (2, 32, 32), (1, 16, 16)]
-    res_list = [1., 2., 4., 8., 16.]  # in nm
-    density = 20
-    phi = pi/2
-
-    x = []
-    y = []
-
-    # Analytic solution:
-    L = dim_list[0][1]  # in px/nm
-    Lz = 0.5 * dim_list[0][0]  # in px/nm
-    R = 0.25 * L  # in px/nm
-    x0 = L / 2  # in px/nm
-
-    def F(x):
-        coeff = - pi * Lz / (2*PHI_0)
-        result = coeff * (- (x - x0) * np.sin(phi))
-        result *= np.where(np.abs(x - x0) <= R, 1, (R / (x - x0)) ** 2)
-        return result
-    x_an = np.linspace(0, L, 1000)
-    y_an = F(x_an)
-
-    # Starting magnetic distribution:
-    dim_start = (2*dim_list[0][0], 2*dim_list[0][1], 2*dim_list[0][2])
-    res_start = res_list[0] / 2
-    center = (dim_start[0]/2 - 0.5, dim_start[1]/2 - 0.5, dim_start[2]/2 - 0.5)
-    radius = 0.25 * dim_start[1]
-    height = 0.5 * dim_start[0]
-    mag_shape = mc.Shapes.disc(dim_start, center, radius, height)
-    mag_data = MagData(res_start, mc.create_mag_dist(mag_shape, phi))
-
-    for i, (dim, res) in enumerate(zip(dim_list, res_list)):
-        # Create coarser grid for the magnetization:
-        print 'dim = ', dim, 'res = ', res
-        z_mag = mag_data.magnitude[0].reshape(dim[0], 2, dim[1], 2, dim[2], 2)
-        y_mag = mag_data.magnitude[1].reshape(dim[0], 2, dim[1], 2, dim[2], 2)
-        x_mag = mag_data.magnitude[2].reshape(dim[0], 2, dim[1], 2, dim[2], 2)
-        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))
-        mag_data = MagData(res, magnitude)
-        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, res = {}'.format(res))
-        x.append(np.linspace(0, dim[1]*res, dim[1]))
-        y.append(phase_map.phase[dim[1]/2, :])
-
-    # Plot:
-    fig = plt.figure()
-    axis = fig.add_subplot(1, 1, 1)
-    axis.plot(x[0], y[0], 'r',
-              x[1], y[1], 'm',
-              x[2], y[2], 'y',
-              x[3], y[3], 'g',
-              x[4], y[4], 'c',
-              x_an, y_an, 'k')
-    axis.set_xlabel('x [nm]')
-    axis.set_ylabel('phase')
-
-
-if __name__ == "__main__":
-    try:
-        compare_vortices()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)

scripts/compare_method_errors_fourier_padding.py

-#! 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 matplotlib.pyplot as plt
-
-import pyramid.magcreator as mc
-import pyramid.projector as pj
-import pyramid.phasemapper as pm
-import pyramid.analytic as an
-from pyramid.magdata import MagData
-from pyramid.phasemap import PhaseMap
-import shelve
-
-
-def phase_from_mag():
-    '''Calculate and display the phase map from a given magnetization.
-    Arguments:
-        None
-    Returns:
-        None
-
-    '''
-    # Create / Open databank:
-    data_shelve = shelve.open('../output/method_errors_shelve')
-
-    '''FOURIER PADDING->RMS|DURATION'''
-    # Parameters:
-    b_0 = 1  # in T
-    res = 10.0  # in nm
-    dim = (1, 128, 128)
-    phi = -pi/2
-    padding_list = [0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
-                    16, 17, 18, 19, 20, 21, 22]
-    geometry = 'disc'
-    # Create magnetic shape:
-    if geometry == 'slab':
-        center = (0, dim[1]/2-0.5, dim[2]/2-0.5)  # in px (z, y, x) index starts with 0!
-        width = (1, 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 = (0, dim[1]/2-0.5, dim[2]/2-0.5)  # in px (z, y, x) index starts with 0!
-        radius = dim[1] / 4  # in px
-        height = 1  # in px
-        mag_shape = mc.Shapes.disc(dim, center, radius, height)
-        phase_ana = an.phase_mag_disc(dim, res, phi, center, radius, b_0)
-    # Project the magnetization data:
-    mag_data = MagData(res, mc.create_mag_dist(mag_shape, phi))
-    projection = pj.simple_axis_projection(mag_data)
-    # Create data:
-    data = np.zeros((3, len(padding_list)))
-    data[0, :] = padding_list
-    for (i, padding) in enumerate(padding_list):
-        print 'padding =', padding_list[i]
-        # Key:
-        key = ', '.join(['Padding->RMS|duration', 'Fourier', 'padding={}'.format(padding_list[i]),
-                        'B0={}'.format(b_0), 'res={}'.format(res), 'dim={}'.format(dim),
-                        'phi={}'.format(phi), 'geometry={}'.format(geometry)])
-        if key in data_shelve:
-            data[:, i] = data_shelve[key]
-        else:
-            start_time = time.time()
-            phase_num = pm.phase_mag_fourier(res, projection, b_0, padding_list[i])
-            data[2, i] = time.time() - start_time
-            phase_diff = phase_ana - phase_num
-            PhaseMap(res, phase_diff).display()
-            data[1, i] = np.std(phase_diff)
-            data_shelve[key] = data[:, i]
-    # Plot duration against padding:
-    fig = plt.figure()
-    axis = fig.add_subplot(1, 1, 1)
-    axis.set_yscale('log')
-    axis.plot(data[0], data[1])
-    axis.set_title('Fourier Space Approach: Variation of the Padding')
-    axis.set_xlabel('padding')
-    axis.set_ylabel('RMS')
-    # Plot RMS against padding:
-    fig = plt.figure()
-    axis = fig.add_subplot(1, 1, 1)
-    axis.plot(data[0], data[2])
-    axis.set_title('Fourier Space Approach: Variation of the Padding')
-    axis.set_xlabel('padding')
-    axis.set_ylabel('duration [s]')
-    # Close shelve:
-    data_shelve.close()
-
-
-if __name__ == "__main__":
-    try:
-        phase_from_mag()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)

scripts/compare_methods.py

-#! 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    # in T
-    res = 10.0  # in nm
-    phi = pi/4
-    padding = 12
-    density = 1
-    dim = (16, 128, 128)  # in px (z, y, x)
-    # Create magnetic shape:
-    geometry = 'disc'
-    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 with 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 with 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(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, b_0, padding))
-    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
-    # 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)')
-    # 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)
-    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)
-    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))
-
-
-if __name__ == "__main__":
-    try:
-        compare_methods()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)

scripts/compare_pixel_fields.py

-#! python
-# -*- coding: utf-8 -*-
-"""Compare the phase map of one pixel for different real space approaches."""
-
-
-import pdb
-import traceback
-import sys
-from numpy import pi
-
-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
-
-    '''
-    # Input parameters:
-    res = 10.0  # in nm
-    phi = pi/2  # in rad
-    dim = (1, 5, 5)
-    pixel = (0,  2,  2)
-    # Create magnetic data, project it, get the phase map and display the holography image:
-    mag_data = MagData(res, mc.create_mag_dist(mc.Shapes.pixel(dim, pixel), phi))
-    mag_data.save_to_llg('../output/mag_dist_single_pixel.txt')
-    projection = pj.simple_axis_projection(mag_data)
-    phase_map_slab = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
-    phase_map_slab.display('Phase of one Pixel (Slab)')
-    phase_map_disc = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc'))
-    phase_map_disc.display('Phase of one Pixel (Disc)')
-    phase_map_diff = PhaseMap(res, phase_map_disc.phase - phase_map_slab.phase)
-    phase_map_diff.display('Phase difference of one Pixel (Disc - Slab)')
-
-if __name__ == "__main__":
-    try:
-        compare_pixel_fields()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)

scripts/compare_vortices.py

-#! python
-# -*- coding: utf-8 -*-
-"""Create the Pyramid-Logo."""
-
-
-import pdb
-import traceback
-import sys
-import numpy as np
-from numpy import pi
-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
-
-
-PHI_0 = -2067.83  # magnetic flux in T*nm²
-
-
-def compare_vortices():
-    '''Calculate and display the Pyramid-Logo.
-    Arguments:
-        None
-    Returns:
-        None
-
-    '''
-    # Input parameters:
-    dim_list = [(16, 256, 256), (8, 128, 128), (4, 64, 64), (2, 32, 32), (1, 16, 16)]
-    res_list = [1., 2., 4., 8., 16.]  # in nm
-    density = 20
-
-    x = []
-    y = []
-
-    # Starting magnetic distribution:
-    dim_start = (2*dim_list[0][0], 2*dim_list[0][1], 2*dim_list[0][2])
-    res_start = res_list[0] / 2
-    center = (dim_start[0]/2 - 0.5, dim_start[1]/2 - 0.5, dim_start[2]/2 - 0.5)
-    radius = 0.25 * dim_start[1]
-    height = 0.5 * dim_start[0]
-    mag_shape = mc.Shapes.disc(dim_start, center, radius, height)
-    mag_data = MagData(res_start, mc.create_mag_dist_vortex(mag_shape))
-
-    # Analytic solution:
-    L = dim_list[0][1]  # in px/nm
-    Lz = 0.5 * dim_list[0][0]  # in px/nm
-    R = 0.25 * L  # in px/nm
-    x0 = L / 2  # in px/nm
-
-    def F(x):
-        coeff = pi*Lz/PHI_0
-        result = coeff * np.where(np.abs(x - x0) <= R, (np.abs(x-x0)-R), 0)
-        return result
-    x_an = np.linspace(0, L, 1001)
-    y_an = F(x_an)
-
-    for i, (dim, res) in enumerate(zip(dim_list, res_list)):
-        # Create coarser grid for the magnetization:
-        print 'dim = ', dim, 'res = ', res
-        z_mag = mag_data.magnitude[0].reshape(dim[0], 2, dim[1], 2, dim[2], 2)
-        y_mag = mag_data.magnitude[1].reshape(dim[0], 2, dim[1], 2, dim[2], 2)
-        x_mag = mag_data.magnitude[2].reshape(dim[0], 2, dim[1], 2, dim[2], 2)
-        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))
-        mag_data = MagData(res, magnitude)
-        #mag_data.quiver_plot()
-        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, res = {}'.format(res))
-        x.append(np.linspace(0, dim[1]*res, dim[1]))
-        y.append(phase_map.phase[dim[1]/2, :])
-
-    # Plot:
-    fig = plt.figure()
-    axis = fig.add_subplot(1, 1, 1)
-    axis.plot(x[0], y[0], 'r',
-              x[1], y[1], 'm',
-              x[2], y[2], 'y',
-              x[3], y[3], 'g',
-              x[4], y[4], 'c',
-              x_an, y_an, 'k')
-    axis.set_xlabel('x [nm]')
-    axis.set_ylabel('phase')
-
-
-if __name__ == "__main__":
-    try:
-        compare_vortices()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)

scripts/create_alternating_filaments.py

-#! python
-# -*- coding: utf-8 -*-
-"""Create magnetic distribution of alternating filaments"""
-
-
-import pdb
-import traceback
-import sys
-import numpy as np
-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
-
-    '''
-    # Input parameters:
-    filename = '../output/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 lists for magnetic objects:
-    count = int((dim[1]-1) / spacing) + 1
-    mag_shape_list = np.zeros((count,) + dim)
-    phi_list = np.zeros(count)
-    for i in range(count):
-        pos = i * spacing
-        mag_shape_list[i] = mc.Shapes.filament(dim, (0, pos))
-        phi_list[i] = (1-2*(int(pos/spacing) % 2)) * phi
-    # Create magnetic distribution
-    magnitude = mc.create_mag_dist_comb(mag_shape_list, phi_list)
-    mag_data = MagData(res, magnitude)
-    mag_data.quiver_plot()
-    mag_data.save_to_llg(filename)
-
-
-if __name__ == "__main__":
-    try:
-        create_alternating_filaments()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)

scripts/create_logo.py

-#! 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(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_multiple_samples.py

-# -*- coding: utf-8 -*-
-"""Create random magnetic distributions."""
-
-import random as rnd
-import pdb
-import traceback
-import sys
-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_multiple_samples():
-    '''Calculate, display and save a random magnetic distribution to file.
-    Arguments:
-        None
-    Returns:
-        None
-
-    '''
-    filename = '../output/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 lists for magnetic objects:
-    mag_shape_list = [mag_shape_slab, mag_shape_disc, mag_shape_sphere]
-    phi_list = [pi/4, pi/2, pi]
-    # Create magnetic distribution:
-    magnitude = mc.create_mag_dist_comb(mag_shape_list, phi_list)
-    mag_data = MagData(res, magnitude)
-    mag_data.quiver_plot('z', dim[0]/2)
-    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_random_pixels.py

-# -*- coding: utf-8 -*-
-"""Create random magnetic distributions."""
-
-import random as rnd
-import pdb
-import traceback
-import sys
-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
-
-    '''
-    # Input parameters:
-    count = 10
-    dim = (1, 128, 128)
-    res = 10  # in nm
-    rnd.seed(12)
-    # Create lists for magnetic objects:
-    mag_shape_list = np.zeros((count,) + dim)
-    phi_list = np.zeros(count)
-    magnitude_list = np.zeros(count)
-    for i in range(count):
-        pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
-        mag_shape_list[i, ...] = mc.Shapes.pixel(dim, pixel)
-        phi_list[i] = 2 * pi * rnd.random()
-        magnitude_list[i] = 1  # TODO: rnd.random()
-    # Create magnetic distribution:
-    magnitude = mc.create_mag_dist_comb(mag_shape_list, phi_list, magnitude_list)
-    mag_data = MagData(res, magnitude)
-    mag_data.quiver_plot()
-    mag_data.save_to_llg('../output/mag_dist_random_pixels.txt')
-    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_random_slabs.py

-# -*- coding: utf-8 -*-
-"""Create random magnetic distributions."""
-
-import random as rnd
-import pdb
-import traceback
-import sys
-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
-
-    '''
-    # Input parameters:
-    count = 10
-    dim = (1, 128, 128)
-    res = 10  # in nm
-    rnd.seed(42)
-    w_max = 10
-    # Create lists for magnetic objects:
-    mag_shape_list = np.zeros((count,) + dim)
-    phi_list = np.zeros(count)
-    magnitude_list = np.zeros(count)
-    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_list[i, ...] = mc.Shapes.slab(dim, center, width)
-        phi_list[i] = 2 * pi * rnd.random()
-        magnitude_list[i] = 1  # TODO: rnd.random()
-    # Create magnetic distribution:
-    magnitude = mc.create_mag_dist_comb(mag_shape_list, phi_list, magnitude_list)
-    mag_data = MagData(res, magnitude)
-    mag_data.quiver_plot()
-    mag_data.save_to_llg('../output/mag_dist_random_slabs.txt')
-    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)