Completed new program structure!

Scripts are nearly adapted to the new structure and easier to read.
(Except analytic module and scripts get_jacobi and inverse)

pyramid/analytic.py

@@ -27,7 +27,7 @@ def plot_phase(phase, res, name):
 
 def phasemap_slab(dim, res, beta, center, width, b_0):
     '''INPUT VARIABLES'''
-    y_dim, x_dim = dim
+    z_dim, y_dim, x_dim = dim
     # y0, x0 have to be in the center of a pixel, hence: cellindex + 0.5
     y0 = res * (center[0] + 0.5)
     x0 = res * (center[1] + 0.5)

pyramid/holoimage.py

@@ -5,121 +5,122 @@
 import numpy as np
 import matplotlib as mpl
 import matplotlib.pyplot as plt
+from pyramid.phasemap import PhaseMap
 from numpy import pi
 from PIL import Image
-
-
-class HoloImage:
-    
-    '''An object storing magnetization data.'''
     
 
-    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)],
+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)]}
+         '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)
+HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
     
     
-    def __init__(self, phase, res, density):
-        '''Returns a holography image with color-encoded gradient direction.
-        Arguments:
-            phase   - the phasemap that should be displayed
-            res     - the resolution of the phasemap
-            density - the factor for determining the number of contour lines
-        Returns:
-            Image
-            
-        '''
-        img_holo = (1 + np.cos(density * phase * pi/2)) /2
-        
-        phase_grad_y, phase_grad_x = np.gradient(phase, res, res)
-        
-        phase_angle = (1 - np.arctan2(phase_grad_y, phase_grad_x)/pi) / 2
-        
-        # TODO: Delete
-    #    import pdb; pdb.set_trace()      
+def holo_image(phase_map, density):
+    '''Returns a holography image with color-encoded gradient direction.
+    Arguments:
+        phase_map - a PhaseMap object storing the phase informations
+        density   - the factor for determining the number of contour lines
+    Returns:
+        holography image
         
-        phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)    
-        phase_magnitude /= phase_magnitude.max()
-        phase_magnitude = np.sin(phase_magnitude * pi / 2)
+    '''
+    assert isinstance(phase_map, PhaseMap), 'phase_map has to be a PhaseMap object!'
+    # Calculate the holography image intensity:
+    img_holo = (1 + np.cos(density * phase_map.phase * pi/2)) /2
+    # Calculate the phase gradients, expressed by magnitude and angle:
+    phase_grad_y, phase_grad_x = np.gradient(phase_map.phase, phase_map.res, phase_map.res)
+    phase_angle = (1 - np.arctan2(phase_grad_y, phase_grad_x)/pi) / 2
+    phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)    
+    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
         
-        cmap = HoloImage.HOLO_CMAP
-        rgba_img = cmap(phase_angle)
-        rgb_img = np.delete(rgba_img, 3, 2)
-        red, green, blue = rgb_img[:, :, 0], rgb_img[:, :, 1], rgb_img[:, :, 2]    
-        red   *= 255.999 * img_holo * phase_magnitude
-        green *= 255.999 * img_holo * phase_magnitude
-        blue  *= 255.999 * img_holo * phase_magnitude
-        rgb = np.dstack((red, green, blue)).astype(np.uint8)
-        # TODO: Which one?
-        rgb = (255.999 * img_holo.T * phase_magnitude.T 
-               * rgba_img[:, :, :3].T).T.astype(np.uint8)
-        img = Image.fromarray(rgb)
-        
-        self.image = img
+    
+def make_color_wheel():
+    '''Display a color wheel for the gradient direction.
+    Arguments:
+        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()
+    ax = fig.add_subplot(111, aspect='equal')
+    ax.imshow(color_wheel)
+    ax.set_title('Color Wheel')
+    ax.set_xlabel('x-axis')
+    ax.set_ylabel('y-axis')
     
-    @classmethod
-    def make_color_wheel(self):
-        '''Display a color wheel for the gradient direction.
-        Arguments:
-            None
-        Returns:
-            None
-            
-        '''
-        x = np.linspace(-256, 256, num=512)
-        y = np.linspace(-256, 256, num=512)
-        xx, yy = np.meshgrid(x, y)
-        color_wheel_angle = (1 - np.arctan2(xx, -yy)/pi) / 2
+    
+def display(holo_image, title='Holography Image', axis=None):
+    '''Display the color coded holography image resulting from a given phase map.
+    Arguments:
+        holo_image - holography image created with the holo_image function of this module
+        title      - the title of the plot (default: 'Holography Image')
+        axis       - the axis on which to display the plot (default: None, a new figure is created)
+    Returns:
+        None
             
-        r = np.sqrt(xx ** 2 + yy ** 2)
-        color_wheel_magnitude = (1 - np.cos(r * pi/360)) / 2
-        color_wheel_magnitude *= 0 * (r > 256) + 1 * (r <= 256)
-        
-        cmap = HoloImage.HOLO_CMAP
-        rgba_img = cmap(color_wheel_angle)
-        rgb_img = np.delete(rgba_img, 3, 2)
-        red, green, blue = rgb_img[:, :, 0], rgb_img[:, :, 1], rgb_img[:, :, 2]    
-        red   *= 255.999 * color_wheel_magnitude
-        green *= 255.999 * color_wheel_magnitude
-        blue  *= 255.999 * color_wheel_magnitude
-        #rgb = np.dstack((red, green, blue)).astype(np.uint8)
-        # TODO Evtl. einfacher:
-        rgb = (255.999 * color_wheel_magnitude.T * rgba_img[:, :, :3].T).T.astype(np.uint8)
-        color_wheel = Image.fromarray(rgb)
-        
+    '''
+    # If no axis is specified, a new figure is created:
+    if axis == None:    
         fig = plt.figure()
-        ax = fig.add_subplot(111, aspect='equal')
-        ax.imshow(color_wheel)
-        ax.set_title('Color Wheel')
-        ax.set_xlabel('x-axis')
-        ax.set_ylabel('y-axis')
-        
+        axis = fig.add_subplot(1,1,1, aspect='equal')
+    # Plot the image and set axes:
+    axis.imshow(holo_image)
+    axis.set_title(title)
+    axis.set_xlabel('x-axis')
+    axis.set_ylabel('y-axis')
     
-    def display_holo(holo, title, axis=None):
-        # TODO: Docstring
-        if axis == None:    
-            fig = plt.figure()
-            axis = fig.add_subplot(1,1,1, aspect='equal')
+    
+def display_combined(phase_map, density, title='Combined Plot'):
+    '''Display a given phase map and the resulting color coded holography image in one plot.
+    Arguments:
+        phase_map - the PhaseMap object from which the holography image is calculated
+        density   - the factor for determining the number of contour lines
+        title     - the title of the combined plot (default: 'Combined Plot')
+    Returns:
+        None
             
