Changed location of the repository (Desktop -> Daten)

analytic: changed analytic solution for discs (still not perfect)
magcreator: added logo function
phasemap: revised the analytical formulas -> more uniform
added new folder for PYRAMID project
added testsuite module

PYRAMID/desktop.ini

+[.ShellClassInfo]
+IconResource=C:\Users\Jan\Daten\FZ-Jlich\Programs\pyramid\pyramid.ico,0
+[ViewState]
+Mode=
+Vid=
+FolderType=Generic

Python/.spyderworkspace

@@ -2,7 +2,7 @@
 S'project_paths'
 p2
 (lp3
-VC:\u005cUsers\u005cJan\u005cDesktop\u005cFZ-Jlich\u005cPrograms\u005cPython\u005cmagneticimaging
+VC:\u005cUsers\u005cJan\u005cDaten\u005cFZ-Jlich\u005cPrograms\u005cPython\u005cmagneticimaging
 p4
 asS'name'
 p5

Python/magneticimaging/analytic.py

@@ -66,17 +66,26 @@ def phasemap_slab(dim, res, beta, center, width, b_0):
 def phasemap_disc(dim, res, beta, center, radius, b_0):
     '''INPUT VARIABLES'''
     y_dim, x_dim = dim
-    y0, x0 = res * center[0], res * center[1]
-    R = res * radius
+    # 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)
+    # TODO: Explanation
+    # Ly, Lx have to be odd, because the slab borders should not lie in the
+    # center of a pixel (so L/2 can't be an integer)
+    R  = res * radius
     
     '''COMPUTATION MAGNETIC PHASE SHIFT (REAL SPACE) DISC'''
       
     coeff = - pi * res * b_0 / ( 2 * PHI_0 )
-      
+    
+#    import pdb; pdb.set_trace()    
+    
     def phiMag(x,y):
-        r = np.hypot(x-x0, y-y0)      
+        r = np.hypot(x-x0, y-y0)
+        r[center[0], center[1]] = 1E-18
         result = coeff * ((y-y0) * np.cos(beta) - (x-x0) * np.sin(beta))
-        in_or_out = 1 * (r <= R) + (R / r) ** 2 * (r > R)    
+        in_or_out = 1 * (r < R) + (R / r) ** 2 * (r > R)
+#        import pdb; pdb.set_trace() 
         result *= in_or_out
         return result
     

Python/magneticimaging/magcreator.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
+from numpy import pi
 
 
 def slab(dim, params):
@@ -34,10 +35,8 @@ def disc(dim, params):
         
     '''
     center, radius = params
-    shape_mag = np.array([[(np.sqrt((x - center[1]) ** 2 +
-                                    (y - center[0]) ** 2) <= radius)
-                                    for x in range(dim[1])] 
-                                    for y in range(dim[0])])    
+    shape_mag = 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
 
     
@@ -140,6 +139,48 @@ def create_hom_mag(dim, res, beta, shape_fun, params,
     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 (True):
+        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', ''))

Python/magneticimaging/phasemap.py

@@ -7,7 +7,11 @@ import matplotlib.pyplot as plt
 from numpy import pi
 
 
-PHI_0 = 2067.83  # magnetic flux in T*nm²
+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
 
 
 def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
@@ -57,7 +61,7 @@ def fourier_space(mag_data, b_0=1, padding=0, v_0=0, v_acc=30000):
     return phase
     
     
-def real_space(mag_data, b_0=1, v_0=0, v_acc=30000):
+def real_space_slab(mag_data, b_0=1, v_0=0, v_acc=30000):
     '''Calculate phasemap from magnetization data (real space approach).
     Arguments:
         mag_data - MagDataLLG object (from magneticimaging.dataloader) storing
@@ -77,76 +81,135 @@ def real_space(mag_data, b_0=1, v_0=0, v_acc=30000):
     
     beta = np.arctan2(y_mag, x_mag)
 
-    abs_mag = np.sqrt(x_mag**2 + y_mag**2)    
+    mag = np.hypot(x_mag, y_mag)    
      
-    coeff = abs_mag * res / ( 4 * PHI_0 ) / res # TODO: Lz = res 
+    coeff = - b_0 * res**2 / ( 2 * PHI_0 ) / res # TODO: why / res
 
-    def F0(n, m):
+    def F_h(n, m):
         a = np.log(res**2 * (n**2 + m**2))
         b = np.arctan(n / m)
-        return res * (n*a - 2*n + 2*m*b)
+        return n*a - 2*n + 2*m*b
      
-    def F_part(n, m):
-        return ( F0(n-0.5, m-0.5) - F0(n+0.5, m-0.5)
-                -F0(n-0.5, m+0.5) + F0(n+0.5, m+0.5) )
-    
-    def phiMag(xx, yy, xi, yj, coeffij, betaij):
-        return coeffij * ( - np.cos(betaij) * F_part(xx-xi, yy-yj)
-                           + np.sin(betaij) * F_part(yy-yj, xx-xi) )
+    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) )
+                
+    def phi_mag(i, j):  # TODO: rename
+        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)
     
