From 78d70d7c1632d4aa7ebd11d69cda670b30dcd6df Mon Sep 17 00:00:00 2001
From: Jan Caron <caron@fz-juelich.de>
Date: Mon, 6 May 2013 11:32:30 +0200
Subject: [PATCH] scripts: Splitted the main in several smaller scripts test:
 changed the dataloader TestCase pyramex: first experiments with Cython and
 SCons dataloader: changed the order in which the dimensions are saved,       
  z_len is no longer part of the magnetization (this was the reason why       
  we had to divide by res, which is equal to z_len for calculating coeff)
 magcreator: shapes are now separate functions, mag_shape is now an argument
 for         one unified create_hom_mag-function phasemap: unification of the
 phasemap-functions in real space (just the pixel-         field is calculated
 different now)

---
 SConstruct                         |  25 +++++
 scripts/compare_methods.py         | 127 +++++++--------------
 scripts/compare_pixel_fields.py    |  44 ++++++++
 scripts/create_logo.py             |  26 ++++-
 scripts/create_sample.py           |  59 ++++++++++
 scripts/create_samples.py          | 172 -----------------------------
 scripts/cython_example.py          |  55 +++++++++
 scripts/get_jacobi.py              | 138 +++--------------------
 scripts/main.py                    | 172 -----------------------------
 src/pyramex/.sconsign.dblite       | Bin 0 -> 2379 bytes
 src/pyramex/SConstruct             |  16 +++
 src/pyramex/SConstructBACKUP       |  70 ++++++++++++
 src/pyramex/c1.pyx                 |  13 +++
 src/pyramex/c2.pyx                 |  16 +++
 src/pyramex/c3.pyx                 |  23 ++++
 src/pyramex/hello.pyx              |   2 +
 src/pyramex/setup.py               |  36 ++++++
 src/pyramid/dataloader.py          |  12 +-
 src/pyramid/magcreator.py          | 107 +++++-------------
 src/pyramid/phasemap.py            | 171 +++++++++-------------------
 test/run_tests.py                  |  74 +++++++++++++
 test/test_compliance.py            |  82 ++++++++++++++
 test/test_dataloader.py            |  71 ++++++++++--
 test/test_dataloader/test_data.txt |  17 +++
 test/testsuite.py                  |  29 -----
 25 files changed, 750 insertions(+), 807 deletions(-)
 create mode 100644 SConstruct
 create mode 100644 scripts/compare_pixel_fields.py
 create mode 100644 scripts/create_sample.py
 delete mode 100644 scripts/create_samples.py
 create mode 100644 scripts/cython_example.py
 delete mode 100644 scripts/main.py
 create mode 100644 src/pyramex/.sconsign.dblite
 create mode 100644 src/pyramex/SConstruct
 create mode 100644 src/pyramex/SConstructBACKUP
 create mode 100644 src/pyramex/c1.pyx
 create mode 100644 src/pyramex/c2.pyx
 create mode 100644 src/pyramex/c3.pyx
 create mode 100644 src/pyramex/hello.pyx
 create mode 100644 src/pyramex/setup.py
 create mode 100644 test/run_tests.py
 create mode 100644 test/test_compliance.py
 create mode 100644 test/test_dataloader/test_data.txt
 delete mode 100644 test/testsuite.py