-        axis.imshow(holo)
-        
-        axis.set_title(title)
-        axis.set_xlabel('x-axis')
-        axis.set_ylabel('y-axis')
+    '''
+    # Create combined plot and set title:
+    fig = plt.figure(figsize=(14, 7))
+    fig.suptitle(title, fontsize=20)
+    # Plot holography image:
\ No newline at end of file
+    holo_axis = fig.add_subplot(1,2,1, aspect='equal')
+    display(holo_image(phase_map, density), axis=holo_axis)
+    # Plot phase map:
+    phase_axis = fig.add_subplot(1,2,2, aspect='equal')
+    phase_map.display(axis=phase_axis)
\ No newline at end of file

pyramid/magcreator.py

@@ -3,130 +3,184 @@
 
 
 import numpy as np
-from magdata import MagData
+from pyramid.magdata import MagData
 
 
-def shape_slab(dim, center, width):
-    '''Get the magnetic shape of a slab.
-    Arguments:
-        dim    - the dimensions of the grid, shape(y, x)
-        center - center of the slab in pixel coordinates, shape(y, x)
-        width  - width of the slab in pixel coordinates, shape(y, x)
-    Returns:
-        the magnetic shape as a 2D-boolean-array.
-        
-    '''
-    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
-
-
-def shape_disc(dim, center, radius, height):
-    '''Get the magnetic shape of a disc.
-    Arguments:
-        dim    - the dimensions of the grid, shape(y, x)
-        center - center of the disc in pixel coordinates, shape(y, x)
-        radius - radius of the disc in pixel coordinates, shape(y, x)
-    Returns:
-        the magnetic shape as a 2D-boolean-array.
-        
-    '''# TODO: up till now only in z-direction
-    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])])    
-    return mag_shape
-
+class Shapes:
     
-def shape_filament(dim, pos, x_or_y):
-    '''Get the magnetic shape of a single filament in x or y direction.
-    Arguments:
-        dim    - the dimensions of the grid, shape(y, x)
-        pos    - position of the filament (pixel coordinate)
-        x_or_y - string that determines the orientation of the filament
-                 ('y' or 'x')
-    Returns:
-        the magnetic shape as a 2D-boolean-array.
-        
-    '''# TODO: up till now no z-direction
-    mag_shape = np.zeros(dim)
-    if x_or_y == 'y':
-        mag_shape[:, :, pos] = 1
-    else:
-        mag_shape[:, pos, :] = 1
-    return mag_shape
-
-
-def shape_alternating_filaments(dim, spacing, x_or_y):
-    '''Get the magnetic shape of alternating filaments in x or y direction.
-    Arguments:
-        dim     - the dimensions of the grid, shape(y, x)
-        spacing - the distance between two filaments (pixel coordinate)
-        x_or_y  - string that determines the orientation of the filament
-                  ('y' or 'x')
-    Returns:
-        the magnetic shape as a 2D-boolean-array.
-        
-    '''
-    mag_shape = np.zeros(dim)
-    if x_or_y == 'y':
-        for i in range(0, dim[1], spacing):
-            mag_shape[:, :, i] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1
-    else:
-        for i in range(0, dim[0], spacing):
-            mag_shape[:, i, :] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1            
-    return mag_shape
-
+    '''Class containing functions for generating shapes.'''
+       
+    @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
     
-def shape_single_pixel(dim, pixel):
-    '''Get the magnetic shape of a single magnetized pixel.
+    @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 single_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, beta, magnitude=1):
+    '''Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
     Arguments:
-        dim   - the dimensions of the grid, shape(y, x)
-        pixel - the coordinates of the magnetized pixel, shape(y, x)
+        mag_shape - the magnetic shapes (numpy arrays, see Shapes.* for examples)
+        beta      - the angle, describing the direction of the magnetized object
+        magnitude - the relative magnitudes for the magnetic shape (optional, one if not specified)
     Returns:
-        the magnetic shape as a 2D-boolean-array.
+        the 3D magnetic distribution as a MagData object (see pyramid.magdata for reference)
         
     '''
-    mag_shape = np.zeros(dim)
-    mag_shape[pixel] = 1
-    return mag_shape
+    dim = np.shape(mag_shape)
+    assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
+    z_mag = np.array(np.zeros(dim))  # TODO: Implement another angle!
+    y_mag = np.array(np.ones(dim)) * np.sin(beta) * mag_shape * magnitude
+    x_mag = np.array(np.ones(dim)) * np.cos(beta) * mag_shape * magnitude
+    return z_mag, y_mag, x_mag
 
 
-def hom_mag(dim, res, mag_shape, beta, magnitude=1):
-    '''Create homog. magnetization data, saved in a file with LLG-format.
+def create_mag_dist_comb(mag_shape_list, beta_list, magnitude_list=None):
+    '''Create a 3-dimensional magnetic distribution from a list of homogeneous magnetized objects.
     Arguments:
-        dim       - the dimensions of the grid, shape(y, x)
-        res       - the resolution of the grid 
-                    (real space distance between two points)
-        beta      - the angle of the magnetization
-        filename  - the name of the file in which to save the data
-        mag_shape - an array of shape dim, representing the shape of the magnetic object,
-                    a few are supplied in this module
+        mag_shape_list - a list of magnetic shapes (numpy arrays, see Shapes.* for examples)
+        beta_list      - a list of angles, describing the direction of the magnetized objects
+        magnitude_list - a list of relative magnitudes for the magnetic shapes
+                         (optional, if not specified, every relative magnitude is set to one)
     Returns:
-        the the magnetic distribution as a 2D-boolean-array.
+        the 3D magnetic distribution as a MagData object (see pyramid.magdata for reference)
         
-    ''' # TODO: renew Docstring
-    z_mag = np.array(np.zeros(dim))        
-    y_mag = np.array(np.ones(dim)) * np.sin(beta) * mag_shape * magnitude
-    x_mag = np.array(np.ones(dim)) * np.cos(beta) * mag_shape * magnitude
-    return z_mag, y_mag, x_mag
-    
-                                      
-def create_mag_dist(dim, res, mag_shape_list, beta_list, magnitude_list=None):
-    # TODO: Docstring: can take just one or a list of objects!  OR JUST A LIST!!!  
+    '''
     # If no relative magnitude is specified, 1 is assumed for every homog. object:
     if magnitude_list is None:
