diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index cb0f7eb..c646f25 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -1,6 +1,6 @@
 set(CMAKE_SHARED_MODULE_PREFIX)
 
-include_directories(${NUMPY_INCLUDE_DIRS} ${PYTHON_INCLUDE_PATH} ${CMAKE_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/src)
+include_directories(${NUMPY_INCLUDE_DIRS} ${PYTHON_INCLUDE_PATH} ${CMAKE_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/src ${CMAKE_SOURCE_DIR}/python)
 
 IF(CYTHON)
 add_custom_command(
diff --git a/python/_project.pyx b/python/_project.pyx
index 61ca0dc..af65590 100644
--- a/python/_project.pyx
+++ b/python/_project.pyx
@@ -9,7 +9,7 @@ DTYPE=np.float64
 
 __all__=["project_cic","line_of_sight_projection","spherical_projection","DTYPE","interp3d","interp2d"]
 
-cdef extern from "project_tools.hpp" namespace "":
+cdef extern from "project_tool.hpp" namespace "":
 
   DTYPE_t compute_projection(DTYPE_t *vertex_value, DTYPE_t *u, DTYPE_t *u0, DTYPE_t rho)
 
@@ -497,7 +497,7 @@ def tophat_fourier(x not None):
 
 @cython.boundscheck(False)
 @cython.cdivision(True)
-cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]):
+cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]) nogil:
     cdef DTYPE_t alpha_max
     cdef DTYPE_t tmp_a
     cdef DTYPE_t v[3]
@@ -505,7 +505,7 @@ cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1]):
 
     alpha_max = 10.0 # A big number
 
-    for i in range(3):
+    for i in xrange(3):
         if u[i] == 0.:
             continue
 
@@ -539,7 +539,7 @@ cdef DTYPE_t mysum(DTYPE_t *v, int q) nogil:
 
 @cython.boundscheck(False)
 @cython.cdivision(True)
-cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8]):
+cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8]) nogil:
     cdef DTYPE_t alpha_max
     cdef DTYPE_t tmp_a
     cdef DTYPE_t v[3], term[4]
@@ -565,12 +565,33 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t
     # we integrate between 0 and alpha_max (curvilinear coordinates)
     
     return compute_projection(vertex_value, u, u0, alpha_max)
+
+
+@cython.boundscheck(False)
+cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
+                         DTYPE_t u[3], DTYPE_t u0[3], int iu0[3], int jumper[1]) nogil:
+    cdef DTYPE_t d
+    
+    d = density[iu0[0], iu0[1], iu0[2]]
+    
+    return cube_integral(u, u0, jumper)*d
+
+@cython.boundscheck(False)
+cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
+                         DTYPE_t u[3], DTYPE_t u0[3], int iu0[3], int jumper[1]) nogil:
+    cdef DTYPE_t vertex_value[8]
+    cdef DTYPE_t d
+    
+    d = density[iu0[0], iu0[1], iu0[2]]
+    return cube_integral(u, u0, jumper)*d
+    cube_integral_trilin(u, u0, jumper, vertex_value)
+
     
 @cython.boundscheck(False)
 def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
                              npx.ndarray[DTYPE_t] a_u,
                              DTYPE_t min_distance,
-                             DTYPE_t max_distance):
+                             DTYPE_t max_distance, int integrator_id=0):
 
     cdef DTYPE_t u[3], ifu0[3], u0[3], utot[3]
     cdef int iu0[3]
@@ -580,6 +601,14 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
     cdef int completed
     cdef DTYPE_t I0, d, dist2, delta, s, max_distance2
     cdef int jumper[1]
+    
+    cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:],
+                               DTYPE_t u[3], DTYPE_t u0[3], int iu0[3], int jumper[1]) nogil
+
+    if integrator_id == 0:
+      integrator = integrator0
+    else:
+      integrator = integrator1
 
     max_distance2 = max_distance**2
 
@@ -606,15 +635,15 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
     jumper[0] = 0
     while completed == 0:
         
-        d = density[iu0[0], iu0[1], iu0[2]]
-        s = cube_integral(u, u0, jumper)
-        I0 += s*d
+        I0 += integrator(density, u, u0, iu0, jumper)
 
         if u[jumper[0]] < 0:
             iu0[jumper[0]] -= 1
+            direction = -1
             u0[jumper[0]] = 1
         else:
             iu0[jumper[0]] += 1
+            direction = 1
             u0[jumper[0]] = 0
 
             
@@ -630,8 +659,8 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
 
           if (dist2 > max_distance2):
             # Remove the last portion of the integral
-            delta = sqrt(dist2) - max_distance
-            I0 -= d*delta
+            #delta = sqrt(dist2) - max_distance
+            #I0 -= d*delta
             completed = 1
             
     return I0