diff --git a/SConstruct b/SConstruct
new file mode 100644
index 0000000..2b285f9
--- /dev/null
+++ b/SConstruct
@@ -0,0 +1,25 @@
+import os
+import distutils.sysconfig
+import sys
+
+PYTHON_LIBPATH = distutils.sysconfig.get_python_lib()
+PYTHON_INCPATH = distutils.sysconfig.get_python_inc()
+PYTHON_LIBRARY = "python" + sys.version[:3]
+
+print 'PYTHON_LIBPATH: ' + PYTHON_LIBPATH
+print 'PYTHON_INCPATH: ' + PYTHON_INCPATH
+print 'PYTHON_LIBRARY: ' + PYTHON_LIBRARY
+
+env = Environment(ENV = os.environ)
+
+if ARGUMENTS.get('VERBOSE') != '1':
+    env['CCCOMSTR'] = "Compiling $TARGET"
+    env['LINKCOMSTR'] = "Linking $TARGET"
+
+Progress('$TARGET\r', overwrite=True)
+
+env.AppendUnique(LIBPATH=[PYTHON_LIBPATH], CPPPATH=[PYTHON_INCPATH])
+
+env.Program('hello', 'hellotest.c', CPPPATH='.')
+
+env.Decider('MD5-timestamp')
\ No newline at end of file
diff --git a/scripts/compare_methods.py b/scripts/compare_methods.py
index cd8a8cd..6f76bd1 100644
--- a/scripts/compare_methods.py
+++ b/scripts/compare_methods.py
@@ -25,75 +25,42 @@ def phase_from_mag():
         None
     
     '''
-    
-    '''INPUT'''
     # TODO: Input via GUI
-    filename = '../output/output.txt'
-    b_0 = 10.0  # in T
+    b_0 = 1.0  # in T
     v_0 = 0  # TODO: units?
     v_acc = 30000  # in V
     padding = 20
-    density = 100
+    density = 10
     
     dim = (50, 50)  # in px (y,x)
-    res = 1.0  # in nm
+    res = 10.0  # in nm
     beta = pi/4
     
-    plot_mag_distr = True
-    
-    # Slab:
-    shape_fun = mc.slab
     center = (24, 24)  # in px (y,x) index starts with 0!
     width  = (25, 25)  # in px (y,x)
-    params = (center, width)
-#    # Disc:
-#    shape_fun = mc.disc
-#    center = (4, 4)  # in px (y,x)
-#    radius = 2.5
-#    params = (center, radius)
-#    # Filament:
-#    shape_fun = mc.filament
-#    pos = 5
-#    x_or_y = 'y'
-#    params = (pos, x_or_y)
-#    # Alternating Filaments:
-#    shape_fun = mc.alternating_filaments
-#    spacing = 5
-#    x_or_y = 'y'
-#    params = (spacing, x_or_y)
-#    # Single Pixel:
-#    shape_fun = mc.single_pixel
-#    pixel = (5, 5)
-#    params = pixel
-    
-#    '''CREATE LOGO'''    
-#    mc.create_logo(128, res, beta, filename, plot_mag_distr)
-#    mag_data = dl.MagDataLLG(filename)
-#    phase= pm.real_space_slab(mag_data, b_0)  
-#    holo = hi.holo_image(phase, mag_data.res, density)
-#    hi.display_holo(holo, '')
+    radius = 12.5
+    
+    filename = '../output/output.txt'
+    geometry = 'slab'        
+    
+    if geometry == 'slab':
+        mag_shape = mc.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)
+        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, shape_fun, params,
-                      filename, plot_mag_distr)
+    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=1, v_acc=200000)    
-    
-    
-    phi_cos_real_slab = np.zeros((2*dim[0]-1, 2*dim[1]-1))
-    pm.real_space_slab(mag_data, b_0, pixel_map=phi_cos_real_slab)
-    pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
-    phi_cos_real_disc = np.zeros((2*dim[0]-1, 2*dim[1]-1))
-    pm.real_space_disc(mag_data, b_0, pixel_map=phi_cos_real_disc)
-    pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
-    phi_cos_diff = phi_cos_real_slab - phi_cos_real_disc
-    pm.display_phase(phi_cos_diff, mag_data.res, 'Difference: One Pixel Slab - Disc')
-    
+    phase_el = pm.phase_elec(mag_data, v_0, v_acc)    
     
     '''NUMERICAL SOLUTION'''
     # numerical solution Fourier Space:
@@ -102,54 +69,34 @@ def phase_from_mag():
     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')
+    display_combined(phase_fft, mag_data.res, holo_fft, 'Fourier Space Approach')
+    
     # numerical solution Real Space (Slab):
-    jacobi = np.zeros((2*dim[0]*dim[1], dim[0]*dim[1]))
     tic = time.clock()
-    phase_real_slab = pm.real_space_slab(mag_data, b_0, jacobi=None)
+    phase_slab = pm.real_space(mag_data, 'slab', b_0)
     toc = time.clock()
-    np.savetxt('../output/jacobi.npy', jacobi)
     print 'Time for Real Space Approach (Slab): ' + str(toc - tic)
-    holo_real_slab = hi.holo_image(phase_real_slab, mag_data.res, density)
-    display_combined(phase_real_slab, mag_data.res, holo_real_slab, 
-                     'Real Space Approach (Slab)')
+    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_real_disc = pm.real_space_disc(mag_data, b_0)
+    phase_disc = pm.real_space(mag_data, 'disc', b_0)
     toc = time.clock()
     print 'Time for Real Space Approach (Disc): ' + str(toc - tic)
-    holo_real_disc = hi.holo_image(phase_real_disc, mag_data.res, density)
-    display_combined(phase_real_disc, mag_data.res, holo_real_disc,
-                     'Real Space Approach (Disc)')
-    
-    '''ANALYTIC SOLUTION'''
-    # analytic solution slab:
-    phase = an.phasemap_slab(dim, res, beta, center, width, b_0)
-    holo  = hi.holo_image(phase, res, density)
-    display_combined(phase, res, holo, 'Analytic Solution - Slab')
-#    # analytic solution disc:
-#    phase = an.phasemap_disc(dim, res, beta, center, radius, b_0)
-#    holo  = hi.holo_image(phase, res, density)
-#    display_combined(phase, res, holo, 'Analytic Solution - Disc')
-#    # analytic solution sphere:
-#    phase = an.phasemap_sphere(dim, res, beta, center, radius, b_0)
-#    holo  = hi.holo_image(phase, res, density)
-#    display_combined(phase, res, holo, 'Analytic Solution - Sphere')
-
+    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_to_ana       = phase_fft       - phase 
-    diff_real_slab_to_ana = phase_real_slab - phase
-    diff_real_disc_to_ana = phase_real_disc - phase
-    diff_slab_to_disc     = phase_real_disc - phase_real_slab
-    pm.display_phase(diff_fft_to_ana,       res, 'Difference: FFT - Analytic')
-    pm.display_phase(diff_real_slab_to_ana, res, 'Difference: Slab - Analytic')
-    pm.display_phase(diff_real_disc_to_ana, res, 'Difference: Disc - Analytic')
-    pm.display_phase(diff_slab_to_disc,     res, 'Difference: Disc - Slab')
-    
-    # TODO: Delete
-#    import pdb; pdb.set_trace()
+    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):
diff --git a/scripts/compare_pixel_fields.py b/scripts/compare_pixel_fields.py
new file mode 100644
index 0000000..8d0f175
--- /dev/null
+++ b/scripts/compare_pixel_fields.py
@@ -0,0 +1,44 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr 03 11:15:38 2013
+
+@author: Jan
+"""
+
+import numpy as np
+import pyramid.phasemap as pm
+import pdb, traceback, sys
+
+
+def compare_pixel_fields():
+    '''Calculate and display the phase map from a given magnetization.
+    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)')
+    
+    
+if __name__ == "__main__":
+    try:
+        compare_pixel_fields()
+    except:
+        type, value, tb = sys.exc_info()
+        traceback.print_exc()
+        pdb.post_mortem(tb)
\ No newline at end of file
diff --git a/scripts/create_logo.py b/scripts/create_logo.py
index 0b3c0e9..0bbebce 100644
--- a/scripts/create_logo.py
+++ b/scripts/create_logo.py
@@ -10,6 +10,7 @@ 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
 
 
@@ -21,21 +22,34 @@ def create_logo():
         None
     
     '''
-
-    filename = '../output/logo_magnetization.txt'
+    filename = '../output/mag_distr_logo.txt'
     b_0 = 1.0  # in T
     res = 10.0  # in nm
-    beta = pi/2
+    beta = pi/2  # in rad
     density = 10
+    dim = (128, 128)    
+    
+    x = range(dim[1])
+    y = range(dim[0])
+    xx, yy = np.meshgrid(x, y)    
+    bottom = (yy >= 0.25*dim[0])
+    left   = (yy <= 0.75/0.5 * dim[0]/dim[1] * xx)
+    right  = np.fliplr(left)
+    mag_shape = np.logical_and(np.logical_and(left, right), bottom)
     
-    mc.create_logo(128, res, beta, filename)
+    mc.create_hom_mag(dim, res, beta, mag_shape, filename)
     mag_data = dl.MagDataLLG(filename)
-    phase= pm.real_space_slab(mag_data, b_0)  
+    phase= pm.real_space(mag_data, 'slab', b_0)  
     holo = hi.holo_image(phase, mag_data.res, density)
-    hi.display_holo(holo, '')
+    hi.display_holo(holo, '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:
diff --git a/scripts/create_sample.py b/scripts/create_sample.py
new file mode 100644
index 0000000..89f0620
--- /dev/null
+++ b/scripts/create_sample.py
@@ -0,0 +1,59 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr 03 11:15:38 2013
+
+@author: Jan
+"""
+
+import pyramid.magcreator as mc
+import pdb, traceback, sys
+from numpy import pi
+
+
+def create_sample():
+    '''Calculate and display the phase map from a given magnetization.
+    Arguments:
+        None
+    Returns:
+        None
+    
+    '''
+    
+    # TODO: Input via GUI
+    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    
+    
+    if   key == 'slab':
+        mag_shape = mc.slab(dim, center, width)
+    elif key == 'disc':
+        mag_shape = mc.disc(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)
+    elif key == 'pixel':
+        mag_shape = mc.single_pixel(dim, pixel)
+    
+    mc.create_hom_mag(dim, res, beta, mag_shape, filename, plot_mag_distr)
+    
+    
+if __name__ == "__main__":
+    try:
+        create_sample()
+    except:
+        type, value, tb = sys.exc_info()
+        traceback.print_exc()
+        pdb.post_mortem(tb)
\ No newline at end of file
diff --git a/scripts/create_samples.py b/scripts/create_samples.py
deleted file mode 100644
index cd8a8cd..0000000
--- a/scripts/create_samples.py
+++ /dev/null
@@ -1,172 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Apr 03 11:15:38 2013
-
-@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
-from numpy import pi
-
-
-def phase_from_mag():
-    '''Calculate and display the phase map from a given magnetization.
-    Arguments:
-        None
-    Returns:
-        None
-    
-    '''
-    
-    '''INPUT'''
-    # TODO: Input via GUI
-    filename = '../output/output.txt'
-    b_0 = 10.0  # in T
-    v_0 = 0  # TODO: units?
-    v_acc = 30000  # in V
-    padding = 20
-    density = 100
-    
-    dim = (50, 50)  # in px (y,x)
-    res = 1.0  # in nm
-    beta = pi/4
-    
-    plot_mag_distr = True
-    
-    # Slab:
-    shape_fun = mc.slab
-    center = (24, 24)  # in px (y,x) index starts with 0!
-    width  = (25, 25)  # in px (y,x)
-    params = (center, width)
-#    # Disc:
-#    shape_fun = mc.disc
-#    center = (4, 4)  # in px (y,x)
-#    radius = 2.5
-#    params = (center, radius)
-#    # Filament:
-#    shape_fun = mc.filament
-#    pos = 5
-#    x_or_y = 'y'
-#    params = (pos, x_or_y)
-#    # Alternating Filaments:
-#    shape_fun = mc.alternating_filaments
-#    spacing = 5
-#    x_or_y = 'y'
-#    params = (spacing, x_or_y)
-#    # Single Pixel:
-#    shape_fun = mc.single_pixel
-#    pixel = (5, 5)
-#    params = pixel
-    
-#    '''CREATE LOGO'''    
-#    mc.create_logo(128, res, beta, filename, plot_mag_distr)
-#    mag_data = dl.MagDataLLG(filename)
-#    phase= pm.real_space_slab(mag_data, b_0)  
-#    holo = hi.holo_image(phase, mag_data.res, density)
-#    hi.display_holo(holo, '')
-    
-    '''CREATE MAGNETIC DISTRIBUTION'''
-    mc.create_hom_mag(dim, res, beta, shape_fun, params,
-                      filename, plot_mag_distr)
-    
-    '''LOAD MAGNETIC DISTRIBUTION'''
-    mag_data = dl.MagDataLLG(filename)
-    
-    
-    # TODO: get it to work:
-    phase_el = pm.phase_elec(mag_data, v_0=1, v_acc=200000)    
-    
-    
-    phi_cos_real_slab = np.zeros((2*dim[0]-1, 2*dim[1]-1))
-    pm.real_space_slab(mag_data, b_0, pixel_map=phi_cos_real_slab)
-    pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
-    phi_cos_real_disc = np.zeros((2*dim[0]-1, 2*dim[1]-1))
-    pm.real_space_disc(mag_data, b_0, pixel_map=phi_cos_real_disc)
-    pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
-    phi_cos_diff = phi_cos_real_slab - phi_cos_real_disc
-    pm.display_phase(phi_cos_diff, mag_data.res, 'Difference: One Pixel Slab - Disc')
-    
-    
-    '''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):
-    jacobi = np.zeros((2*dim[0]*dim[1], dim[0]*dim[1]))
-    tic = time.clock()
-    phase_real_slab = pm.real_space_slab(mag_data, b_0, jacobi=None)
-    toc = time.clock()
-    np.savetxt('../output/jacobi.npy', jacobi)
-    print 'Time for Real Space Approach (Slab): ' + str(toc - tic)
-    holo_real_slab = hi.holo_image(phase_real_slab, mag_data.res, density)
-    display_combined(phase_real_slab, mag_data.res, holo_real_slab, 
-                     'Real Space Approach (Slab)')
-    # numerical solution Real Space (Disc):
-    tic = time.clock()
-    phase_real_disc = pm.real_space_disc(mag_data, b_0)
-    toc = time.clock()
-    print 'Time for Real Space Approach (Disc): ' + str(toc - tic)
-    holo_real_disc = hi.holo_image(phase_real_disc, mag_data.res, density)
-    display_combined(phase_real_disc, mag_data.res, holo_real_disc,
-                     'Real Space Approach (Disc)')
-    
-    '''ANALYTIC SOLUTION'''
-    # analytic solution slab:
-    phase = an.phasemap_slab(dim, res, beta, center, width, b_0)
-    holo  = hi.holo_image(phase, res, density)
-    display_combined(phase, res, holo, 'Analytic Solution - Slab')
-#    # analytic solution disc:
-#    phase = an.phasemap_disc(dim, res, beta, center, radius, b_0)
-#    holo  = hi.holo_image(phase, res, density)
-#    display_combined(phase, res, holo, 'Analytic Solution - Disc')
-#    # analytic solution sphere:
-#    phase = an.phasemap_sphere(dim, res, beta, center, radius, b_0)
-#    holo  = hi.holo_image(phase, res, density)
-#    display_combined(phase, res, holo, 'Analytic Solution - Sphere')
-
-    
-    '''DIFFERENCES'''
-    diff_fft_to_ana       = phase_fft       - phase 
-    diff_real_slab_to_ana = phase_real_slab - phase
-    diff_real_disc_to_ana = phase_real_disc - phase
-    diff_slab_to_disc     = phase_real_disc - phase_real_slab
-    pm.display_phase(diff_fft_to_ana,       res, 'Difference: FFT - Analytic')
-    pm.display_phase(diff_real_slab_to_ana, res, 'Difference: Slab - Analytic')
-    pm.display_phase(diff_real_disc_to_ana, res, 'Difference: Disc - Analytic')
-    pm.display_phase(diff_slab_to_disc,     res, 'Difference: Disc - Slab')
-    
-    # TODO: Delete
-#    import pdb; pdb.set_trace()
- 
-    
-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)
-    
-    
-if __name__ == "__main__":
-    try:
-        phase_from_mag()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)
\ No newline at end of file
diff --git a/scripts/cython_example.py b/scripts/cython_example.py
new file mode 100644
index 0000000..6da48b0
--- /dev/null
+++ b/scripts/cython_example.py
@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri May 03 09:20:51 2013
+
+@author: Jan
+"""
+
+import cython
+import numpy as np
+
+
+n= 100
+A = np.array(range(n**2)).reshape((n,n))
+
+
+def naive_convolve(f, g):
+    # f is an image and is indexed by (v, w)
+    # g is a filter kernel and is indexed by (s, t),
+    #   it needs odd dimensions
+    # h is the output image and is indexed by (x, y),
+    #   it is not cropped
+    if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
+        raise ValueError("Only odd dimensions on filter supported")
+    # smid and tmid are number of pixels between the center pixel
+    # and the edge, ie for a 5x5 filter they will be 2.
+    #
+    # The output size is calculated by adding smid, tmid to each
+    # side of the dimensions of the input image.
+    vmax = f.shape[0]
+    wmax = f.shape[1]
+    smax = g.shape[0]
+    tmax = g.shape[1]
+    smid = smax // 2
+    tmid = tmax // 2
+    xmax = vmax + 2*smid
+    ymax = wmax + 2*tmid
+    # Allocate result image.
+    h = np.zeros([xmax, ymax], dtype=f.dtype)
+    # Do convolution
+    for x in range(xmax):
+        for y in range(ymax):
+            # Calculate pixel value for h at (x,y). Sum one component
+            # for each pixel (s, t) of the filter g.
+            s_from = max(smid - x, -smid)
+            s_to = min((xmax - x) - smid, smid + 1)
+            t_from = max(tmid - y, -tmid)
+            t_to = min((ymax - y) - tmid, tmid + 1)
+            value = 0
+            for s in range(s_from, s_to):
+                for t in range(t_from, t_to):
+                    v = x - smid + s
+                    w = y - tmid + t
+                    value += g[smid - s, tmid - t] * f[v, w]
+            h[x, y] = value
+    return h
\ No newline at end of file
diff --git a/scripts/get_jacobi.py b/scripts/get_jacobi.py
index d6b04c7..ad3091b 100644
--- a/scripts/get_jacobi.py
+++ b/scripts/get_jacobi.py
@@ -6,12 +6,9 @@ Created on Wed Apr 03 11:15:38 2013
 """
 
 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
 from numpy import pi