-        magnitude_list = np.ones(np.size(beta_list))    
+        magnitude_list = np.ones(len(beta_list))    
     # For every shape of a homogeneous object a relative magnitude and angle have to be set:
-    assert (np.shape(mag_shape_list)[0],) == np.shape(beta_list) == np.shape(magnitude_list), \
-        'Lists have not the same length!'
+    assert np.shape(mag_shape_list)[0] == len(beta_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)    
+    z_mag = np.zeros(dim)
+    # Add every specified homogeneous object:
     for i in range(np.size(beta_list)):
-        pixel_mag = hom_mag(dim, res, mag_shape_list[i], beta_list[i], magnitude_list[i])    
-        z_mag += pixel_mag[0]
-        y_mag += pixel_mag[1]
-        x_mag += pixel_mag[2]
-    return MagData(res, z_mag, y_mag, x_mag)
+        mag_object = create_mag_dist(mag_shape_list[i], beta_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
\ No newline at end of file

pyramid/magdata.py

 # -*- coding: utf-8 -*-
-"""Load magnetization data from LLG files."""
+"""Class for creating objects to store magnetizatin data."""
 
 
 import numpy as np
 import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
 from matplotlib.patches import FancyArrowPatch
 from mpl_toolkits.mplot3d import proj3d
 
@@ -13,64 +12,66 @@ class MagData:
     
     '''An object storing magnetization data.'''
     
-    
-    def __init__(self, res, z_mag, y_mag, x_mag):  # TODO: electrostatic component!
-        '''Load magnetization in LLG-file format.
+    def __init__(self, res, magnitude):  # TODO: electrostatic component!
+        '''Constructor for a MagData object for storing magnetization data.
         Arguments:
-            filename - the name of the file where the data are stored
+            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:
-            None
-        
-        '''# TODO: Docstring
-        assert np.shape(x_mag) == np.shape(y_mag) == np.shape(z_mag), 'Dimensions do not match!'
-        self.res = res
-        self.dim = np.shape(x_mag)
-        self.magnitude = (z_mag, y_mag, x_mag)
+            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
         
     @classmethod
-    def load_from_llg(self, filename):
-        # TODO: Docstring
-        scale = 1.0E-9 / 1.0E-2  #from cm to nm
+    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, 
+        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)]
-        #Reshape in Python and Igor is different, Python fills rows first, Igor columns!
-        return MagData(res, z_mag, y_mag, x_mag)
+        return MagData(res, (z_mag, y_mag, x_mag))
 
 
-    def save_to_llg(self, filename='output.txt'):
-        '''Create homog. magnetization data, saved in a file with LLG-format.
+    def save_to_llg(self, filename='magdata_output.txt'):
+        '''Save magnetization data in a file with LLG-format.
         Arguments:
-            dim       - the dimensions of the grid, shape(y, x)
-            res       - the resolution of the grid 
-                        (real space distance between two points)
-            beta      - the angle of the magnetization
-            filename  - the name of the file in which to save the data
-            mag_shape - an array of shape dim, representing the shape of the magnetic object,
-                        a few are supplied in this module
+            filename - the name of the LLG-file in which to store the magnetization data
+                       (default: 'magdata_output.txt')
         Returns:
-            the the magnetic distribution as a 2D-boolean-array.
+            None
             
-        ''' # TODO: Renew Docstring
+        '''
         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 = np.reshape(xx,(-1))
         yy = np.reshape(yy,(-1))
         zz = np.reshape(zz,(-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('LLGFileCreator2D: %s\n' % filename.replace('.txt', ''))
@@ -78,44 +79,55 @@ class MagData:
             mag_file.writelines('\n'.join('   '.join('{:7.6e}'.format(cell) 
                                           for cell in row) for row in data) )
          
-                                                             
     def quiver_plot(self, axis='z', ax_slice=0):
-        # TODO: Docstring
-        assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as string).'
-        if axis == 'z':
-            mag_slice_1 = self.magnitude[2][ax_slice,...]
-            mag_slice_2 = self.magnitude[1][ax_slice,...]
-        elif axis == 'y':
-            mag_slice_1 = self.magnitude[2][:,ax_slice,:]
-            mag_slice_2 = self.magnitude[0][:,ax_slice,:]
-        elif axis == 'x':
-            mag_slice_1 = self.magnitude[1][...,ax_slice]
-            mag_slice_2 = self.magnitude[0][...,ax_slice]
+        '''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
+        Returns:
+            None
             
+        '''
+        assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as string).'
+        if axis == 'z':  # Slice of the xy-plane with z = ax_slice
+            mag_slice_1 = self.magnitude[2][ax_slice, ...]
+            mag_slice_2 = self.magnitude[1][ax_slice, ...]
+        elif axis == 'y': # Slice of the xz-plane with y = ax_slice
+            mag_slice_1 = self.magnitude[2][:, ax_slice, :]
+            mag_slice_2 = self.magnitude[0][:, ax_slice, :]
+        elif axis == 'x': # Slice of the yz-plane with x = ax_slice
+            mag_slice_1 = self.magnitude[1][..., ax_slice]
+            mag_slice_2 = self.magnitude[0][..., ax_slice]
+        # Plot the magnetization vectors:
         fig = plt.figure()
         fig.add_subplot(111, aspect='equal')
         plt.quiver(mag_slice_1, mag_slice_2, pivot='middle', angles='xy', scale_units='xy', 
                    scale=1, headwidth=6, headlength=7)
-                   
-                   
-    def quiver_plot_3D(self):
-        # TODO: Docstring 
-        res = self.res
-        dim = self.dim
-        fig = plt.figure()
-        ax = fig.add_subplot(111, projection='3d')        
-        ax.set_aspect("equal")
-
+                                     
+    def quiver_plot3d(self):  # XXX: Still buggy, use only for small distributions!
+        '''3D-Quiver-Plot of the magnetization as vectors.
+        Arguments:
+            None
+        Returns:
+            None
+            
+        '''
         class Arrow3D(FancyArrowPatch):
+            '''Class representing one magnetization vector.'''
             def __init__(self, xs, ys, zs, *args, **kwargs):
-                FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
+                FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs)
                 self._verts3d = xs, ys, zs
             def draw(self, renderer):
                 xs3d, ys3d, zs3d = self._verts3d
-                xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
-                self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
+                xs, ys = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)[:-1]
+                self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
                 FancyArrowPatch.draw(self, renderer)
-                
+        res = self.res
+        dim = self.dim
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection='3d')        
+        ax.set_aspect("equal")
+        # 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]
@@ -125,13 +137,14 @@ class MagData:
         x_mag = np.reshape(self.magnitude[2], (-1))        
         y_mag = np.reshape(self.magnitude[1], (-1))
         z_mag = np.reshape(self.magnitude[0], (-1))
-        
+        # Add every individual magnetization vector:
         for i in range(np.size(xx)):
             xs = [xx[i] - x_mag[i]*res/2, xx[i] + x_mag[i]*res/2]
             ys = [yy[i] - y_mag[i]*res/2, yy[i] + y_mag[i]*res/2]
             zs = [zz[i] - z_mag[i]*res/2, zz[i] + z_mag[i]*res/2]
             a = Arrow3D(xs, ys, zs, mutation_scale=10, lw=1, arrowstyle="-|>", color="k")
             ax.add_artist(a)
+        # Rescale the axes and show plot:
         ax.set_xlim3d(0, xx[-1]+res/2)
         ax.set_ylim3d(0, yy[-1]+res/2)
         ax.set_zlim3d(0, zz[-1]+res/2)

pyramid/phasemap.py

 # -*- coding: utf-8 -*-
-"""
-Created on Mon May 13 16:00:57 2013
+"""Class for creating objects to store phase maps."""
+
+
+import numpy as np
+import matplotlib.pyplot as plt
 
-@author: Jan
-"""
 
 class PhaseMap:
     
     '''An object storing magnetization data.'''
     
-    
-    def __init__(self, phase):  # TODO: more arguments?
-        '''Load magnetization in LLG-file format.
+    def __init__(self, res, phase):
+        '''Constructor for a MagData object for storing magnetization data.
         Arguments:
-            filename - the name of the file where the data are stored
+            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:
-            None
+            MagData object
         
-        '''# TODO: Docstring
+        '''
+        dim = np.shape(phase)
+        assert len(dim) == 2, 'Phasemap has to be 2-dimensional!'
+        self.res   = res
+        self.dim   = dim
         self.phase = phase
         
-        
-    def load_from_file(self, filename): # TODO: Implement
-        pass
-#        # TODO: Docstring
-#        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
-#        x_mag, y_mag, z_mag = [data[:, i].reshape(z_dim, y_dim, x_dim) 
-#                               for i in range(3,6)]
-#        #Reshape in Python and Igor is different, Python fills rows first, Igor columns!
-#        return MagData(res, z_mag, y_mag, x_mag)
-
+    @classmethod
+    def load_from_file(cls, filename):
+        '''Construct PhaseMap object from 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='    ', skiprows=2)
+        return PhaseMap(res, phase)
 
