diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 797a14e024631b37c965adc739534ed1a6365431..61221c7dc68d79f46b8b2abc50a6175a3c8d95e7 100644
--- a/pyramid/magdata.py
+++ b/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:
diff --git a/pyramid/numcore/__init__.py b/pyramid/numcore/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..c4a30520d889c50066a74ef98ade1e1b30dad777 100644
--- a/pyramid/numcore/__init__.py
+++ b/pyramid/numcore/__init__.py
@@ -0,0 +1 @@
+from phase_mag_real import *
diff --git a/pyramid/numcore/phase_mag_real.pyx b/pyramid/numcore/phase_mag_real.pyx
new file mode 100644
index 0000000000000000000000000000000000000000..efb838f3d76ae2382b2e9cfb24e5a9f3b4a485c5
--- /dev/null
+++ b/pyramid/numcore/phase_mag_real.pyx
@@ -0,0 +1,32 @@
+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]
+
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 2939b35d2b2322aa25a2e0278302e4d254e0f1a7..7d4a3f06d538fe7bc5b503fe6003a41570d15cbe 100644
--- a/pyramid/phasemapper.py
+++ b/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
 
diff --git a/setup.py b/setup.py
index 48a121a1d7f23566837c29c8b1b8784376c1978d..b4a8f41b60874848ecf819267f71aba2257b7475 100644
--- a/setup.py
+++ b/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')))
 )