@@ -25,147 +22,38 @@ def phase_from_mag():
         None
     
     '''
-    
-    '''INPUT'''
     # TODO: Input via GUI
     filename = '../output/output.txt'
-    b_0 = 1.0  # in T
-    v_0 = 0  # TODO: units?
-    v_acc = 30000  # in V
-    padding = 20
-    density = 100
-    
-    dim = (50, 50)  # in px (y,x)
-    res = 1.0  # in nm
-    beta = pi/4
-    
-    plot_mag_distr = True
+    b_0 = 1.0  # in T    
+    dim = (3, 3)  # in px (y,x)
+    res = 10.0  # in nm
+    beta = pi/4    
     
-    # Slab:
-    shape_fun = mc.slab
-    center = (24, 24)  # in px (y,x) index starts with 0!
-    width  = (25, 25)  # in px (y,x)
-    params = (center, width)
-#    # Disc:
-#    shape_fun = mc.disc
-#    center = (4, 4)  # in px (y,x)
-#    radius = 2.5
-#    params = (center, radius)
-#    # Filament:
-#    shape_fun = mc.filament
-#    pos = 5
-#    x_or_y = 'y'
-#    params = (pos, x_or_y)
-#    # Alternating Filaments:
-#    shape_fun = mc.alternating_filaments
-#    spacing = 5
-#    x_or_y = 'y'
-#    params = (spacing, x_or_y)
-#    # Single Pixel:
-#    shape_fun = mc.single_pixel
-#    pixel = (5, 5)
-#    params = pixel
-    
-#    '''CREATE LOGO'''    
-#    mc.create_logo(128, res, beta, filename, plot_mag_distr)
-#    mag_data = dl.MagDataLLG(filename)
-#    phase= pm.real_space_slab(mag_data, b_0)  
-#    holo = hi.holo_image(phase, mag_data.res, density)
-#    hi.display_holo(holo, '')
+    center = (1, 1)  # in px (y,x) index starts with 0!
+    width  = (1, 1)  # in px (y,x)
+    mag_shape = mc.slab(dim, center, width)
     
     '''CREATE MAGNETIC DISTRIBUTION'''
-    mc.create_hom_mag(dim, res, beta, shape_fun, params,
-                      filename, plot_mag_distr)
+    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=1, v_acc=200000)    
-    
-    
-    phi_cos_real_slab = np.zeros((2*dim[0]-1, 2*dim[1]-1))
-    pm.real_space_slab(mag_data, b_0, pixel_map=phi_cos_real_slab)
-    pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
-    phi_cos_real_disc = np.zeros((2*dim[0]-1, 2*dim[1]-1))
-    pm.real_space_disc(mag_data, b_0, pixel_map=phi_cos_real_disc)
-    pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
-    phi_cos_diff = phi_cos_real_slab - phi_cos_real_disc
-    pm.display_phase(phi_cos_diff, mag_data.res, 'Difference: One Pixel Slab - Disc')
-    
-    
     '''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):
-    jacobi = np.zeros((2*dim[0]*dim[1], dim[0]*dim[1]))
+    jacobi = np.zeros((dim[0]*dim[1], 2*dim[0]*dim[1]))
     tic = time.clock()
-    phase_real_slab = pm.real_space_slab(mag_data, b_0, jacobi=None)
+    pm.real_space(mag_data, 'slab', b_0, jacobi=jacobi)
     toc = time.clock()
     np.savetxt('../output/jacobi.npy', jacobi)
-    print 'Time for Real Space Approach (Slab): ' + str(toc - tic)
-    holo_real_slab = hi.holo_image(phase_real_slab, mag_data.res, density)
-    display_combined(phase_real_slab, mag_data.res, holo_real_slab, 
-                     'Real Space Approach (Slab)')
-    # numerical solution Real Space (Disc):
-    tic = time.clock()
-    phase_real_disc = pm.real_space_disc(mag_data, b_0)
-    toc = time.clock()
-    print 'Time for Real Space Approach (Disc): ' + str(toc - tic)
-    holo_real_disc = hi.holo_image(phase_real_disc, mag_data.res, density)
-    display_combined(phase_real_disc, mag_data.res, holo_real_disc,
-                     'Real Space Approach (Disc)')
-    
-    '''ANALYTIC SOLUTION'''
-    # analytic solution slab:
-    phase = an.phasemap_slab(dim, res, beta, center, width, b_0)
-    holo  = hi.holo_image(phase, res, density)
-    display_combined(phase, res, holo, 'Analytic Solution - Slab')
-#    # analytic solution disc:
-#    phase = an.phasemap_disc(dim, res, beta, center, radius, b_0)
-#    holo  = hi.holo_image(phase, res, density)
-#    display_combined(phase, res, holo, 'Analytic Solution - Disc')
-#    # analytic solution sphere:
-#    phase = an.phasemap_sphere(dim, res, beta, center, radius, b_0)
-#    holo  = hi.holo_image(phase, res, density)
-#    display_combined(phase, res, holo, 'Analytic Solution - Sphere')
-
-    
-    '''DIFFERENCES'''
-    diff_fft_to_ana       = phase_fft       - phase 
-    diff_real_slab_to_ana = phase_real_slab - phase
-    diff_real_disc_to_ana = phase_real_disc - phase
-    diff_slab_to_disc     = phase_real_disc - phase_real_slab
-    pm.display_phase(diff_fft_to_ana,       res, 'Difference: FFT - Analytic')
-    pm.display_phase(diff_real_slab_to_ana, res, 'Difference: Slab - Analytic')
-    pm.display_phase(diff_real_disc_to_ana, res, 'Difference: Disc - Analytic')
-    pm.display_phase(diff_slab_to_disc,     res, 'Difference: Disc - Slab')
-    
-    # TODO: Delete
-#    import pdb; pdb.set_trace()
- 
-    
-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)
+    print 'Time for Real Space Approach with Jacobi-Matrix (Slab): ' + str(toc - tic)
     
