guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

More projection code

Browse Source
This commit is contained in:
Guilhem Lavaux 2014-06-04 18:06:48 +02:00
parent be03931328
commit 7b7d5b050e
2 changed files with 40 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
python/CMakeLists.txt
View File

@ -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(

49
python/_project.pyx
View File

@ -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]
@ -566,11 +566,32 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t
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]
@ -581,6 +602,14 @@ def line_of_sight_projection(npx.ndarray[DTYPE_t, ndim=3] density,
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
for i in range(3):
@ -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