Skip to content

Commit

Permalink
allow both preflattened and non-flattened calls
Browse files Browse the repository at this point in the history
  • Loading branch information
JohannesBuchner committed Feb 17, 2018
1 parent 2400544 commit 6dfeddc
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 27 deletions.
61 changes: 50 additions & 11 deletions raytrace.py
Expand Up @@ -111,7 +111,7 @@ def sphere_raytrace_count_between(xx, yy, zz, RR, rho, a, b, c):
ndpointer(dtype=numpy.float64, ndim=1, flags='C_CONTIGUOUS'),
]

def grid_raytrace(rho, x, y, z, a, b, c):
def grid_raytrace_flat(rho_flat, lenrho, x, y, z, a, b, c):
"""
ray tracing on a grid
Expand All @@ -128,16 +128,34 @@ def grid_raytrace(rho, x, y, z, a, b, c):
NHout double array: output; of size m
"""

lenrho = len(rho)
#rho_flat = numpy.array(rho.flatten())
lena = len(a)
NHout = numpy.zeros(shape=lena) - 1
args = [rho, lenrho, x, y, z, a, b, c, lena, NHout]
r = lib.grid_raytrace(*args)
r = lib.grid_raytrace(rho_flat, lenrho, x, y, z, a, b, c, lena, NHout)
if r != 0:
raise Exception("Calculation failed")
return NHout

def grid_raytrace(rho, x, y, z, a, b, c):
"""
ray tracing on a grid
Parameters regarding the spheres:
rho: double array: density for conversion from length to column density
* n: length of rho
x: double array: start vector
y: double array: start vector
z: double array: start vector
a: double array: direction vector
b: double array: direction vector
c: double array: direction vector
* m: length of a, b, c, x, y, z
NHout double array: output; of size m
"""

lenrho = len(rho)
rho_flat = numpy.array(rho.flatten())
return grid_raytrace_flat(rho_flat, lenrho, x, y, z, a, b, c)


lib.voronoi_raytrace.argtypes = [
ndpointer(dtype=numpy.float64, ndim=1, flags='C_CONTIGUOUS'),
Expand Down Expand Up @@ -364,13 +382,13 @@ def cone_raytrace_finite(thetas, rhos, x, y, z, a, b, c, d):
ndpointer(dtype=numpy.float64, ndim=1, flags='C_CONTIGUOUS'),
]

def grid_raytrace_finite(rho, x, y, z, a, b, c, d):
def grid_raytrace_finite_flat(rho_flat, lenrho, x, y, z, a, b, c, d):
"""
ray tracing on a grid
Parameters regarding the spheres:
rho: double array: density for conversion from length to column density
* n: length of rho
lenrho: length of rho
x: double array: start vector
y: double array: start vector
z: double array: start vector
Expand All @@ -384,8 +402,6 @@ def grid_raytrace_finite(rho, x, y, z, a, b, c, d):
t double array: end position along direction vector. -1 if infinite
"""

lenrho = len(rho)
#rho_flat = numpy.array(rho.flatten())
lena = len(a)
assert len(b) == lena
assert len(c) == lena
Expand All @@ -394,10 +410,33 @@ def grid_raytrace_finite(rho, x, y, z, a, b, c, d):
assert len(z) == lena
assert len(d) == lena
t = numpy.zeros(shape=lena)
args = [rho, lenrho, x, y, z, a, b, c, lena, d, t]
r = lib.grid_raytrace_finite(*args)
r = lib.grid_raytrace_finite(rho_flat, lenrho, x, y, z, a, b, c, lena, d, t)
if r != 0:
raise Exception("Calculation failed")
return t

def grid_raytrace_finite(rho, x, y, z, a, b, c, d):
"""
ray tracing on a grid
Parameters regarding the spheres:
rho: double array: density for conversion from length to column density
* n: length of rho
x: double array: start vector
y: double array: start vector
z: double array: start vector
a: double array: direction vector
b: double array: direction vector
c: double array: direction vector
* m: length of a, b, c, x, y, z
NHmax double array: stop at this NH
Returns:
t double array: end position along direction vector. -1 if infinite
"""

lenrho = len(rho)
rho_flat = numpy.array(rho.flatten())
return grid_raytrace_finite_flat(rho_flat, lenrho, x, y, z, a, b, c, d)