-    phase_axis = fig.add_subplot(1,2,2, aspect='equal')
-    pm.display_phase(phase, res, 'Phasemap', phase_axis)
+    return jacobi
     
     
 if __name__ == "__main__":
     try:
-        phase_from_mag()
+        jacobi = phase_from_mag()
     except:
         type, value, tb = sys.exc_info()
         traceback.print_exc()
diff --git a/scripts/main.py b/scripts/main.py
deleted file mode 100644
index 6ccadb5..0000000
--- a/scripts/main.py
+++ /dev/null
@@ -1,172 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Apr 03 11:15:38 2013
-
-@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
-from numpy import pi
-
-
-def phase_from_mag():
-    '''Calculate and display the phase map from a given magnetization.
-    Arguments:
-        None
-    Returns:
-        None
-    
-    '''
-    
-    '''INPUT'''
-    # TODO: Input via GUI
-    filename = '../output/output.txt'
-    b_0 = 1.0  # in T
-    v_0 = 0  # TODO: units?
-    v_acc = 30000  # in V
-    padding = 20
-    density = 10
-    
-    dim = (50, 50)  # in px (y,x)
-    res = 10.0  # in nm
-    beta = pi/4
-    
-    plot_mag_distr = True
-    
-    # Slab:
-    shape_fun = mc.slab
-    center = (24, 24)  # in px (y,x) index starts with 0!
-    width  = (25, 25)  # in px (y,x)
-    params = (center, width)
-#    # Disc:
-#    shape_fun = mc.disc
-#    center = (4, 4)  # in px (y,x)
-#    radius = 2.5
-#    params = (center, radius)
-#    # Filament:
-#    shape_fun = mc.filament
-#    pos = 5
-#    x_or_y = 'y'
-#    params = (pos, x_or_y)
-#    # Alternating Filaments:
-#    shape_fun = mc.alternating_filaments
-#    spacing = 5
-#    x_or_y = 'y'
-#    params = (spacing, x_or_y)
-#    # Single Pixel:
-#    shape_fun = mc.single_pixel
-#    pixel = (5, 5)
-#    params = pixel
-    
-#    '''CREATE LOGO'''    
-#    mc.create_logo(128, res, beta, filename, plot_mag_distr)
-#    mag_data = dl.MagDataLLG(filename)
-#    phase= pm.real_space_slab(mag_data, b_0)  
-#    holo = hi.holo_image(phase, mag_data.res, density)
-#    hi.display_holo(holo, '')
-    
-    '''CREATE MAGNETIC DISTRIBUTION'''
-    mc.create_hom_mag(dim, res, beta, shape_fun, params,
-                      filename, plot_mag_distr)
-    
-    '''LOAD MAGNETIC DISTRIBUTION'''
-    mag_data = dl.MagDataLLG(filename)
-    
-    
-    # TODO: get it to work:
-    phase_el = pm.phase_elec(mag_data, v_0=1, v_acc=200000)    
-    
-    
-    phi_cos_real_slab = np.zeros((2*dim[0]-1, 2*dim[1]-1))
-    pm.real_space_slab(mag_data, b_0, pixel_map=phi_cos_real_slab)
-    pm.display_phase(phi_cos_real_slab, res, 'Phase of one Pixel-Slab (Cos - Part)')
-    phi_cos_real_disc = np.zeros((2*dim[0]-1, 2*dim[1]-1))
-    pm.real_space_disc(mag_data, b_0, pixel_map=phi_cos_real_disc)
-    pm.display_phase(phi_cos_real_disc, res, 'Phase of one Pixel-Disc (Cos - Part)')
-    phi_cos_diff = phi_cos_real_slab - phi_cos_real_disc
-    pm.display_phase(phi_cos_diff, mag_data.res, 'Difference: One Pixel Slab - Disc')
-    
-    
-    '''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):
-    jacobi = np.zeros((2*dim[0]*dim[1], dim[0]*dim[1]))
-    tic = time.clock()
-    phase_real_slab = pm.real_space_slab(mag_data, b_0, jacobi=None)
-    toc = time.clock()
-    np.savetxt('../output/jacobi.npy', jacobi)
-    print 'Time for Real Space Approach (Slab): ' + str(toc - tic)
-    holo_real_slab = hi.holo_image(phase_real_slab, mag_data.res, density)
-    display_combined(phase_real_slab, mag_data.res, holo_real_slab, 
-                     'Real Space Approach (Slab)')
-    # numerical solution Real Space (Disc):
-    tic = time.clock()
-    phase_real_disc = pm.real_space_disc(mag_data, b_0)
-    toc = time.clock()
-    print 'Time for Real Space Approach (Disc): ' + str(toc - tic)
-    holo_real_disc = hi.holo_image(phase_real_disc, mag_data.res, density)
-    display_combined(phase_real_disc, mag_data.res, holo_real_disc,
-                     'Real Space Approach (Disc)')
-    
-    '''ANALYTIC SOLUTION'''
-    # analytic solution slab:
-    phase = an.phasemap_slab(dim, res, beta, center, width, b_0)
-    holo  = hi.holo_image(phase, res, density)
-    display_combined(phase, res, holo, 'Analytic Solution - Slab')
-#    # analytic solution disc:
-#    phase = an.phasemap_disc(dim, res, beta, center, radius, b_0)
-#    holo  = hi.holo_image(phase, res, density)
-#    display_combined(phase, res, holo, 'Analytic Solution - Disc')
-#    # analytic solution sphere:
-#    phase = an.phasemap_sphere(dim, res, beta, center, radius, b_0)
-#    holo  = hi.holo_image(phase, res, density)
-#    display_combined(phase, res, holo, 'Analytic Solution - Sphere')
-
-    
-    '''DIFFERENCES'''
-    diff_fft_to_ana       = phase_fft       - phase 
-    diff_real_slab_to_ana = phase_real_slab - phase
-    diff_real_disc_to_ana = phase_real_disc - phase
-    diff_slab_to_disc     = phase_real_disc - phase_real_slab
-    pm.display_phase(diff_fft_to_ana,       res, 'Difference: FFT - Analytic')
-    pm.display_phase(diff_real_slab_to_ana, res, 'Difference: Slab - Analytic')
-    pm.display_phase(diff_real_disc_to_ana, res, 'Difference: Disc - Analytic')
-    pm.display_phase(diff_slab_to_disc,     res, 'Difference: Disc - Slab')
-    
-    # TODO: Delete
-#    import pdb; pdb.set_trace()
- 
-    
-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)
-    
-    
-if __name__ == "__main__":
-    try:
-        phase_from_mag()
-    except:
-        type, value, tb = sys.exc_info()
-        traceback.print_exc()
-        pdb.post_mortem(tb)
\ No newline at end of file
diff --git a/src/pyramex/.sconsign.dblite b/src/pyramex/.sconsign.dblite
new file mode 100644
index 0000000000000000000000000000000000000000..0e8ccb44f58bbb31c7780f36f464f52b6d17a122
GIT binary patch
literal 2379
zcmdUxe{<V37{^y=*VSFSb?v%#uo0j@U^KBUOSa4S(S@`dOFOQq5B{hxXUTG8u1)g9
z4y@eyz3?XB^&#0w>Wtm>Uv`)C>Eu}YKF|02k)6R@)mvwKXM@XB^n0d#mat(wNR!il
z=5E}`V9wK6ey7&_wF&Pr)*oDUhP<oE3+wlibU0#Zf2hIxRn~6EiCTV*yFG0(X!g?K
zr3Uk?)1jmIJdJyuj@GBYah7VZaGajN%?uW+-a;?#MSZv>rb`|7LtZ;*YX|W;SK}A*
zz6MJoe@M$<x#}%<e&@p^?)N%z2)AF(c}r=0!ILz&7{E%)K-jW%0~t<B^X6F+cVX45
zAO}%|kU^On1VlG2?$~-jImM<CI6(w=y!j;llf&H>aS1Y8+K74!G*`Fwc6HBWVPPTf
zULkKigN>?pizfZiknsfW``!|bF9zo^i&NP2<*kqpcrQ%gfiINu<;g_cO<>E1ZEukV
zELEZ(9(ya(hvh^9ANcU0d}6%1nt9~I;}IRvq~@(jKQH25_YY+5zP>z@%1{3uySdn9
z>>agIM>Zz*8&<klk4zV1n?-u$TGR|89O*WucH|Obm8>-6-;iy+b7s!fPXAiSQLK6e
zQ#o`0Ll#Tu-InMP+jcD3bBW#zr7JS)1z)Qz%f@9x?p#aP-GZ*Q4DOZbTA!h7qfFO*
zU#J~t=-Tw*!4+Lw|DkK!hljItDc%%a@)R*@eJNAS7fK#w@VF|}p(Vx#(?=Ops$$Xs
zgCL`&<r0bm=wq*fbz+8rJX%OZ42}ZbiQFi}5l0p`BVviwmTAh?sn!V96EQCPlMFts
zigh1$avcml3&e&GpL^Slz4}qyJNzXlzAlNcPeA`75Jd~brZ2YS=8<I*Ofd}!vKV6y
z!$dc@W8jcFEVP6yzFI2h%M89M%DDnil_O*!=o<WFAtmerZ(bc}@Wgw%(|+D)RCfE7
zozwl3!{%vamsQl<L2Li8osTBNPNUIybtHds_LX<iXtss);_H&Y-3)4P>$ehU5ZlEL
zL)?il)saq|fa?~;oKQQA2$ljJL$`I+=^*#_ru27k`C6?%cDSB_KI^cNfmzLaFm_qT
zwrfk5<&L3&!oajq7&;j1Om5>c$CN|eB7w_wKypvy%1pfBt(5^*-bzES3<Dyzj+G;}
zs@j<*kbzrLxR=2<vkISP@a>F3Ng|O|DuqEP-8MOoEJAe5EbK;@1U93FiA@v`_^#ZQ
znf!g9!4G-<_PueR;8~u)hWA(|@aOS2M7>>nL)1&34nLM+2u&44%8kfkw&4=X32|si
zNh~x7lbV>TQwj&F@WYbR&og*YIDPe6pGNVzm-<h#s_E-I5Kp|-Lfs@+@KUw7Rod<-
VgX3wiCKl((7bO$(*o<nU+P^hG&hh{N

literal 0
HcmV?d00001

diff --git a/src/pyramex/SConstruct b/src/pyramex/SConstruct
new file mode 100644
index 0000000..dbcb3e9
--- /dev/null
+++ b/src/pyramex/SConstruct
@@ -0,0 +1,16 @@
+import distutils.sysconfig
+import os
+import sys
+
+PYTHON_LIBPATH = distutils.sysconfig.get_python_lib()
+PYTHON_INCPATH = distutils.sysconfig.get_python_inc()
+PYTHON_LIBRARY = "python" + sys.version[:3]
+
+cythonBuilder = Builder(action ='cython $SOURCE',
+                        suffix = '.c',
+                        src_suffix = '.pyx')
+
+env = Environment(ENV      =   os.environ,
+                  BUILDERS = {'Py2C':cythonBuilder})
+
+env.Py2C('c1.c', 'c1.pyx')
\ No newline at end of file
diff --git a/src/pyramex/SConstructBACKUP b/src/pyramex/SConstructBACKUP
new file mode 100644
index 0000000..58a920c
--- /dev/null
+++ b/src/pyramex/SConstructBACKUP
@@ -0,0 +1,70 @@
+import distutils.sysconfig
+import os
+import sys
+
+PYTHON_LIBPATH = distutils.sysconfig.get_python_lib()
+PYTHON_INCPATH = distutils.sysconfig.get_python_inc()
+PYTHON_LIBRARY = "python" + sys.version[:3]
+
+print 'PYTHON_LIBPATH: ' + PYTHON_LIBPATH
+print 'PYTHON_INCPATH: ' + PYTHON_INCPATH
+print 'PYTHON_LIBRARY: ' + PYTHON_LIBRARY
+
+env = Environment(ENV = os.environ)
+
+env.AppendUnique(LIBPATH=[PYTHON_LIBPATH], CPPPATH=[PYTHON_INCPATH])
+
+if ARGUMENTS.get('VERBOSE') != '1':
+    env['CCCOMSTR']   = "Compiling $TARGET"
+    env['LINKCOMSTR'] = "Linking $TARGET"
+
+Progress('$TARGET\r', overwrite=True)
+
+env.Decider('MD5-timestamp')
+
+
+
+
+
+
+
+
+
+
+
+# This builder takes the cython source code and generates c code.
+cythonBuilder = Builder(action ='cython $SOURCE',
+                        suffix = '.c',
+                        src_suffix = '.pyx')
+
+# Setup the compilation environment.
+env = Environment(PYEXT_USE_DISTUTILS=True,
+            ENV = os.environ,
+            tools=['mingw'],
+            BUILDERS={'Py2C':cythonBuilder},
+            CCFLAGS='-LC:\Python27\libs -IC:\Python27\include',
+            LINKFLAGS='-lpython2.7 -shared -o $SOURCE')
+
+# Method that calls routines neccessary to create shared library.
+def cythonPseudoBuilder(env,source):
+    cCode = env.Py2C(source)
+    #env.SharedLibrary(source,cCode,LIBPREFIX='')
+
+env.AddMethod(cythonPseudoBuilder,'Cython')
+
+env.Cython('c1')
+
+#env.SharedLibrary('c1.pyd', ['c1.c'])
+
+
+c1_module = env.LoadableModule('c1.pyd', ['c1.c'])
+
+
+
+
+#import setup
+
+#for ext in setup.get_extensions():
+#	print e.sources[:]
+
+#env.Command('hello.pyd', '', 'python setup.py build_ext --compiler=mingw32 --inplace')
\ No newline at end of file
diff --git a/src/pyramex/c1.pyx b/src/pyramex/c1.pyx
new file mode 100644
index 0000000..bfa8239
--- /dev/null
+++ b/src/pyramex/c1.pyx
@@ -0,0 +1,13 @@
+import math
+
+def great_circle(float lon1,float lat1,float lon2,float lat2):
+    cdef float radius = 3956.0 
+    cdef float pi = 3.14159265
+    cdef float x = pi/180.0
+    cdef float a,b,theta,c
+
+    a = (90.0-lat1)*(x)
+    b = (90.0-lat2)*(x)
+    theta = (lon2-lon1)*(x)
+    c = math.acos((math.cos(a)*math.cos(b)) + (math.sin(a)*math.sin(b)*math.cos(theta)))
+    return radius*c
\ No newline at end of file
diff --git a/src/pyramex/c2.pyx b/src/pyramex/c2.pyx
new file mode 100644
index 0000000..c01cfd0
--- /dev/null
+++ b/src/pyramex/c2.pyx
@@ -0,0 +1,16 @@
+cdef extern from "math.h":
+    float cosf(float theta)
+    float sinf(float theta)
+    float acosf(float theta)
+
+def great_circle(float lon1,float lat1,float lon2,float lat2):
+    cdef float radius = 3956.0 
+    cdef float pi = 3.14159265
+    cdef float x = pi/180.0
+    cdef float a,b,theta,c
+
+    a = (90.0-lat1)*(x)
+    b = (90.0-lat2)*(x)
+    theta = (lon2-lon1)*(x)
+    c = acosf((cosf(a)*cosf(b)) + (sinf(a)*sinf(b)*cosf(theta)))
+    return radius*c
\ No newline at end of file
diff --git a/src/pyramex/c3.pyx b/src/pyramex/c3.pyx
new file mode 100644
index 0000000..1ffbe36
--- /dev/null
+++ b/src/pyramex/c3.pyx
@@ -0,0 +1,23 @@
+cdef extern from "math.h":
+    float cosf(float theta)
+    float sinf(float theta)
+    float acosf(float theta)
+
+cdef float _great_circle(float lon1,float lat1,float lon2,float lat2):
+    cdef float radius = 3956.0 
+    cdef float pi = 3.14159265
+    cdef float x = pi/180.0
+    cdef float a,b,theta,c
+
+    a = (90.0-lat1)*(x)
+    b = (90.0-lat2)*(x)
+    theta = (lon2-lon1)*(x)
+    c = acosf((cosf(a)*cosf(b)) + (sinf(a)*sinf(b)*cosf(theta)))
+    return radius*c
+
+def great_circle(float lon1,float lat1,float lon2,float lat2,int num):
+    cdef int i
+    cdef float x
+    for i from 0 <= i < num:
+        x = _great_circle(lon1,lat1,lon2,lat2)
+    return x
\ No newline at end of file
diff --git a/src/pyramex/hello.pyx b/src/pyramex/hello.pyx
new file mode 100644
index 0000000..3ca5a38
--- /dev/null
+++ b/src/pyramex/hello.pyx
@@ -0,0 +1,2 @@
+def say_hello_to(name):
+    print 'Hello %s!' % name
\ No newline at end of file
diff --git a/src/pyramex/setup.py b/src/pyramex/setup.py
new file mode 100644
index 0000000..b802b79
--- /dev/null
+++ b/src/pyramex/setup.py
@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri May 03 10:27:04 2013
+
+@author: Jan
+"""
+
+#from distutils.core import setup, Extension
+#
+#module1 = Extension('c1',
+#                    sources = ['c1.pyx'])
+#
+#setup (name = 'pyramex',
+#       version = '1.0',
+#       description = 'This is a demo package',
+#       ext_modules = [module1])
+
+
+
+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Distutils import build_ext
+
+def get_extensions():
+	return [
+		Extension('hello', ['hello.pyx']),
+		Extension('c1',    ['c1.pyx']),
+		Extension('c2',    ['c2.pyx']),
+		Extension('c3',    ['c3.pyx'])
+		]
+
+setup(
+  name = 'pyramex',
+  cmdclass = {'build_ext': build_ext},
+  ext_modules = get_extensions()
+)
\ No newline at end of file
diff --git a/src/pyramid/dataloader.py b/src/pyramid/dataloader.py
index 248c2b4..5900c56 100644
--- a/src/pyramid/dataloader.py
+++ b/src/pyramid/dataloader.py
@@ -18,21 +18,21 @@ class MagDataLLG:
         
         '''
         scale = 1.0E-9 / 1.0E-2  #from cm to nm
