Merge

pyramid/analytic.py

@@ -64,10 +64,10 @@ def phase_mag_disc(dim, res, beta, center, radius, height, b_0=1):
     '''
     # Function for the phase:
     def phiMag(x, y):
-        r = np.hypot(x-x0, y-y0)
+        r = np.hypot(x - x0, y - y0)
         r[center[1], center[2]] = 1E-30
-        result = coeff * Lz * ((y-y0) * np.cos(beta) - (x-x0) * np.sin(beta))
-        result *= np.where(r <= R, 1, (R/r)**2)
+        result = coeff * Lz * ((y - y0) * np.cos(beta) - (x - x0) * np.sin(beta))
+        result *= np.where(r <= R, 1, (R / r) ** 2)
         return result
     # Process input parameters:
     z_dim, y_dim, x_dim = dim
@@ -97,10 +97,10 @@ def phase_mag_sphere(dim, res, beta, center, radius, b_0=1):
     '''
     # Function for the phase:
     def phiMag(x, y):
-        r = np.hypot(x-x0, y-y0)
+        r = np.hypot(x - x0, y - y0)
         r[center[1], center[2]] = 1E-30
-        result = coeff * R**3/r**2 * ((y-y0) * np.cos(beta) - (x-x0) * np.sin(beta))
-        result *= np.where(r > R, 1, (1-(1-(r/R)**2)**(3./2.)))
+        result = coeff * R ** 3 / r ** 2 * ((y - y0) * np.cos(beta) - (x - x0) * np.sin(beta))
+        result *= np.where(r > R, 1, (1 - (1 - (r / R) ** 2) ** (3. / 2.)))
         return result
     # Process input parameters:
     z_dim, y_dim, x_dim = dim
@@ -109,8 +109,8 @@ def phase_mag_sphere(dim, res, beta, center, radius, b_0=1):
     R = res * radius
     coeff = - 2./3. * pi * b_0 / PHI_0
     # Create grid:
-    x = np.linspace(res/2, x_dim*res-res/2, num=x_dim)
-    y = np.linspace(res/2, y_dim*res-res/2, num=y_dim)
+    x = np.linspace(res / 2, x_dim * res - res / 2, num=x_dim)
+    y = np.linspace(res / 2, y_dim * res - res / 2, num=y_dim)
     xx, yy = np.meshgrid(x, y)
     # Return phase:
     return phiMag(xx, yy)

pyramid/magdata.py

@@ -5,7 +5,6 @@
 import numpy as np
 import tables.netcdf3 as nc
 import matplotlib.pyplot as plt
-from mayavi import mlab
 
 
 class MagData:
@@ -164,6 +163,8 @@ class MagData:
             None
 
         '''
+        from mayavi import mlab
+
         res = self.res
         dim = self.dim
         # Create points and vector components as lists:

pyramid/numcore/__init__.py

+from phase_mag_real import *

pyramid/numcore/c1.pyx

 import math
 
-def great_circle(float lon1,float lat1,float lon2,float lat2):
+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
+    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)
+    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
+    return radius*c

pyramid/numcore/phase_mag_real.pyx

+import numpy as np
+import math
+cimport cython
+cimport numpy as np
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def phase_mag_real_helper_1(
+        unsigned int v_dim, unsigned int u_dim,
+        double[:, :] phi_u,
+        double[:, :] phi_v,
+        double[:, :] u_mag,
+        double[:, :] v_mag,
+        double[:, :] phase, float threshold):
+    cdef unsigned int i, j, ii, jj, iii, jjj
+    cdef double u, v
+    for j in range(v_dim):
+        for i in range(u_dim):
+            u = u_mag[j, i]
+            v = v_mag[j, i]
+            iii = u_dim - 1 - i
+            jjj = v_dim - 1 - j
+            if abs(u) > threshold:
+               for jj in range(phase.shape[0]):
+                    for ii in range(phase.shape[1]):
+                        phase[jj, ii] += u * phi_u[jjj + jj, iii + ii]
+            if abs(v) > threshold:
+               for jj in range(phase.shape[0]):
+                    for ii in range(phase.shape[1]):
+                        phase[jj, ii] -= v * phi_v[jjj + jj, iii + ii]
+

pyramid/phasemapper.py

@@ -3,7 +3,8 @@
 
 
 import numpy as np
-
+import numcore
+import time
 
 # Physical constants
 PHI_0 = -2067.83    # magnetic flux in T*nm²
@@ -102,12 +103,19 @@ def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
                     phase -= v_mag[j, i] * phase_v
         ############################### TODO: NUMERICAL CORE  #####################################
     else:  # Without Jacobi matrix (faster)
+##        phasecopy = phase.copy()
+##        start_time = time.time()
+##        numcore.phase_mag_real_helper_1(v_dim, u_dim, phi_u, phi_v, u_mag, v_mag, phasecopy, threshold)
+##        print time.time() - start_time
+##        start_time = time.time()
         for j in range(v_dim):
             for i in range(u_dim):
                 if abs(u_mag[j, i]) > threshold:
                     phase += u_mag[j, i] * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
                 if abs(v_mag[j, i]) > threshold:
                     phase -= v_mag[j, i] * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+##        print time.time() - start_time
+##        print ((phase - phasecopy) ** 2).sum()
     # Return the phase:
     return phase
 

setup.py

@@ -15,10 +15,10 @@ from Cython.Build import cythonize
 
 setup(
     name = 'Pyramid',
-    version = '0.1',    
+    version = '0.1',
     description = 'PYthon based Reconstruction Algorithm for MagnetIc Distributions',
     author = 'Jan Caron',
     author_email = 'j.caron@fz-juelich.de',
-    packages = ['pyramid'],
-    ext_modules = cythonize(glob.glob(os.path.join('pyramid','numcore','*.pyx')))
+    packages = ['pyramid', 'pyramid.numcore'],
+    ext_modules = cythonize(glob.glob(os.path.join('pyramid', 'numcore', '*.pyx')))
 )