-    def save_to_file(self, filename='output.txt'): # TODO: Implement
-        pass
-#        '''Create homog. magnetization data, saved in a file with LLG-format.
-#        Arguments:
-#            dim       - the dimensions of the grid, shape(y, x)
-#            res       - the resolution of the grid 
-#                        (real space distance between two points)
-#            beta      - the angle of the magnetization
-#            filename  - the name of the file in which to save the data
-#            mag_shape - an array of shape dim, representing the shape of the magnetic object,
-#                        a few are supplied in this module
-#        Returns:
-#            the the magnetic distribution as a 2D-boolean-array.
-#            
-#        ''' # TODO: Renew Docstring
-#        dim = self.dim
-#        res = self.res * 1.0E-9 / 1.0E-2  # from nm to cm     
-#        
-#        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 = np.reshape(xx,(-1))
-#        yy = np.reshape(yy,(-1))
-#        zz = np.reshape(zz,(-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))
-#        
-#        data = np.array([xx, yy, zz, x_mag, y_mag, z_mag]).T
-#        with open(filename,'w') as mag_file:
-#            mag_file.write('LLGFileCreator2D: %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) )
+    def save_to_file(self, filename='phasemap_output.txt'):
+        '''Save magnetization data in a file with LLG-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='    ')        
 
-    def display_phase(self, res, title, axis=None):
+    def display(self, title='Phase Map', axis=None, cmap='gray'):
         '''Display the phasemap as a colormesh.
         Arguments:
-            phase - the phasemap that should be displayed
-            res   - the resolution of the phasemap
-            title - the title of the plot
+            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
             
         '''
+        # If no axis is specified, a new figure is created:
         if axis == None:
             fig = plt.figure()
             axis = fig.add_subplot(1,1,1, aspect='equal')
-        
-        im = plt.pcolormesh(phase, cmap='gray')
-    
-        ticks = axis.get_xticks()*res
+        im = plt.pcolormesh(self.phase, cmap=cmap)
+        # Set the axes ticks and labels:
+        ticks = axis.get_xticks()*self.res
         axis.set_xticklabels(ticks)
-        ticks = axis.get_yticks()*res
+        ticks = axis.get_yticks()*self.res
         axis.set_yticklabels(ticks)
-    
         axis.set_title(title)
         axis.set_xlabel('x-axis [nm]')
         axis.set_ylabel('y-axis [nm]')
-        
+        # Plot the phase map:
         fig = plt.gcf()
         fig.subplots_adjust(right=0.85)
         cbar_ax = fig.add_axes([0.9, 0.15, 0.02, 0.7])
         fig.colorbar(im, cax=cbar_ax)
-        
         plt.show()
\ No newline at end of file

pyramid/phasemapper.py

@@ -3,148 +3,131 @@
 
 
 import numpy as np
-import matplotlib.pyplot as plt
-import pyramid.projector as pj
-from numpy import pi
-from phasemap import PhaseMap
 
 
+# Physical constants
 PHI_0 = 2067.83    # magnetic flux in T*nm²
-H_BAR = 6.626E-34  # TODO: unit
-M_E   = 9.109E-31  # TODO: unit
-Q_E   = 1.602E-19  # TODO: unit
-C     = 2.998E8    # TODO: unit
+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 fourier_space(mag_data, b_0=1, padding=0):
-    '''Calculate phasemap from magnetization data (fourier approach).
+def phase_mag_fourier(res, projection, b_0=1, padding=0):
+    '''Calculate phasemap from magnetization data (Fourier space approach).
     Arguments:
-        mag_data - MagDataLLG object (from magneticimaging.dataloader) storing
-                   the filename, coordinates and magnetization in x, y and z
-        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
-        v_0      - average potential of the sample in V (default: 0)
-        v_acc    - acceleration voltage of the microscop in V (default: 30000)
+        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
     
-    '''    
-    res  = mag_data.res
-    z_dim, y_dim, x_dim = mag_data.dim
-    z_mag, y_mag, x_mag = pj.simple_axis_projection(mag_data)#mag_data.magnitude  
-    
-    # TODO: interpolation (redefine xDim,yDim,zDim) thickness ramp
-    
-    x_pad = x_dim/2 * padding
-    y_pad = y_dim/2 * padding
-    x_mag_big = np.zeros(((1 + padding) * y_dim, (1 + padding) * x_dim))
-    y_mag_big = np.zeros(((1 + padding) * y_dim, (1 + padding) * x_dim))
-    # TODO: padding so that x_dim and y_dim = 2^n
-    x_mag_big[y_pad : y_pad + y_dim, x_pad : x_pad + x_dim] = x_mag
-    y_mag_big[y_pad : y_pad + y_dim, x_pad : x_pad + x_dim] = y_mag
-    
-    x_mag_fft = np.fft.fftshift(np.fft.rfft2(x_mag_big), axes=0)
-    y_mag_fft = np.fft.fftshift(np.fft.rfft2(y_mag_big), axes=0)
+    '''
+    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
-    x = np.linspace(      0, nyq/2, x_mag_fft.shape[1])
-    y = np.linspace( -nyq/2, nyq/2, x_mag_fft.shape[0]+1)[:-1]
-    xx, yy = np.meshgrid(x, y)
-                         
-    phase_fft = (1j * res * b_0) / (2 * PHI_0) * ((x_mag_fft * yy - y_mag_fft * xx) 
-                                           / (xx ** 2 + yy ** 2 + 1e-18))
+    u = np.linspace(      0, nyq/2, u_mag_fft.shape[1])
+    v = np.linspace( -nyq/2, nyq/2, u_mag_fft.shape[0]+1)[:-1]
+    uu, vv = np.meshgrid(u, v)
+    coeff = (1j * res * b_0) / (2 * PHI_0)              
+    phase_fft = coeff * (u_mag_fft*vv - v_mag_fft*uu) / (uu**2 + 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[y_pad : y_pad + y_dim, x_pad : x_pad + x_dim]
-    
-    # TODO: Electrostatic Component
-
-    return PhaseMap(phase)
-      
+    phase = phase_big[v_pad : v_pad + v_dim, u_pad : u_pad + u_dim]
+    return phase
       
-def phi_pixel(method, n, m, res, b_0):  # TODO: rename
-    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            
-        coeff = -b_0 * res**2 / ( 2 * PHI_0 )
-        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':
-        coeff = - b_0 * res**2 / ( 2 * PHI_0 )
-        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
 
-
-def real_space(mag_data, method, b_0=1, jacobi=None):
+def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
     '''Calculate phasemap from magnetization data (real space approach).
     Arguments:
-        mag_data - MagDataLLG object (from magneticimaging.dataloader) storing the filename, 
-                   coordinates and magnetization in x, y and z
+        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)
+        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
-        
-    '''    
-    # TODO: Expand docstring!  
-
-    res  = mag_data.res
-    z_dim, y_dim, x_dim = mag_data.dim
-    z_mag, y_mag, x_mag = mag_data.magnitude
     
-    z_mag, y_mag, x_mag = pj.simple_axis_projection(mag_data)
-    
-    beta = np.arctan2(y_mag, x_mag)
-
-    mag = np.hypot(x_mag, y_mag) 
+    '''# TODO: Docstring!
+    def phi_lookup(method, n, m, res, b_0):  # TODO: rename
+        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            
+            coeff = -b_0 * res**2 / ( 2 * PHI_0 )
+            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':
+            coeff = - b_0 * res**2 / ( 2 * PHI_0 )
+            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
+    
+    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) 
                 
     '''CREATE COORDINATE GRIDS'''
-    x = np.linspace(0,(x_dim-1),num=x_dim)
-    y = np.linspace(0,(y_dim-1),num=y_dim)
-    xx, yy = np.meshgrid(x,y)
+    u = np.linspace(0,(u_dim-1),num=u_dim)
+    v = np.linspace(0,(v_dim-1),num=v_dim)
+    uu, vv = np.meshgrid(u,v)
      
-    x_big = np.linspace(-(x_dim-1), x_dim-1, num=2*x_dim-1)
-    y_big = np.linspace(-(y_dim-1), y_dim-1, num=2*y_dim-1)
-    xx_big, yy_big = np.meshgrid(x_big, y_big)
+    u_lookup = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
+    v_lookup = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
+    uu_lookup, vv_lookup = np.meshgrid(u_lookup, v_lookup)
     
-    phi_cos = phi_pixel(method, xx_big, yy_big, res, b_0)
-    phi_sin = phi_pixel(method, yy_big, xx_big, res, b_0)
+    phi_cos = phi_lookup(method, uu_lookup, vv_lookup, res, b_0)
+    phi_sin = phi_lookup(method, vv_lookup, uu_lookup, res, b_0)
             
     def phi_mag(i, j):  # TODO: rename