-    xF = np.linspace(-(x_dim-1), x_dim-1, num=2*x_dim-1)
-    yF = np.linspace(-(y_dim-1), y_dim-1, num=2*y_dim-1)
-    xxF, yyF = np.meshgrid(xF,yF)
+    phi_cos = phi_pixel(xx_big, yy_big)
+    phi_sin = phi_pixel(yy_big, xx_big)
     
-    F_part_cos = F_part(xxF, yyF)
-    F_part_sin = F_part(yyF, xxF)
+    display_phase(phi_cos, res, 'Phase of one Pixel (Cos - Part)')
+#    display_phase(phi_sin, res, 'Phase of one Pixel (Sin - Part)')
     
-    display_phase(F_part_cos, res, 'F_part_cos')
-    display_phase(F_part_sin, res, 'F_part_sin')      
+    '''CALCULATE THE PHASE'''
+    phase = np.zeros((y_dim, x_dim))
     
-    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):
+            phase += phi_mag(i, j)
     
-    for j in y:
-        for i in x:
-            phase += phiMag(xx, yy, i, j, coeff[j,i], beta[j,i])
-            
-    return phase
+    return (phase, phi_cos)
     
     
-    xF = np.linspace(-(x_dim-1), x_dim-1, num=2*x_dim-1)
-    yF = np.linspace(-(y_dim-1), y_dim-1, num=2*y_dim-1)
-    xxF, yyF = np.meshgrid(xF,yF)
+def real_space_disc(mag_data, b_0=1, v_0=0, v_acc=30000):
+    '''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    
     
-    F_part_cos = F_part(xxF, yyF)
-    F_part_sin = F_part(yyF, xxF)
+    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])
     
-    display_phase(F_part_cos, res, 'F_part_cos')
-    display_phase(F_part_sin, res, 'F_part_sin')    
+    '''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)
     
-    def phiMag2(xx, yy, i, j):
-        #x_ind = xxF[yy]
-        
-        return coeff[j,i] * ( - np.cos(beta[j,i]) * F_part_cos[yy.min()-j:(2*yy.max()-1)-j, xx.min()-i:(2*xx.max()-1)-i]
-                              + np.sin(beta[j,i]) * F_part_sin[xx.min()-i:(2*xx.max()-1)-i, yy.min()-j:(2*yy.max()-1)-j] )
+    phi_cos = phi_pixel(xx_big, yy_big)
+    phi_sin = phi_pixel(yy_big, xx_big)
     
+    display_phase(phi_cos, res, 'Phase of one Pixel (Cos - Part)')
+#    display_phase(phi_sin, res, 'Phase of one Pixel (Sin - Part)')
     