-        data = np.genfromtxt(filename, skip_header=2)
+        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_len, y_len, z_len = [data[-1, i]/scale+res/2 for i in range(3)]
-        x_mag, y_mag, z_mag = [data[:, i].reshape(z_dim, y_dim, x_dim).mean(0) 
-                               *z_len for i in range(3,6)]
+        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!
         self.filename = filename
         self.res = res
-        self.dim = (x_dim, y_dim, z_dim)
-        self.length = (x_len, y_len, z_len)
-        self.magnitude = (x_mag, y_mag, z_mag)
+        self.dim = (z_dim, y_dim, x_dim)
+        self.length = (z_len, y_len, x_len)
+        self.magnitude = (z_mag, y_mag, x_mag)
     
     def __str__(self):
         '''Return the filename of the loaded file.
diff --git a/src/pyramid/magcreator.py b/src/pyramid/magcreator.py
index ac00677..a497325 100644
--- a/src/pyramid/magcreator.py
+++ b/src/pyramid/magcreator.py
@@ -1,13 +1,11 @@
 # -*- coding: utf-8 -*-
 """Create simple LLG Files which describe magnetization in 2D (z-Dim=1)."""
 
-
 import numpy as np
 import matplotlib.pyplot as plt
-from numpy import pi
 
 
-def slab(dim, params):
+def slab(dim, center, width):
     '''Get the magnetic shape of a slab.
     Arguments:
         dim    - the dimensions of the grid, shape(y, x)