-        return (np.cos(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i]
-               -np.sin(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i])
+        return (np.cos(beta[j,i])*phi_cos[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+               -np.sin(beta[j,i])*phi_sin[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
                                           
     def phi_mag_deriv(i, j):  # TODO: rename
-        return -(np.sin(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i]
-                +np.cos(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i])
+        return -(np.sin(beta[j,i])*phi_cos[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+                +np.cos(beta[j,i])*phi_sin[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
                                            
     def phi_mag_fd(i, j, h):  # TODO: rename
         return ((np.cos(beta[j,i]+h) - np.cos(beta[j,i])) / h 
-                      * phi_cos[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i]
+                      * phi_cos[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
                -(np.sin(beta[j,i]+h) - np.sin(beta[j,i])) / h 
-                      * phi_sin[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i])
+                      * phi_sin[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
     
     '''CALCULATE THE PHASE'''
-    phase = np.zeros((y_dim, x_dim))
+    phase = np.zeros((v_dim, u_dim))
     
     # TODO: only iterate over pixels that have a magn. > threshold (first >0)
     if jacobi is not None:
         jacobi_fd = jacobi.copy()
     h = 0.0001
     
-    for j in range(y_dim):
-        for i in range(x_dim):
+    for j in range(v_dim):
+        for i in range(u_dim):
             #if (mag[j,i] != 0 ):#or jacobi is not None): # TODO: same result with or without?
                 phi_mag_cache = phi_mag(i, j)
                 phase += mag[j,i] * phi_mag_cache
                 if jacobi is not None:
-                    jacobi[:,i+x_dim*j] = phi_mag_cache.reshape(-1)
-                    jacobi[:,x_dim*y_dim+i+x_dim*j] = (mag[j,i]*phi_mag_deriv(i,j)).reshape(-1)
+                    jacobi[:,i+u_dim*j] = phi_mag_cache.reshape(-1)
+                    jacobi[:,u_dim*v_dim+i+u_dim*j] = (mag[j,i]*phi_mag_deriv(i,j)).reshape(-1)
                     
-                    jacobi_fd[:,i+x_dim*j] = phi_mag_cache.reshape(-1)
-                    jacobi_fd[:,x_dim*y_dim+i+x_dim*j] = (mag[j,i]*phi_mag_fd(i,j,h)).reshape(-1)  
+                    jacobi_fd[:,i+u_dim*j] = phi_mag_cache.reshape(-1)
+                    jacobi_fd[:,u_dim*v_dim+i+u_dim*j] = (mag[j,i]*phi_mag_fd(i,j,h)).reshape(-1)  
                     
                     
     
@@ -155,23 +138,20 @@ def real_space(mag_data, method, b_0=1, jacobi=None):
     return phase
     
 
-def phase_elec(mag_data, b_0=1, v_0=0, v_acc=30000):
-    # TODO: Docstring
-
-    # TODO: Delete    
-#    import pdb; pdb.set_trace()
-    
-    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
+#def phase_elec(mag_data, v_0=0, v_acc=30000):
+#    # TODO: Docstring
+#
+#    res  = mag_data.res
+#    z_dim, y_dim, x_dim = mag_data.dim
+#    z_mag, y_mag, x_mag = mag_data.magnitude  
+#    
+#    phase = np.logical_or(x_mag, y_mag, z_mag)    
+#    
+#    lam = H_BAR / np.sqrt(2 * M_E * v_acc * (1 + Q_E*v_acc / (2*M_E*C**2)))
+#    
+#    Ce = (2*pi*Q_E/lam * (Q_E*v_acc +   M_E*C**2) /
+#            (Q_E*v_acc * (Q_E*v_acc + 2*M_E*C**2)))
+#    
+#    phase *= res * v_0 * Ce
+#    
+#    return phase
\ No newline at end of file
-    
-    return phase
\ No newline at end of file

pyramid/projector.py

 # -*- coding: utf-8 -*-
-"""
-Created on Tue May 14 14:44:36 2013
+"""Planar projection of the magnetization distribution of a MagData object."""
 
-@author: Jan
-"""# TODO: Docstring
 
-from magdata import MagData
+from pyramid.magdata import MagData
 
-def simple_axis_projection(mag_data, axis=0):
-    # TODO: assert isinstance(mag_data, MagData), 'Object is no instance of MagData!'
-    return (mag_data.magnitude[0].mean(axis),  # x_mag
-            mag_data.magnitude[1].mean(axis),  # y_mag
-            mag_data.magnitude[2].mean(axis))  # z_mag
+
+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].mean(0) * mag_data.dim[0],  # y_mag -> v_mag
+                      mag_data.magnitude[2].mean(0) * mag_data.dim[0])  # x_mag -> u_mag
+    elif axis == 'y':
+        projection = (mag_data.magnitude[0].mean(1) * mag_data.dim[1],  # z_mag -> v_mag
+                      mag_data.magnitude[2].mean(1) * mag_data.dim[1])  # x_mag -> u_mag
+    elif axis == 'x':
+        projection = (mag_data.magnitude[0].mean(2) * mag_data.dim[2],  # y_mag -> v_mag
+                      mag_data.magnitude[1].mean(2) * mag_data.dim[2])  # x_mag -> u_mag
+    return projection
     
     
-# TODO: proper projection algorithm with two angles and such!
+# TODO: proper projection algorithm with two angles and such!
\ No newline at end of file
+# CAUTION: the res for the projection does not have to be the res of the 3D-magnetization!
+# Just for a first approach
\ No newline at end of file

pyramid/reconstructor.py

@@ -5,3 +5,4 @@ Created on Tue May 14 15:19:33 2013
 @author: Jan
 """
 
+# TODO: Implement
\ No newline at end of file

scripts/compare_methods.py

 # -*- coding: utf-8 -*-
-"""
-Created on Wed Apr 03 11:15:38 2013
+"""Compare the different methods to create phase maps."""
 
-@author: Jan
-"""
 
-import numpy as np
-import matplotlib.pyplot as plt
-import pyramid.magcreator as mc
-import pyramid.dataloader as dl
-import pyramid.phasemap as pm
-import pyramid.holoimage as hi
-import pyramid.analytic as an
 import time
 import pdb, traceback, 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 phase_from_mag():
     '''Calculate and display the phase map from a given magnetization.
@@ -25,89 +24,53 @@ def phase_from_mag():
         None
     
     '''
-    # TODO: Input via GUI
-    b_0 = 1.0  # in T
-    v_0 = 0  # TODO: units?
-    v_acc = 30000  # in V
+    # Input parameters:
+    b_0     =  1    # in T
+    res     = 10.0  # in nm
+    beta    = pi/4
     padding = 20
     density = 10
-    
-    dim = (50, 50)  # in px (y,x)
-    res = 10.0  # in nm
-    beta = pi/4
-    
-    center = (24, 24)  # in px (y,x) index starts with 0!
-    width  = (25, 25)  # in px (y,x)
-    radius = 12.5
-    
-    filename = '../output/output.txt'
+    dim = (1, 50, 50)  # in px (z, y, x)    
+    # Create magnetic shape:
     geometry = 'slab'        
-    
     if geometry == 'slab':
-        mag_shape = mc.slab(dim, center, width)
+        center = (0, 24, 24)  # in px (z, y, x) index starts with 0!
+        width  = (1, 25, 25)  # in px (z, y, x)
+        mag_shape = mc.Shapes.slab(dim, center, width)
         phase_ana = an.phasemap_slab(dim, res, beta, center, width, b_0)
-    elif geometry == 'slab':
-        mag_shape = mc.disc(dim, center, radius)
+    elif geometry == 'disc':
+        radius = 12.5  # in px 
+        height =  1    # in px
+        mag_shape = mc.Shapes.disc(dim, center, radius, height)
         phase_ana = an.phasemap_disc(dim, res, beta, center, radius, b_0)
-    
-    holo_ana  = hi.holo_image(phase_ana, res, density)
-    display_combined(phase_ana, res, holo_ana, 'Analytic Solution')
-    
-    '''CREATE MAGNETIC DISTRIBUTION'''
-    mc.create_hom_mag(dim, res, beta, mag_shape, filename)
-    
-    '''LOAD MAGNETIC DISTRIBUTION'''
-    mag_data = dl.MagDataLLG(filename)
-    
-    # TODO: get it to work:
-    phase_el = pm.phase_elec(mag_data, v_0, v_acc)    
-    
-    '''NUMERICAL SOLUTION'''
-    # numerical solution Fourier Space:
-    tic = time.clock()
-    phase_fft = pm.fourier_space(mag_data, b_0, padding)   
-    toc = time.clock()
-    print 'Time for Fourier Space Approach:     ' + str(toc - tic)
-    holo_fft = hi.holo_image(phase_fft, mag_data.res, density)
-    display_combined(phase_fft, mag_data.res, holo_fft, 'Fourier Space Approach')
-    
-    # numerical solution Real Space (Slab):
-    tic = time.clock()
-    phase_slab = pm.real_space(mag_data, 'slab', b_0)
-    toc = time.clock()
-    print 'Time for Real Space Approach (Slab): ' + str(toc - tic)
-    holo_slab = hi.holo_image(phase_slab, mag_data.res, density)
-    display_combined(phase_slab, mag_data.res, holo_slab, 'Real Space Approach (Slab)')
-    
-    # numerical solution Real Space (Disc):
-    tic = time.clock()
-    phase_disc = pm.real_space(mag_data, 'disc', b_0)
-    toc = time.clock()
-    print 'Time for Real Space Approach (Disc): ' + str(toc - tic)
-    holo_disc = hi.holo_image(phase_disc, mag_data.res, density)
-    display_combined(phase_disc, mag_data.res, holo_disc, 'Real Space Approach (Disc)')
-    
-    '''DIFFERENCES'''
-    diff_fft  = phase_fft  - phase_ana 
-    diff_slab = phase_slab - phase_ana
-    diff_disc = phase_disc - phase_ana
-    rms_fft  = np.sqrt((diff_fft**2).mean())
-    rms_slab = np.sqrt((diff_slab**2).mean())
-    rms_disc = np.sqrt((diff_disc**2).mean())
-    pm.display_phase(diff_fft,  res, 'FFT - Analytic (RMS = ' + '{:3.2e}'.format(rms_fft) + ')')
-    pm.display_phase(diff_slab, res, 'Slab - Analytic (RMS = ' +'{:3.2e}'.format(rms_slab) + ')')
-    pm.display_phase(diff_disc, res, 'Disc - Analytic (RMS = ' + '{:3.2e}'.format(rms_disc) + ')')
- 
-    
-def display_combined(phase, res, holo, title):
-    fig = plt.figure(figsize=(14, 7))
-    fig.suptitle(title, fontsize=20)
-    
-    holo_axis = fig.add_subplot(1,2,1, aspect='equal')
-    hi.display_holo(holo, 'Holography Image', holo_axis)
-    
-    phase_axis = fig.add_subplot(1,2,2, aspect='equal')
-    pm.display_phase(phase, res, 'Phasemap', phase_axis)
+    # Project the magnetization data:    
+    mag_data = MagData(res, mc.create_mag_dist(mag_shape, beta))    
+    projection = pj.simple_axis_projection(mag_data)
+    # Construct phase maps:
+    phase_map_ana  = PhaseMap(res, phase_ana)
+    phase_map_fft  = PhaseMap(res, pm.phase_mag_fourier(res, projection, b_0, padding))
+    phase_map_slab = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab', b_0))
+    phase_map_disc = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc', b_0))
+    # 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)')
+    # Display all phase maps:
+    phase_map_ana.display('Analytic Solution')
+    phase_map_fft.display('Fourier Space')
+    phase_map_slab.display('Real Space (Slab)')
+    phase_map_disc.display('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  = phase_map_diff_fft.phase
+    RMS_slab = phase_map_diff_slab.phase
+    RMS_disc = phase_map_diff_disc.phase
+    phase_map_diff_fft.display('Fourier Space (RMS = {})'.format(np.std(RMS_fft)))
+    phase_map_diff_slab.display('Real Space (Slab) (RMS = {})'.format(np.std(RMS_slab)))
+    phase_map_diff_disc.display('Real Space (Disc) (RMS = {})'.format(np.std(RMS_disc)))
     
     
 if __name__ == "__main__":

scripts/compare_pixel_fields.py

 # -*- coding: utf-8 -*-
-"""
-Created on Wed Apr 03 11:15:38 2013
+"""Compare the phase map of one pixel for different real space approaches."""
 
-@author: Jan
-"""
 
-import numpy as np
-import pyramid.phasemap as pm
 import pdb, traceback, 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 from a given magnetization.
+    '''Calculate and display the phase map for different real space approaches.
     Arguments:
         None
     Returns:
         None
     
     '''
-    # TODO: Input via GUI
-    b_0 = 1.0  # in T
-    res = 10.0   
-    dim = (10, 10)
-    
-    x_big = np.linspace(-(dim[1]-1), dim[1]-1, num=2*dim[1]-1)
-    y_big = np.linspace(-(dim[0]-1), dim[0]-1, num=2*dim[0]-1)
-    xx_big, yy_big = np.meshgrid(x_big, y_big)    
-    
-    phi_cos_real_slab = pm.phi_pixel('slab', xx_big, yy_big, res, b_0)
-    pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
-    phi_cos_real_disc = pm.phi_pixel('disc', xx_big, yy_big, res, b_0)
-    pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
-    phi_cos_diff = phi_cos_real_disc - phi_cos_real_slab
-    pm.display_phase(phi_cos_diff, res, 'Phase of one Pixel-Disc (Cos - Part)')
+    # Input parameters:    
+    res   = 10.0  # in nm
+    beta  = pi/2  # in rad
+    dim   = (1, 11, 11)
+    pixel = (0,  5,  5) 
+    # Create magnetic data, project it, get the phase map and display the holography image:    
+    mag_data   = MagData(res, mc.create_mag_dist(mc.Shapes.single_pixel(dim, pixel), beta)) 
+    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__":

scripts/create_alternating_filaments.py

+# -*- coding: utf-8 -*-
+"""Create magnetic distribution of alternating filaments"""
+
+
+import pdb, traceback, 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
+    beta = pi/2       
+    spacing = 5    
+    # Create lists for magnetic objects:
+    count = int((dim[1]-1) / spacing) + 1
+    mag_shape_list = np.zeros((count,) + dim)
+    beta_list      = np.zeros(count)
+    for i in range(count):  
+        pos = i * spacing
+        mag_shape_list[i] = mc.Shapes.filament(dim, (0, pos))
+        beta_list[i] = (1-2*(int(pos/spacing)%2)) * beta
+    # Create magnetic distribution
+    magnitude = mc.create_mag_dist_comb(mag_shape_list, beta_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)
\ No newline at end of file

scripts/create_logo.py

 # -*- coding: utf-8 -*-
-"""
-Created on Wed Apr 03 11:15:38 2013
+"""Create the Pyramid-Logo."""
 
-@author: Jan
-"""
 
-import pyramid.magcreator as mc
-import pyramid.dataloader as dl
-import pyramid.phasemap as pm
-import pyramid.holoimage as hi
 import pdb, traceback, 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 phase map from a given magnetization.
+    '''Calculate and display the Pyramid-Logo.
     Arguments:
         None
     Returns:
         None
     
     '''
-    filename = '../output/mag_distr_logo.txt'
-    b_0 = 1.0  # in T
+    # Input parameters:    
     res = 10.0  # in nm
     beta = pi/2  # in rad
     density = 10
-    dim = (128, 128)    
-    
-    x = range(dim[1])
-    y = range(dim[0])
+    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[0])
-    left   = (yy <= 0.75/0.5 * dim[0]/dim[1] * xx)
+    bottom = (yy >= 0.25*dim[1])
+    left   = (yy <= 0.75/0.5 * dim[1]/dim[2] * xx)
     right  = np.fliplr(left)
-    mag_shape = np.logical_and(np.logical_and(left, right), bottom)
-    
-    mc.create_hom_mag(dim, res, beta, mag_shape, filename)
-    mag_data = dl.MagDataLLG(filename)
-    phase= pm.real_space(mag_data, 'slab', b_0)  
-    holo = hi.holo_image(phase, mag_data.res, density)
-    hi.display_holo(holo, 'PYRAMID - LOGO')
+    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, beta)) 
+    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')
     
     
 if __name__ == "__main__":
-#    parser = argparse.ArgumentParser(description='Create the PYRAMID logo')
-#    parser.add_argument('-d','--dimensions', help='Logo dimensions.', required=False)
-#    args = parser.parse_args()
-#    if args.dimensions is None:
-#        args.dimensions = (128,128)
     try:
         create_logo()
     except:

scripts/create_random_distribution.py

 # -*- coding: utf-8 -*-
-"""
-Created on Mon May 13 13:05:40 2013
-
-@author: Jan
-"""
+"""Create random magnetic distributions."""
 
 import random as rnd
-import numpy as np
-import pyramid.magcreator as mc
 import pdb, traceback, 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_distribution():
+    '''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)
-    
-    mag_shape_list = np.zeros((count, dim[0], dim[1], dim[2]))
+    rnd.seed(12)
+    # Create lists for magnetic objects:
+    mag_shape_list = np.zeros((count,) + dim)
     beta_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.shape_single_pixel(dim, pixel)
+        mag_shape_list[i,...] = mc.Shapes.single_pixel(dim, pixel)
         beta_list[i] = 2*pi*rnd.random()
         magnitude_list[i] = 1#rnd.random()
-        
-    mag_data = mc.create_mag_dist(dim, res, mag_shape_list, beta_list, magnitude_list)
+    # Create magnetic distribution
+    magnitude = mc.create_mag_dist_comb(mag_shape_list, beta_list, magnitude_list) 
+    mag_data = MagData(res, magnitude)
     mag_data.quiver_plot()
-    #mag_data.quiver_plot_3D()
     mag_data.save_to_llg('../output/mag_dist_random_pixel.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__":

scripts/create_sample.py

 # -*- coding: utf-8 -*-
-"""
-Created on Wed Apr 03 11:15:38 2013
+"""Create magnetic distributions with simple geometries."""
 
-@author: Jan
-"""
 
-import pyramid.magcreator as mc
 import pdb, traceback, sys
 from numpy import pi
 
+import pyramid.magcreator as mc
+from pyramid.magdata import MagData
+
 
 def create_sample():
-    '''Calculate and display the phase map from a given magnetization.
+    '''Calculate, display and save simple magnetic distributions to file.
     Arguments:
         None
     Returns:
         None
     
     '''
-    
-    # TODO: Input via GUI
+    # Input parameters:
     key = 'slab'
-    
-    filename = '../output/mag_distr_' + key + '.txt'    
-    dim = (50, 50)  # in px (y,x)
-    res = 1.0  # in nm
-    beta = pi/4
-    plot_mag_distr = True    
-    
-    center = (24, 24)  # in px (y,x) index starts with 0!
-    width  = (25, 25)  # in px (y,x)
-    radius = 12.5  # in px
-    pos = 24  # in px
-    spacing = 5  # in px
-    x_or_y = 'y'  
-    pixel = (24, 24) # in px    
-    
+    filename = '../output/mag_dist_' + key + '.txt'    
+    dim = (1, 128, 128)  # in px (z, y, x)
+    res = 10.0  # in nm
+    beta = pi/4   
+    # Geometry parameters:    
+    center = (0, 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.slab(dim, center, width)
+        mag_shape = mc.Shapes.slab(dim, center, width)
     elif key == 'disc':
-        mag_shape = mc.disc(dim, center, radius)
+        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.filament(dim, pos, x_or_y)
-    elif key == 'alternating_filaments':
-        mag_shape = mc.alternating_filaments(dim, spacing, x_or_y)
+        mag_shape = mc.Shapes.filament(dim, pos)
     elif key == 'pixel':
-        mag_shape = mc.single_pixel(dim, pixel)
-    
-    mc.create_hom_mag(dim, res, beta, mag_shape, filename, plot_mag_distr)
-    
+        mag_shape = mc.Shapes.single_pixel(dim, pixel)
+    # Create magnetic distribution
+    magnitude = mc.create_mag_dist(mag_shape, beta) 
+    mag_data = MagData(res, magnitude)
+    mag_data.quiver_plot()
+    mag_data.save_to_llg(filename)
+
     
 if __name__ == "__main__":
     try: