Further Changes of the project structure

jacobi matrix now contains the derivatives instead of descriptive strings

main.py → scripts/main.py

@@ -35,7 +35,7 @@ def phase_from_mag():
     padding = 20
     density = 100
     
-    dim = (3, 3)  # in px (y,x)
+    dim = (50, 50)  # in px (y,x)
     res = 1.0  # in nm
     beta = pi/4
     
@@ -43,8 +43,8 @@ def phase_from_mag():
     
     # Slab:
     shape_fun = mc.slab
-    center = (2, 2)  # in px (y,x) index starts with 0!
-    width  = (1, 1)  # in px (y,x)
+    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
@@ -105,13 +105,11 @@ def phase_from_mag():
     display_combined(phase_fft, mag_data.res, holo_fft, 
                      'Fourier Space Approach')
     # numerical solution Real Space (Slab):
-    jacobi = np.chararray((2*dim[0]*dim[1], dim[0]*dim[1]), itemsize=10)
+#    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=jacobi)
+    phase_real_slab = pm.real_space_slab(mag_data, b_0, jacobi=None)
     toc = time.clock()
-    np.save('jacobi.npy', jacobi)
-    jacobi = np.load('jacobi.npy')
-    enter_pdb # <- raises error
+#    np.savetxt('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, 

pyramid/phasemap.py → src/pyramid/phasemap.py

@@ -96,27 +96,16 @@ def real_space_slab(mag_data, b_0=1, v_0=0, v_acc=30000,
                               -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])
-
-    Y = np.chararray((y_dim,x_dim))
-    Y[:] = 'y'
-    Div = np.chararray((y_dim,x_dim))
-    Div[:] = '/'
-    Index = np.chararray((y_dim,x_dim), itemsize=3)
-    Index[:] = np.array(range(y_dim*x_dim),dtype='|S3').reshape((y_dim,x_dim))
-    Mag = np.chararray((y_dim,x_dim))
-    Mag[:] = 'm'
-    Beta = np.chararray((y_dim,x_dim))
-    Beta[:] = 'b'
-                                                   
-    def phi_del_mag(i, j):  # TODO: rename
-        return Y + Index + Div + Mag + (i+x_dim*j).astype('|S5')
-                                                   
-    def phi_del_beta(i, j):  # TODO: rename      
-        return Y + Index + Div + Beta + (i+x_dim*j).astype('|S5')
+        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)
@@ -139,14 +128,15 @@ def real_space_slab(mag_data, b_0=1, v_0=0, v_acc=30000,
     # 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):
-            phase += phi_mag(i, j)
-            if jacobi is not None:
-                jacobi[i+x_dim*j]             = phi_del_mag(i, j).reshape(-1)
-                jacobi[y_dim*x_dim+i+x_dim*j] = phi_del_beta(i, j).reshape(-1)
-                
+            if (mag[j, i] != 0): # 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)
     
     return phase
     
@@ -211,7 +201,8 @@ def real_space_disc(mag_data, b_0=1, v_0=0, v_acc=30000,
     # 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)
+#            if (mag[j, i] != 0): # TODO: same result with or without?
+                phase += phi_mag(i, j)
     
     return phase