@@ -17,14 +15,13 @@ def slab(dim, params):
         the magnetic shape as a 2D-boolean-array.
         
     '''
-    center, width = params
-    shape_mag = np.array([[abs(x - center[1]) <= width[1] / 2
+    mag_shape = np.array([[abs(x - center[1]) <= width[1] / 2
                        and abs(y - center[0]) <= width[0] / 2
                        for x in range(dim[1])] for y in range(dim[0])])
-    return shape_mag
+    return mag_shape
 
 
-def disc(dim, params):
+def disc(dim, center, radius):
     '''Get the magnetic shape of a disc.
     Arguments:
         dim    - the dimensions of the grid, shape(y, x)
@@ -34,13 +31,12 @@ def disc(dim, params):
         the magnetic shape as a 2D-boolean-array.
         
     '''
-    center, radius = params
-    shape_mag = np.array([[(np.hypot(x-center[1], y-center[0]) <= radius)
+    mag_shape = np.array([[(np.hypot(x-center[1], y-center[0]) <= radius)
                            for x in range(dim[1])] for y in range(dim[0])])    
-    return shape_mag
+    return mag_shape
 
     
-def filament(dim, params):
+def 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)
@@ -51,16 +47,15 @@ def filament(dim, params):
         the magnetic shape as a 2D-boolean-array.
         
     '''
-    pos, x_or_y = params
-    shape_mag = np.zeros(dim)
+    mag_shape = np.zeros(dim)
     if x_or_y == 'y':
-        shape_mag[:, pos] = 1
+        mag_shape[:, pos] = 1
     else:
-        shape_mag[pos, :] = 1
-    return shape_mag
+        mag_shape[pos, :] = 1
+    return mag_shape
 
 
-def alternating_filaments(dim, params):
+def 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)
@@ -71,19 +66,18 @@ def alternating_filaments(dim, params):
         the magnetic shape as a 2D-boolean-array.
         
     '''
-    spacing, x_or_y = params
-    shape_mag = np.zeros(dim)
+    mag_shape = np.zeros(dim)
     if x_or_y == 'y':
         # TODO: List comprehension:
         for i in range(0, dim[1], spacing):
-            shape_mag[:, i] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1
+            mag_shape[:, i] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1
     else:
         for i in range(0, dim[0], spacing):
-            shape_mag[i, :] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1            
-    return shape_mag
+            mag_shape[i, :] = 1 - 2 * (int(i / spacing) % 2)  # 1 or -1            
+    return mag_shape
 
     
-def single_pixel(dim, params):
+def single_pixel(dim, pixel):
     '''Get the magnetic shape of a single magnetized pixel.
     Arguments:
         dim   - the dimensions of the grid, shape(y, x)
@@ -92,14 +86,12 @@ def single_pixel(dim, params):
         the magnetic shape as a 2D-boolean-array.
         
     '''
-    pixel = params
-    shape_mag = np.zeros(dim)
-    shape_mag[pixel[0], pixel[1]] = 1
-    return shape_mag
+    mag_shape = np.zeros(dim)
+    mag_shape[pixel[0], pixel[1]] = 1
+    return mag_shape
 
 
-def create_hom_mag(dim, res, beta, shape_fun, params,
-                   filename='output.txt', plot_mag_distr=False):
+def create_hom_mag(dim, res, beta, mag_shape, filename='output.txt', plot_mag_distr=False):
     '''Create homog. magnetization data, saved in a file with LLG-format.
     Arguments:
         dim       - the dimensions of the grid, shape(y, x)
