add linear_extrude() function to make.py. rotate_extrude() now takes the number of rotational steps to extrude instead of the angle step size. updated documention in make.py.

9 files changed, 135 insertions, 91 deletions

diff --git a/detectors/lbne.py b/detectors/lbne.py
index d3c12da..8942cb2 100644
--- a/detectors/lbne.py
+++ b/detectors/lbne.py
@@ -21,7 +21,7 @@ def build_lbne(radius, height, nstrings, pmts_per_string, endcap_spacing, physic
     lbne = Geometry()
 
     # outer cylinder
-    cylinder_mesh = segmented_cylinder(radius, height+height/(pmts_per_string-1), theta=(2*np.pi/nstrings)/8, n=200)
+    cylinder_mesh = segmented_cylinder(radius, height+height/(pmts_per_string-1), nsteps=16*nstrings, nsegments=1200)
     cylinder_mesh.vertices = rotate(cylinder_mesh.vertices, np.pi/2, (-1,0,0))
     lbne.add_solid(Solid(cylinder_mesh, water_wcsim, vacuum, black_surface, 0xff0000ff))
 
diff --git a/detectors/sno.py b/detectors/sno.py
index 156dcbe..3e748b4 100644
--- a/detectors/sno.py
+++ b/detectors/sno.py
@@ -9,7 +9,7 @@ from solids import build_8inch_pmt_with_lc
 from transform import rotate, make_rotation_matrix
 import os
 
-def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20):
+def sno_vessel(sphere_radius, neck_radius, neck_top, nsteps=40):
     '''Compute the 2D coordinates of the profile of one side of a 
     SNO-style acrylic vessel.  The center of the sphere is at (0,0), with
     the neck extending along the positive y direction.
@@ -17,7 +17,7 @@ def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20):
        sphere_radius: Radius of spherical part of profile
        neck_radius: Radius of neck part
        neck_top: y coordinate of top of neck
-       angle_step: angular step size (radius) between points on the sphere
+       nsteps: number of points around the circumference of the sphere
 
        Returns: Tuple of x and y coordinate numpy arrays.
     '''
@@ -31,7 +31,7 @@ def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20):
         raise ValueError('neck_top must be greater than the y-value where the sphere and cylinder intersect')
 
     # Start with point at bottom
-    angles = np.arange(-math.pi/2, max_angle, angle_step)
+    angles = np.linspace(-math.pi/2, max_angle, nsteps)
     x = list(np.cos(angles) * sphere_radius)
     y = list(np.sin(angles) * sphere_radius)
     x[0] = 0.0 # Round-off error might make cos(-pi/2) not exactly zero
@@ -52,19 +52,19 @@ def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20):
 
 ##### SNO Parts
 
-angle_step = math.pi/20
+nsteps = 40
 av_outside_profile = sno_vessel(sphere_radius=6.0604, neck_radius=0.79375, 
-                                neck_top=10.50, angle_step=angle_step)
+                                neck_top=10.50, nsteps=nsteps)
 # For simplicity, cap the top of the AV with acrylic
 av_inside_profile  = sno_vessel(sphere_radius=6.0053, neck_radius=0.72898,
-                                neck_top=10.00, angle_step=angle_step)
+                                neck_top=10.00, nsteps=nsteps)
 
 av_outside_mesh = rotate_extrude(av_outside_profile[0], av_outside_profile[1],
-                                 angle_step)
+                                 nsteps)
 av_outside_mesh.vertices = rotate(av_outside_mesh.vertices, np.pi/2, (-1,0,0))
 
 av_inside_mesh = rotate_extrude(av_inside_profile[0], av_inside_profile[1],
-                                 angle_step)
+                                 nsteps)
 av_inside_mesh.vertices = rotate(av_inside_mesh.vertices, np.pi/2, (-1,0,0))
 
 dir = os.path.split(os.path.realpath(__file__))[0]
diff --git a/itertoolset.py b/itertoolset.py
index 9d5200c..a20d03b 100644
--- a/itertoolset.py
+++ b/itertoolset.py
@@ -47,3 +47,7 @@ def consume(iterator, n=None):
     else:
         # advance to the empty slice starting at position n
         next(islice(iterator, n, n), None)
+
+def flatten(listOfLists):
+    "Flatten one level of nesting"
+    return chain.from_iterable(listOfLists)
diff --git a/make.py b/make.py
index ed02cb1..1da5d0e 100644
--- a/make.py
+++ b/make.py
@@ -1,67 +1,109 @@
 import numpy as np
 from geometry import Mesh
 from transform import rotate
+from itertoolset import *
 
-def rotate_extrude(x, y, theta=np.pi/32, remove_duplicate_vertices=True):
-    x, y, = np.asarray(x), np.asarray(y)
+def mesh_grid(grid):
+    return np.vstack(zip(grid[:-1].flatten(),grid[1:].flatten(),np.roll(grid[1:],-1,1).flatten()) + zip(grid[:-1].flatten(),np.roll(grid[1:],-1,1).flatten(),np.roll(grid[:-1],-1,1).flatten()))
 
-    if len(x.shape) != 1 or len(y.shape) != 1 or x.size != y.size:
-        raise ValueError('shape mismatch')
+def linear_extrude(x1, y1, height, x2=None, y2=None):
+    """
+    Return the solid mesh formed by linearly extruding the polygon formed by
+    the x and y points `x1` and `y1` by a distance `height`. If `x2` and `y2`
+    are given extrude by connecting the points `x1` and `y1` to `x2` and `y2`;
+    this allows the creation of tapered solids.
 
-    points = np.zeros((len(x),3), dtype=np.float32)
-    points[:,0] = x
-    points[:,1] = y
+    .. note::
+        The path traced by the points `x` and `y` should go counter-clockwise,
+        otherwise the mesh will be inside out.
 
-    angles = np.arange(0, 2*np.pi, theta)
+    Example:
+        >>> # create a hexagon prism
+        >>> angles = np.linspace(0, 2*np.pi, 6, endpoint=False)
+        >>> m = linear_extrude(np.cos(angles), np.sin(angles), 5.0)
+    """
+    if len(x1) != len(y1):
+        raise Exception('`x` and `y` arrays must have the same length.')
 
-    vertices = np.empty((len(points)*len(angles), 3), dtype=np.float32)
-    triangles = np.zeros((len(angles)*(len(points)-1)*2, 3), dtype=np.int32)
+    if x2 is None:
+        x2 = x1
 
-    step = len(points) - 1
+    if y2 is None:
+        y2 = y1
 
-    for i, angle in enumerate(angles):
-        this_slice = i*len(points)
-        vertices[this_slice:this_slice+len(points)] = rotate(points, angle, (0,-1,0))
+    if len(x2) != len(y2) or len(x2) != len(x1):
+        raise Exception('`x` and `y` arrays must have the same length.')
 
-        start = 2*i*step
-        next_slice = ((i+1) % len(angles))*len(points)
+    n = len(x1)
 
-        triangles[start:start+step,0] = np.arange(this_slice, this_slice+len(points)-1)
-        triangles[start:start+step,1] = np.arange(next_slice, next_slice+len(points)-1)
-        triangles[start:start+step,2] = np.arange(this_slice+1, this_slice+len(points))
+    vertex_iterators = [izip(repeat(0,n),repeat(0,n),repeat(-height/2.0,n)),izip(x1,y1,repeat(-height/2.0,n)),izip(x2,y2,repeat(height/2.0,n)),izip(repeat(0,n),repeat(0,n),repeat(height/2.0,n))]
 