32 changes: 16 additions & 16 deletions test_grid.py
Expand Up @@ -9,7 +9,7 @@ def raytrace_grid_finite(x, v, d):
x, y, z = numpy.array([x0*1. + 128]), numpy.array([y0*1.+128]), numpy.array([z0*1.+128])
a, b, c = numpy.array([dx*1.]), numpy.array([dy*1.]), numpy.array([dz*1.])
NHmax = numpy.array([d*1.])
r = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
r = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
t = r[0]
print('distance:', t, (x, y, z), NHmax)
return x0 + dx * t, y0 + dy * t, z0 + dz * t
Expand Down Expand Up @@ -232,7 +232,7 @@ def test_random_single():
print('direction:', a**2 + b**2 + c**2)
NHmax = numpy.array([63.1791943441])
rho = numpy.ones((256, 256, 256))
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
print('distance:', t)
assert (t > 0).all(), t
assert numpy.isclose(t, NHmax), (t, NHmax)
Expand All @@ -246,11 +246,11 @@ def test_single_small1d_pos():
print('direction:', a**2 + b**2 + c**2)
NHmax = numpy.array([1000.])
rho = numpy.ones((256, 256, 256))
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
print('distance:', t)
assert t == -1, t
NHmax = numpy.array([10.])
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
assert (t > 0).all(), t
assert numpy.isclose(t, NHmax), (t, NHmax)

Expand All @@ -260,11 +260,11 @@ def test_single_small1d_neg():
print('direction:', a**2 + b**2 + c**2)
NHmax = numpy.array([1000.])
rho = numpy.ones((256, 256, 256))
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
print('distance:', t)
assert t == -1, t
NHmax = numpy.array([10.])
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
assert (t > 0).all(), t
assert numpy.isclose(t, NHmax), (t, NHmax)

Expand All @@ -277,7 +277,7 @@ def test_inhomogeneous_single():
a, b, c = numpy.array([-0.321689]), numpy.array([-0.901964]), numpy.array([-0.288056])
print('direction:', a**2 + b**2 + c**2)
NHmax = numpy.array([(10 + 1 + 1) * 3**0.5])
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
print('distance:', t)
assert (t == -1).all(), t

Expand All @@ -290,37 +290,37 @@ def test_density_coord():
b = numpy.array([0., 1., 0., 0., -1., 0.])
c = numpy.array([0., 0., 1., 0., 0., -1.])
NHmax = numpy.array([3.]*6)
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
t_expect = numpy.array([-1,-1,3., -1,-1,-1])
assert numpy.allclose(t, t_expect), (t, t_expect)

rho = numpy.zeros((256, 256, 256))
rho[128,128:140,128] = 1
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
t_expect = numpy.array([-1,3.,-1, -1,-1,-1])
assert numpy.allclose(t, t_expect), (t, t_expect)

rho = numpy.zeros((256, 256, 256))
rho[128:140,128,128] = 1
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
t_expect = numpy.array([3.,-1,-1, -1,-1,-1])
assert numpy.allclose(t, t_expect), (t, t_expect)

rho = numpy.zeros((256, 256, 256))
rho[128,128,120:129] = 1
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
t_expect = numpy.array([-1,-1,-1, -1,-1,3.])
assert numpy.allclose(t, t_expect), (t, t_expect)

rho = numpy.zeros((256, 256, 256))
rho[128,120:129,128] = 1
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
t_expect = numpy.array([-1,-1,-1, -1,3.,-1])
assert numpy.allclose(t, t_expect), (t, t_expect)

rho = numpy.zeros((256, 256, 256))
rho[120:129,128,128] = 1
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
t_expect = numpy.array([-1,-1,-1, 3.,-1,-1])
assert numpy.allclose(t, t_expect), (t, t_expect)

Expand All @@ -338,15 +338,15 @@ def test_inhomogeneous():
a, b, c = numpy.array(a), numpy.array(b), numpy.array(c)
# call raytrace_grid_finite_c()
NHmax = 4.99 + numpy.zeros((N))
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
print('distance:', t)
assert numpy.allclose(t, 0.499), t

# go one further
rho[124:132,124:132,124:132] = 1
rho[128,128,128] = 10
NHmax = 5.99 + numpy.zeros((N))
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
print('distance:', t)
assert (t <= 1.499).all(), t
assert (t >= 0.599).all(), t
Expand All @@ -357,7 +357,7 @@ def test_inhomogeneous():
rho[128,128,128] = 10
NHmax = (10 + 1 + 1) * 3**0.5 + numpy.zeros((N))
#NHmax = 400 + numpy.zeros((20))
t = grid_raytrace_finite(numpy.array(rho.flatten()), x, y, z, a, b, c, NHmax)
t = grid_raytrace_finite(rho, x, y, z, a, b, c, NHmax)
print('distance:', t)
assert (t == -1).all(), t

Expand Down

0 comments on commit 6dfeddc

Please sign in to comment.