@@ -107,23 +99,20 @@ def create_hom_mag(dim, res, beta, shape_fun, params,
                     (real space distance between two points)
         beta      - the angle of the magnetization
         filename  - the name of the file in which to save the data
-        shape_fun - the function to determine the shape of the magnetization,
-                    several are defined in this module
-        params   - a tuple of parameters used in shape_fun
+        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 shape as a 2D-boolean-array.
+        the the magnetic distribution as a 2D-boolean-array.
         
     '''
     res *= 1.0E-9 / 1.0E-2  # from nm to cm     
     
     x = np.linspace(res / 2, dim[1] * res - res / 2, num=dim[1])
     y = np.linspace(res / 2, dim[0] * res - res / 2, num=dim[0])
-    xx, yy = np.meshgrid(x, y)
-    
-    shape_mag = shape_fun(dim, params)                        
+    xx, yy = np.meshgrid(x, y)                       
                         
-    x_mag = np.array(np.ones(dim)) * np.cos(beta) * shape_mag
-    y_mag = np.array(np.ones(dim)) * np.sin(beta) * shape_mag
+    x_mag = np.array(np.ones(dim)) * np.cos(beta) * mag_shape
+    y_mag = np.array(np.ones(dim)) * np.sin(beta) * mag_shape
     z_mag = np.array(np.zeros(dim))
     
     if (plot_mag_distr):
@@ -138,49 +127,7 @@ def create_hom_mag(dim, res, beta, shape_fun, params,
     x_mag   = np.reshape(x_mag,(-1))
     y_mag   = np.reshape(y_mag,(-1))
     z_mag   = np.array(np.zeros(dim[0] * dim[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[1], dim[0], 1))
-        mag_file.writelines('\n'.join('   '.join('{:7.6e}'.format(cell) 
-                                      for cell in row) for row in data) )
-                                                       
-                                      
-def create_logo(edge, res, beta = pi/2, filename='logo.txt', plot_mag_distr=False):
-    
-    # TODO: rewrite so that every possible shape_mag can be given as argument
-    # TODO: write a script for creating, displaying and saving the logo
-    
-    dim = (edge, edge) 
-    res *= 1.0E-9 / 1.0E-2  # from nm to cm     
-    
-    x = np.linspace(res / 2, dim[1] * res - res / 2, num=dim[1])
-    y = np.linspace(res / 2, dim[0] * res - res / 2, num=dim[0])
-    xx, yy = np.meshgrid(x, y)
-    
-    bottom = (yy >= 0.25*edge*res)
-    left   = (yy <=            0.75/0.5 *  xx)
-    right  = np.fliplr(left)#(yy <= (edge-1)*res - 0.75/0.5 * (xx - 0.5*(edge-1)*res))
-    shape_mag = np.logical_and(np.logical_and(left, right), bottom)
     
-    x_mag = np.array(np.ones(dim)) * np.cos(beta) * shape_mag
-    y_mag = np.array(np.ones(dim)) * np.sin(beta) * shape_mag
-    z_mag = np.array(np.zeros(dim))
-    
-    if (plot_mag_distr):
-        fig = plt.figure()
-        fig.add_subplot(111, aspect='equal')
-        plt.quiver(x_mag, y_mag, pivot='middle', angles='xy', scale_units='xy', 
-                   scale=1, headwidth=6, headlength=7)    
-    
-    xx = np.reshape(xx,(-1))
-    yy = np.reshape(yy,(-1))
-    zz = np.array(np.ones(dim[0] * dim[1]) * res / 2)
-    x_mag   = np.reshape(x_mag,(-1))
-    y_mag   = np.reshape(y_mag,(-1))
-    z_mag   = np.array(np.zeros(dim[0] * dim[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', ''))
diff --git a/src/pyramid/phasemap.py b/src/pyramid/phasemap.py
index 96cdccb..6096799 100644
--- a/src/pyramid/phasemap.py
+++ b/src/pyramid/phasemap.py
@@ -14,7 +14,7 @@ Q_E   = 1.602E-19  # TODO: unit
 C     = 2.998E8    # TODO: unit
 
 
-def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
+def fourier_space(mag_data, b_0=1, padding=0):
     '''Calculate phasemap from magnetization data (fourier approach).
     Arguments:
         mag_data - MagDataLLG object (from magneticimaging.dataloader) storing
@@ -31,8 +31,8 @@ def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
     
     '''    
     res  = mag_data.res
-    x_dim, y_dim, z_dim = mag_data.dim
-    x_mag, y_mag, z_mag = mag_data.magnitude
+    z_dim, y_dim, x_dim = mag_data.dim
+    z_mag, y_mag, x_mag = mag_data.magnitude  
     
     # TODO: interpolation (redefine xDim,yDim,zDim) thickness ramp
     
@@ -40,6 +40,7 @@ def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
     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
     
@@ -50,7 +51,7 @@ def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
     y = np.linspace( -nyq/2, nyq/2, x_mag_fft.shape[0]+1)[:-1]
     xx, yy = np.meshgrid(x, y)
                          
-    phase_fft = (1j * b_0) / (2 * PHI_0) * ((x_mag_fft * yy - y_mag_fft * xx) 
+    phase_fft = (1j * res * b_0) / (2 * PHI_0) * ((x_mag_fft * yy - y_mag_fft * xx) 
                                            / (xx ** 2 + yy ** 2 + 1e-18))
     phase_big = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
     
@@ -59,56 +60,45 @@ def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
     # TODO: Electrostatic Component
 
     return phase
-    
-    
-def real_space_slab(mag_data, b_0=1, v_0=0, v_acc=30000, 
-                    pixel_map=None, jacobi=None):
+      
+      
+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):
     '''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
+        mag_data - MagDataLLG object (from magneticimaging.dataloader) storing the filename, 
+                   coordinates and magnetization in x, y and z
     Returns:
         the phasemap as a 2 dimensional array
         
-    '''
-    # TODO: combine with disc and sphere approach, pixel field as function
-    
-    # TODO: Expand docstring!
+    '''    
+    # TODO: Expand docstring!  
 
-    # TODO: Delete
-#    import pdb; pdb.set_trace()    
-    
     res  = mag_data.res
-    x_dim, y_dim, z_dim = mag_data.dim
-    x_mag, y_mag, z_mag = mag_data.magnitude    
+    z_dim, y_dim, x_dim = mag_data.dim
+    z_mag, y_mag, x_mag = mag_data.magnitude
+    
+    # TODO: proper projection algorithm       
+    x_mag, y_mag, z_mag = x_mag.mean(0), y_mag.mean(0), z_mag.mean(0)
     
     beta = np.arctan2(y_mag, x_mag)
 
-    mag = np.hypot(x_mag, y_mag)    
-     
-    coeff = - b_0 * res**2 / ( 2 * PHI_0 ) / res # TODO: why / res
-
-    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
-     
-    def phi_pixel(n, m):  # TODO: rename
-        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) )
+    mag = np.hypot(x_mag, y_mag) 
                 
-    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])
-                                          
-    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])
-    
     '''CREATE COORDINATE GRIDS'''
     x = np.linspace(0,(x_dim-1),num=x_dim)
     y = np.linspace(0,(y_dim-1),num=y_dim)
@@ -118,95 +108,38 @@ def real_space_slab(mag_data, b_0=1, v_0=0, v_acc=30000,
     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)
     
-    phi_cos = phi_pixel(xx_big, yy_big)
-    phi_sin = phi_pixel(yy_big, xx_big)
-    
-    if pixel_map is not None:
-        pixel_map[:] = phi_cos
+    phi_cos = phi_pixel(method, xx_big, yy_big, res, b_0)
+    phi_sin = phi_pixel(method, yy_big, xx_big, 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])
+                                          
+    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])   
     
     '''CALCULATE THE PHASE'''
     phase = np.zeros((y_dim, x_dim))
     
     # TODO: only iterate over pixels that have a magn. > threshold (first >0)
-    
-    
-#    import pdb; pdb.set_trace()
     for j in range(y_dim):
         for i in range(x_dim):
-            if (mag[j, i] != 0): # TODO: same result with or without?
+            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+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)         
     
     return phase
-    
-    
-def real_space_disc(mag_data, b_0=1, v_0=0, v_acc=30000,
-                    pixel_map=None, 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
-    Returns:
-        the phasemap as a 2 dimensional array
-        
-    '''
-    # TODO: Expand docstring!
 
-    # TODO: Delete
-#    import pdb; pdb.set_trace()    
     
-    res  = mag_data.res
-    x_dim, y_dim, z_dim = mag_data.dim
-    x_mag, y_mag, z_mag = mag_data.magnitude    
-    
-    beta = np.arctan2(y_mag, x_mag)
-
-    mag = np.hypot(x_mag, y_mag)
-    
-    coeff = - b_0 * res**2 / ( 2 * PHI_0 ) / res # TODO: why / res
 
-    def phi_pixel(n, m):
-        in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
-        result = coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
-        return result
-#        r = np.hypot(n, m)
-#        r[y_dim-1,x_dim-1] = 1E-18  # Prevent div zero error at disc center
-#        return coeff / r**2 * m * (r != 0) # (r > R) = 0 for disc center
-
-    def phi_mag(i, j):
-        return mag[j,i]*(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])
-    
-    '''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)
-     
-    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)
-    
-    phi_cos = phi_pixel(xx_big, yy_big)
-    phi_sin = phi_pixel(yy_big, xx_big)
-    
-    if pixel_map is not None:
-        pixel_map[:] = phi_cos
-    
-    '''CALCULATE THE PHASE'''
-    phase = np.zeros((y_dim, x_dim))
-    
-    # TODO: only iterate over pixels that have a magn. > threshold (first >0)
-    for j in range(y_dim):
-        for i in range(x_dim):
-#            if (mag[j, i] != 0): # TODO: same result with or without?
-                phase += phi_mag(i, j)
-    
-    return phase
     
 
 def phase_elec(mag_data, b_0=1, v_0=0, v_acc=30000):
@@ -216,8 +149,8 @@ def phase_elec(mag_data, b_0=1, v_0=0, v_acc=30000):
 #    import pdb; pdb.set_trace()
     
     res  = mag_data.res
-    x_dim, y_dim, z_dim = mag_data.dim
-    x_mag, y_mag, z_mag = mag_data.magnitude 
+    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)    
     