-        triangles[start+step:start+2*step,0] = np.arange(next_slice, next_slice+len(points)-1)
-        triangles[start+step:start+2*step,2] = np.arange(this_slice+1, this_slice+len(points))
-        triangles[start+step:start+2*step,1] = np.arange(next_slice+1, next_slice+len(points))
+    vertices = np.fromiter(flatten(roundrobin(*vertex_iterators)), float)
+    vertices = vertices.reshape((len(vertices)//3,3))
 
-    return Mesh(vertices, triangles, remove_duplicate_vertices=remove_duplicate_vertices)
+    triangles = mesh_grid(np.arange(len(vertices)).reshape((len(x1),len(vertices)//len(x1))).transpose()[::-1])
+
+    return Mesh(vertices, triangles, remove_duplicate_vertices=True)
+
+def rotate_extrude(x, y, nsteps=64):
+    """
+    Return the solid mesh formed by extruding the profile defined by the x and
+    y points `x` and `y` around the y axis.
+
+    .. note::
+        The path traced by the points `x` and `y` should go counter-clockwise,
+        otherwise the mesh will be inside out.
+    """
+    if len(x) != len(y):
+        raise Exception('`x` and `y` arrays must have the same length.')
+
+    points = np.array([x,y,np.zeros(len(x))]).transpose()
+
+    steps = np.linspace(0, 2*np.pi, nsteps, endpoint=False)
+    vertices = np.vstack([rotate(points,angle,(0,-1,0)) for angle in steps])
+    triangles = mesh_grid(np.arange(len(vertices)).reshape((len(steps),len(points))).transpose()[::-1])
+
+    return Mesh(vertices, triangles, remove_duplicate_vertices=True)
 
 def cube(size=1):
-    a = np.sqrt(size**2/2.0)
-    x = [0,a,a,0]
-    y = [-size/2.0,-size/2.0,size/2.0,size/2.0]
-    return rotate_extrude(x, y, np.pi/2)
-
-def cylinder(radius1=1, radius2=1, height=2, theta=np.pi/32):
-    x = [0,radius1,radius2,0]
-    y = [-height/2.0, -height/2.0, height/2.0, height/2.0]
-    return rotate_extrude(x, y, theta)
-
-def segmented_cylinder(radius, height=2, theta=np.pi/32, n=50):
-    x = np.concatenate((np.linspace(0, radius, n, endpoint=False), 
-                        [radius] * n,
-                        np.linspace(radius, 0, n, endpoint=False),
-                       [0.0]))
-    y = np.concatenate(([-height/2.0] * n,
-                        np.linspace(-height/2.0, height/2.0, n, endpoint=False),
-                        [height/2.0] * (n+1)))
-    return rotate_extrude(x, y, theta)
-
-def sphere(radius=1, theta=np.pi/32):
-    profile_angles = np.arange(-np.pi/2, np.pi/2+theta, theta)
-    return rotate_extrude(radius*np.cos(profile_angles), radius*np.sin(profile_angles), theta)
-
-def torus(radius=1, offset=3, phi=np.pi/32, theta=np.pi/32):
-    profile_angles = np.arange(0, 2*np.pi+phi, phi)
-    x, y = np.cos(profile_angles), np.sin(profile_angles)
-    return rotate_extrude(x + offset, y)
+    "Return a cube mesh whose sides have length `size`."
+    return linear_extrude([-size/2.0,size/2.0,size/2.0,-size/2.0],[-size/2.0,-size/2.0,size/2.0,size/2.0], height=size)
+
+def cylinder(radius=1, height=2, radius2=None, nsteps=64):
+    """
+    Return a cylinder mesh with a radius of length `radius`, and a height of
+    length `height`. If `radius2` is specified, return a cone shaped cylinder
+    with bottom radius `radius`, and top radius `radius2`.
+    """
+    if radius2 is None:
+        radius2 = radius
+
+    return rotate_extrude([0,radius,radius2,0], [-height/2.0, -height/2.0, height/2.0, height/2.0], nsteps)
+
+def segmented_cylinder(radius, height=2, nsteps=64, nsegments=100):
+    """
+    Return a cylinder mesh segmented into `nsegments` points along its profile.
+    """
+    nsegments_radius = int((nsegments*radius/(2*radius+height))/2)
+    nsegments_height = int((nsegments*height/(2*radius+height))/2)
+    x = np.concatenate([np.linspace(0,radius,nsegments_radius,endpoint=False),[radius]*nsegments_height,np.linspace(radius,0,nsegments_radius,endpoint=False),[0]])
+    y = np.concatenate([[-height/2.0]*nsegments_radius,np.linspace(-height/2.0,height/2.0,nsegments_height,endpoint=False),[height/2.0]*(nsegments_radius+1)])
+    return rotate_extrude(x, y, nsteps)
+
+def sphere(radius=1, nsteps=64):
+    "Return a sphere mesh."
+    profile_angles = np.linspace(-np.pi/2, np.pi/2, nsteps)
+    return rotate_extrude(radius*np.cos(profile_angles), radius*np.sin(profile_angles), nsteps)
+
+def torus(radius=1, offset=3, nsteps=64, circle_steps=None):
+    """
+    Return a torus mesh. `offset` is the distance from the center of the torus
+    to the center of the torus barrel. `radius` is the radius of the torus
+    barrel. `nsteps` is the number of steps in the rotational extrusion of the
+    circle. `circle_steps` if specified is the number of steps around the
+    circumference of the torus barrel, else it defaults to `nsteps`.
+    """
+    if circle_steps is None:
+        circle_steps = nsteps
+    profile_angles = np.linspace(0, 2*np.pi, circle_steps)
+    return rotate_extrude(np.cos(profile_angles) + offset, np.sin(profile_angles), nsteps)
diff --git a/solids/__init__.py b/solids/__init__.py
index 46f3119..a24ce9c 100644
--- a/solids/__init__.py
+++ b/solids/__init__.py
@@ -12,11 +12,11 @@ from optics import *
 from view import buildable
 
 @buildable('12inch_pmt')
-def build_12inch_pmt(outer_material=water, theta=np.pi/8):
-    return build_pmt(dir + '/hamamatsu_12inch.txt', 0.003, outer_material, theta)
+def build_12inch_pmt(outer_material=water, nsteps=16):
+    return build_pmt(dir + '/hamamatsu_12inch.txt', 0.003, outer_material, nsteps)
 
 @buildable('12inch_pmt_shell')
-def build_12inch_pmt_shell(outer_material=water, theta=np.pi/8):
+def build_12inch_pmt_shell(outer_material=water, nsteps=16):
     return build_pmt_shell(dir + '/hamamatsu_12inch.txt')
 
 # from Jelena Maricic
@@ -27,22 +27,22 @@ lc_12inch_rmin = 0.1524
 lc_12inch_rmax = 0.235072
 
 @buildable('12inch_pmt_with_lc')
-def build_12inch_pmt_with_lc(outer_material=water, theta=np.pi/8):
-    pmt = build_12inch_pmt(outer_material, theta)
+def build_12inch_pmt_with_lc(outer_material=water, nsteps=16):
+    pmt = build_12inch_pmt(outer_material, nsteps)
     return pmt + build_light_collector(pmt, a=lc_12inch_a, b=lc_12inch_b, d=lc_12inch_d, rmin=lc_12inch_rmin, rmax=lc_12inch_rmax)
 
 @buildable('12inch_pmt_with_lc_hd')
-def build_12inch_pmt_with_lc_hd(outer_material=water, theta=np.pi/64):
-    pmt = build_12inch_pmt(outer_material, theta)
+def build_12inch_pmt_with_lc_hd(outer_material=water, nsteps=128):
+    pmt = build_12inch_pmt(outer_material, nsteps)
     return pmt + build_light_collector(pmt, a=lc_12inch_a, b=lc_12inch_b, d=lc_12inch_d, rmin=lc_12inch_rmin, rmax=lc_12inch_rmax, npoints=100)
 
 
 @buildable('8inch_pmt')
-def build_8inch_pmt(outer_material=water, theta=np.pi/12):
-    return build_pmt(dir + '/sno_pmt_reduced.txt', 0.003, outer_material, theta)
+def build_8inch_pmt(outer_material=water, nsteps=24):
+    return build_pmt(dir + '/sno_pmt_reduced.txt', 0.003, outer_material, nsteps)
 
 @buildable('8inch_pmt_with_lc')
-def build_8inch_pmt_with_lc(outer_material=water, theta=np.pi/12):
-    pmt = build_8inch_pmt(outer_material, theta)
-    lc = build_light_collector_from_file(dir + '/sno_cone.txt', outer_material, theta)
+def build_8inch_pmt_with_lc(outer_material=water, nsteps=24):
+    pmt = build_8inch_pmt(outer_material, nsteps)
+    lc = build_light_collector_from_file(dir + '/sno_cone.txt', outer_material, nsteps)
     return pmt + lc
diff --git a/solids/pmts.py b/solids/pmts.py
index a30d007..fcbc8fc 100644
--- a/solids/pmts.py
+++ b/solids/pmts.py
@@ -25,11 +25,11 @@ def build_light_collector(pmt, a, b, d, rmin, rmax, npoints=10):
 
     lc_offset = np.interp(lc_radii[0], list(reversed(pmt_face_profile[:,0])), list(reversed(pmt_face_profile[:,1])))
 
-    lc_mesh = rotate_extrude(lc_radii, lc_profile + lc_offset, pmt.theta)
+    lc_mesh = rotate_extrude(lc_radii, lc_profile + lc_offset, pmt.nsteps)
 
     return Solid(lc_mesh, pmt.outer_material, pmt.outer_material, surface=shiny_surface)
 
-def build_pmt_shell(filename, outer_material=water, theta=np.pi/8):
+def build_pmt_shell(filename, outer_material=water, nsteps=16):
     profile = read_csv(filename)
 
     # slice profile in half
@@ -44,9 +44,9 @@ def build_pmt_shell(filename, outer_material=water, theta=np.pi/8):
     # convert mm -> m
     profile /= 1000.0
 
-    return Solid(rotate_extrude(profile[:,0], profile[:,1], theta), glass, outer_material, color=0xeeffffff)
+    return Solid(rotate_extrude(profile[:,0], profile[:,1], nsteps), glass, outer_material, color=0xeeffffff)
 
-def build_pmt(filename, glass_thickness, outer_material=water, theta=np.pi/8):
+def build_pmt(filename, glass_thickness, outer_material=water, nsteps=16):
     profile = read_csv(filename)
 
     # slice profile in half
@@ -63,8 +63,8 @@ def build_pmt(filename, glass_thickness, outer_material=water, theta=np.pi/8):
 
     offset_profile = offset(profile, -glass_thickness)
 
-    outer_envelope_mesh = rotate_extrude(profile[:,0], profile[:,1], theta)
-    inner_envelope_mesh = rotate_extrude(offset_profile[:,0], offset_profile[:,1], theta)
+    outer_envelope_mesh = rotate_extrude(profile[:,0], profile[:,1], nsteps)
+    inner_envelope_mesh = rotate_extrude(offset_profile[:,0], offset_profile[:,1], nsteps)
 
     outer_envelope = Solid(outer_envelope_mesh, glass, outer_material)
 
@@ -78,16 +78,16 @@ def build_pmt(filename, glass_thickness, outer_material=water, theta=np.pi/8):
     # light collector
     pmt.profile = profile
     pmt.outer_material = outer_material
-    pmt.theta = theta
+    pmt.nsteps = nsteps
 
     return pmt
 
-def build_light_collector_from_file(filename, outer_material, theta=np.pi/24):
+def build_light_collector_from_file(filename, outer_material, nsteps=48):
     profile = read_csv(filename)
     
     # Convert mm to m
     profile /= 1000.0
 
-    mesh = rotate_extrude(profile[:,0], profile[:,1], theta)
+    mesh = rotate_extrude(profile[:,0], profile[:,1], nsteps)
     solid = Solid(mesh, outer_material, outer_material, surface=shiny_surface)
     return solid
diff --git a/src/intersect.h b/src/intersect.h
index d2f54ce..26e1d7e 100644
--- a/src/intersect.h
+++ b/src/intersect.h
@@ -60,17 +60,20 @@ __device__ bool intersect_triangle(const float3 &origin, const float3 &direction
    angle between the ray and the plane normal to determine the brightness.
 
    `direction` must be normalized. */
-__device__ unsigned int get_color(const float3 &direction, const float3 &v0, const float3& v1, const float3 &v2, const unsigned int base_color=0xFFFFFFFF)
+__device__ unsigned int get_color(const float3 &direction, const float3 &v0, const float3& v1, const float3 &v2, unsigned int base_color=0xFFFFFF)
 {
 	float scale = dot(normalize(cross(v1-v0,v2-v1)),-direction);
 
+	if (scale < 0.0f)
+	{
+		base_color = 0xff0000;
+		scale = dot(-normalize(cross(v1-v0,v2-v1)),-direction);
+	}
+
 	unsigned int r = 0xFF & (base_color >> 16);
 	unsigned int g = 0xFF & (base_color >> 8);
 	unsigned int b = 0xFF & base_color;
 
-	if (scale < 0.0f)
-		scale = dot(-normalize(cross(v1-v0,v2-v1)),-direction);
-
 	r = floorf(r*scale);
 	g = floorf(g*scale);
 	b = floorf(b*scale);
diff --git a/src/photon.h b/src/photon.h
index fa3b85b..f76ca54 100644
--- a/src/photon.h
+++ b/src/photon.h
@@ -103,9 +103,6 @@ __device__ void fill_state(State &s, Photon &p)
 	s.absorption_length = interp_property(p.wavelength, material1.absorption_length);
 	s.scattering_length = interp_property(p.wavelength, material1.scattering_length);
 
-	printf("wavelength = %f", p.wavelength);
-	printf("scattering length = %f\n", s.scattering_length);
-
 } // fill_state
 
 __device__ void rayleigh_scatter(Photon &p, curandState &rng)
@@ -156,8 +153,6 @@ __device__ int propagate_to_boundary(Photon &p, State &s, curandState &rng)
 	{
 		if (scattering_distance <= s.distance_to_boundary)
 		{
-			printf("scattering distance = %f\n", scattering_distance);
-
 			p.time += scattering_distance/(SPEED_OF_LIGHT/s.refractive_index1);
 			p.position += scattering_distance*p.direction;
 
diff --git a/view.py b/view.py
index b10fbda..7ae6a27 100755
--- a/view.py
+++ b/view.py
@@ -206,7 +206,7 @@ def view(viewable, size=(800,600), name='', bits=8, load_bvh=False):
         cuda.Context.synchronize()
         elapsed = time.time() - t0
 
-        print 'elapsed %f sec' % elapsed
+        #print 'elapsed %f sec' % elapsed
 
         pygame.surfarray.blit_array(screen, pixels_gpu.get().reshape(size))
         pygame.display.flip()