-    phase2 = np.zeros((y_dim*x_dim, y_dim, x_dim))
+    '''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):
-            phase2 += phiMag2(xx, yy, 0, 0)
-                
-   
+            phase += phi_mag(i, j)
     
-    phase2 = np.sum(phase2, axis=0)
+    return (phase, phi_cos)
     
-    return phase2
+
+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
+    x_dim, y_dim, z_dim = mag_data.dim
+    x_mag, y_mag, z_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
     
 	
 def display_phase(phase, res, title, axis=None):
@@ -166,9 +229,9 @@ def display_phase(phase, res, title, axis=None):
     im = plt.pcolormesh(phase, cmap='gray')
 
     ticks = axis.get_xticks()*res
-    axis.set_xticklabels(ticks.astype(int))
+    axis.set_xticklabels(ticks)
     ticks = axis.get_yticks()*res
-    axis.set_yticklabels(ticks.astype(int))
+    axis.set_yticklabels(ticks)
 
     axis.set_title(title)
     axis.set_xlabel('x-axis [nm]')

Python/magneticimaging/testsuite.py

+# -*- 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

Python/main.py

@@ -14,7 +14,7 @@ import magneticimaging.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.
@@ -28,27 +28,27 @@ def phase_from_mag():
     '''INPUT'''
     # TODO: Input via GUI
     filename = 'output.txt'
-    b_0 = 1.0  # in T
+    b_0 = 10.0  # in T
     v_0 = 0  # TODO: units?
     v_acc = 30000  # in V
     padding = 20
-    density = 10
+    density = 100
     
-    dim = (50, 50)  # in px (y,x)
-    res = 10.0  # in nm
-    beta = pi/4
+    dim = (100, 100)  # in px (y,x)
+    res = 1.0  # in nm
+    beta = pi/2
     
     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)
+    center = (49, 49)  # in px (y,x) index starts with 0!
+    width  = (50, 50)  # in px (y,x)
     params = (center, width)
 #    # Disc:
 #    shape_fun = mc.disc
 #    center = (4, 4)  # in px (y,x)
-#    radius = 2
+#    radius = 2.5
 #    params = (center, radius)
 #    # Filament:
 #    shape_fun = mc.filament
@@ -65,6 +65,13 @@ def phase_from_mag():
 #    pixel = (5, 5)
 #    params = pixel
     
+    '''CREATE LOGO'''    
+    mc.create_logo(128, res, beta, filename, plot_mag_distr)
+    mag_data = dl.MagDataLLG(filename)
+    phase, pixel_stub = 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)
@@ -72,27 +79,47 @@ def phase_from_mag():
     '''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)    
+    
+    
+    
     '''COLOR WHEEL'''
     hi.make_color_wheel()
-
+    
+    
+    phase_stub, phi_cos_real_slab = pm.real_space_slab(mag_data, b_0)
+    phase_stub, phi_cos_real_disc = pm.real_space_disc(mag_data, b_0) 
+    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:
-    ticf = time.clock()
+    tic = time.clock()
     phase_fft = pm.fourier_space(mag_data, b_0, padding)   
-    tocf = time.clock()
-    print 'Time for Fourier Space Approach: ' + str(tocf - ticf)
+    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:
-    ticr = time.clock()
-    phase_real = pm.real_space(mag_data, b_0)
-    tocr = time.clock()
-    print 'Time for Real Space Approach:    ' + str(tocr - ticr)
-    print 'Fourier Approach is ' + str((tocr-ticr) / (tocf-ticf)) + ' faster!'
-    holo_real = hi.holo_image(phase_real, mag_data.res, density)
-    display_combined(phase_real, mag_data.res, holo_real, 
-                     'Real Space Approach')
+    # numerical solution Real Space (Slab):
+    tic = time.clock()
+    phase_real_slab, pixel_stub = pm.real_space_slab(mag_data, b_0)
+    toc = time.clock()
+    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, pixel_stub = 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:
@@ -110,12 +137,17 @@ def phase_from_mag():
 
     
     '''DIFFERENCES'''
-    diff_real_to_ana = phase_real - phase
-    diff_fft_to_ana  = phase_fft  - phase 
-    diff_real_to_fft = phase_fft - phase_real
-    pm.display_phase(diff_real_to_ana, res, 'Difference: Analytic - Real')
-    pm.display_phase(diff_fft_to_ana,  res, 'Difference: Analytic - Fourier')
-    pm.display_phase(diff_real_to_fft, res, 'Difference: Fourier - Real') 
+    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):