diff --git a/test/run_tests.py b/test/run_tests.py
new file mode 100644
index 0000000..6b3b049
--- /dev/null
+++ b/test/run_tests.py
@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+"""
+Unittests for pyramid.
+"""
+
+import unittest
+import sys
+from test_dataloader import *
+from test_compliance import *
+
+
+def run():
+    
+    suite  = unittest.TestSuite()
+    loader = unittest.TestLoader()
+    runner = unittest.TextTestRunner(verbosity=2)
+    
+    suite.addTest(loader.loadTestsFromTestCase(TestCaseDataloader))
+    suite.addTest(loader.loadTestsFromTestCase(TestCaseCompliance))
+    runner.run(suite)
+    
+    #TODO: why that?
+#    result = runner.run(suite)
+#    if result.wasSuccessful():
+#        sys.exit(0)
+#    else:
+#        sys.exit(1)
+
+
+if __name__ == '__main__':
+    run()
+    
+    
+    
+    
+
+#if __name__ == '__main__':
+#    suite = unittest.TestLoader().loadTestsFromTestCase(TestSuite)
+#    unittest.TextTestRunner(verbosity=2).run(suite)
+
+
+
+''' GLORIPY VERSION '''
+##! /usr/bin/python
+#"""
+#Unittests for pyjurassic.
+#"""
+#
+#import unittest
+#import os
+#import sys
+#import tests
+#
+#
+#def main():
+#    all_tests = unittest.TestSuite()
+#    tl = unittest.defaultTestLoader
+#    if os.getenv("BOOST_TESTS_TO_RUN") is not None:
+#        return
+#    elif os.getenv("PYTHON_TESTS_TO_RUN") is None:
+#        all_tests.addTest(tl.loadTestsFromName("tests"))
+#    else:
+#        all_tests.addTest(tl.loadTestsFromName(os.getenv("PYTHON_TESTS_TO_RUN")))
+#
+#    runner = unittest.TextTestRunner(verbosity=2)
+#    res = runner.run(all_tests)
+#    if res.wasSuccessful():
+#        sys.exit(0)
+#    else:
+#        sys.exit(1)
+#
+#
+#if __name__ == '__main__':
+#    main()
\ No newline at end of file
diff --git a/test/test_compliance.py b/test/test_compliance.py
new file mode 100644
index 0000000..9c7d370
--- /dev/null
+++ b/test/test_compliance.py
@@ -0,0 +1,82 @@
+import os
+import unittest
+
+
+class TestCompliance(unittest.TestCase):
+    """
+    Class for checking compliance of pyjurassic.
+    """
+
+    def setUp(self):
+        # self.getPaths()
+        # self.data_dir = "test_jrp_data"
+        # self.src_ctl = "ctl.py"
+        # remove E201 to get a working version of str(- b)
+        # remove E203 to get a working version of a[:, :]
+        # remove E501 as 80 chars are not enough
+        self.ignores = ["E201", "E203", "E501"]
+        self.options = ["dummy_name", "--ignore", ",".join(self.ignores), "--repeat"]
+        self.options += ["--exclude", "collections_python27.py"]
+        try:
+            import pep8
+            self.pep8 = pep8
+        except ImportError:
+            self.pep8 = None
+            print("\nWARNING: You do not have package pep8!")
+
+    def countErrors(self):
+        """
+        Counts the relevant errors from the pep8 counter structure.
+        """
+        return sum([self.pep8.options.counters[key]
+                    for key in self.pep8.options.messages.keys()
+                    if key not in self.ignores])
+
+    def checkDirectory(self, dir_name):
+        """
+        Checks all python files in the supplied directory and subdirectories.
+        """
+        self.pep8.process_options(self.options)
+        self.pep8.input_dir(dir_name)
+        return self.countErrors()
+
+    def test_pep8(self):
+        """
+        Tests all directories containing python files for PEP8 compliance.
+        """
+        self.assertTrue(self.pep8 is not None,
+                        msg="Install Python pep8 module to fully execute testbench!")
+        errors = 0
+        for dir_name in ["test", os.environ["BUILDDIR"]]:
+            errors += self.checkDirectory(dir_name)
+        self.assertEqual(errors, 0)
+
+#    def test_variables(self):
+#        """
+#        Function for checking that all attributes are present.
+#        """
+#        try:
+#            import gloripy
+##            for name in ['X', 'Y', 'B',
+##                         'APR_VAL', 'APR_STD', 'FINAL_VAL', 'FINAL_STD',
+##                         'IG_VAL', 'TRUE_VAL',
+##                         'CARTESIAN', 'GEO', 'SPHERICAL',
+##                         'P', 'T', 'C']:
+##                self.assertTrue(hasattr(j, name), msg=str(name) + " is missing.")
+#        except ImportError:
+#            self.assertTrue(False, msg="could not import gloripy")
+
+#    def test_import(self):
+#        """
+#        Checks that all modules are importable.
+#        """
+#        module_list = ["gloripy",
+#                       "gloripy.pylevel0",
+#                      ]
+#        for module in module_list:
+#            try:
+#                exec("import " + module)
+#                importable = True
+#            except ImportError:
+#                importable = False
+#            self.assertTrue(importable, msg="importing " + module + " failed")
diff --git a/test/test_dataloader.py b/test/test_dataloader.py
index b26bfd4..4ee0fbd 100644
--- a/test/test_dataloader.py
+++ b/test/test_dataloader.py
@@ -7,23 +7,78 @@ Created on Wed Apr 24 07:10:28 2013
 # py.test
 
 import unittest
+import numpy as np
+import pyramid.dataloader as dl
+from numpy import pi
 
 
-class TestSuite(unittest.TestCase):
+# TODO: define test constants somewhere
+# TODO: proper error messages
+# TODO: Docstring
+
+
+class TestCaseDataloader(unittest.TestCase):
     
     def setUp(self):
-        pass
+        self.filename = 'test_dataloader/test_data.txt'
+        self.mag_data = dl.MagDataLLG(self.filename)
     
     def tearDown(self):
-        pass
+        self.filename = None
+        self.mag_data = None
+        
+    def test_filename(self):
+        self.assertEqual(self.mag_data.filename, self.filename)
+        
+    def test_resolution(self):
+        self.assertEqual(self.mag_data.res, 10.0)
+        
+    def test_dimensions(self):
+        self.assertEqual(self.mag_data.dim, (1, 3, 5))
     
-    def test(self):
-        self.assertTrue(True)
+    def test_length(self):
+        self.assertEqual(self.mag_data.length, (10.0, 30.0, 50.0))
         
-    def test_almost(self):
-        self.assertAlmostEqual(0, 0.01, places=1)
+    def test_magnitude(self):
+        test_shape = (1, 3, 5)
+        test_array = np.zeros(test_shape)
+        z_mag = self.mag_data.magnitude[0]
+        y_mag = self.mag_data.magnitude[1]
+        x_mag = self.mag_data.magnitude[2]
+        self.assertEqual(z_mag.shape, test_shape)
+        self.assertEqual(y_mag.shape, test_shape)
+        self.assertEqual(x_mag.shape, test_shape)
+        np.testing.assert_array_equal(z_mag, test_array, 'Testmessage')
+        test_array[:,1,1:4] = np.cos(pi/4)
+        np.testing.assert_array_almost_equal(y_mag, test_array, err_msg='y failure')
+        np.testing.assert_array_almost_equal(x_mag, test_array, err_msg='x failure')
         
         
+        
+        
+        
+    
+#mc.create_hom_mag((3,5),10.0,pi/4,mc.slab,((1,2),(1,3)),'test.txt',plot_mag_distr=True)   
+#        '{:7.6e}'
+#        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_len, y_len, z_len = [data[-1, i]/scale+res/2 for i in range(3)]
+#        x_mag, y_mag, z_mag = [data[:, i].reshape(z_dim, y_dim, x_dim).mean(0) 
+#                               *z_len for i in range(3,6)]
+#        #Reshape in Python and Igor is different, 
+#        #Python fills rows first, Igor columns!
+#        self.filename = filename
+#        self.res = res
+#        self.dim = (x_dim, y_dim, z_dim)
+#        self.length = (x_len, y_len, z_len)
+#        self.magnitude = (x_mag, y_mag, z_mag)
+#        self.assertAlmostEqual(0, 0.01, places=1)
+#        mc.create_hom_mag((10,10),10.0,0,mc.slab,((4,4),(5,5)),'test.txt')
+        
 if __name__ == '__main__':
-    suite = unittest.TestLoader().loadTestsFromTestCase(TestSuite)
+    suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseDataloader)
     unittest.TextTestRunner(verbosity=2).run(suite)
\ No newline at end of file
diff --git a/test/test_dataloader/test_data.txt b/test/test_dataloader/test_data.txt
new file mode 100644
index 0000000..397ffeb
--- /dev/null
+++ b/test/test_dataloader/test_data.txt
@@ -0,0 +1,17 @@
+LLGFileCreator2D: test_data
+    5    3    1
+5.000000e-07   5.000000e-07   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+1.500000e-06   5.000000e-07   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+2.500000e-06   5.000000e-07   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+3.500000e-06   5.000000e-07   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+4.500000e-06   5.000000e-07   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+5.000000e-07   1.500000e-06   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+1.500000e-06   1.500000e-06   5.000000e-07   7.071068e-01   7.071068e-01   0.000000e+00
+2.500000e-06   1.500000e-06   5.000000e-07   7.071068e-01   7.071068e-01   0.000000e+00
+3.500000e-06   1.500000e-06   5.000000e-07   7.071068e-01   7.071068e-01   0.000000e+00
+4.500000e-06   1.500000e-06   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+5.000000e-07   2.500000e-06   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+1.500000e-06   2.500000e-06   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+2.500000e-06   2.500000e-06   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+3.500000e-06   2.500000e-06   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
+4.500000e-06   2.500000e-06   5.000000e-07   0.000000e+00   0.000000e+00   0.000000e+00
\ No newline at end of file
diff --git a/test/testsuite.py b/test/testsuite.py
deleted file mode 100644
index b26bfd4..0000000
--- a/test/testsuite.py
+++ /dev/null
@@ -1,29 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Apr 24 07:10:28 2013
-
-@author: Jan
-"""
-# py.test
-
-import unittest
-
-
-class TestSuite(unittest.TestCase):
-    
-    def setUp(self):
-        pass
-    
-    def tearDown(self):
-        pass
-    
-    def test(self):
-        self.assertTrue(True)
-        
-    def test_almost(self):
-        self.assertAlmostEqual(0, 0.01, places=1)
-        
-        
-if __name__ == '__main__':
-    suite = unittest.TestLoader().loadTestsFromTestCase(TestSuite)
-    unittest.TextTestRunner(verbosity=2).run(suite)
\ No newline at end of file
-- 
GitLab