diff --git a/utils/Makefile.am b/utils/Makefile.am
index beaeeb0..ad2aef0 100644
--- a/utils/Makefile.am
+++ b/utils/Makefile.am
@@ -14,6 +14,13 @@ geomtst_SOURCES = geomtst.c
 geomtst_LDADD = libctlgeom.la
 geomtst_CPPFLAGS = -I$(top_srcdir)/src
 
+check_PROGRAMS = test-prism
+
+test_prism_SOURCES = test-prism.c
+test_prism_LDADD   = libctlgeom.la
+
+TESTS = test-prism
+
 dist_man_MANS = gen-ctl-io.1
 
 BUILT_SOURCES = gen-ctl-io geom-ctl-io.c ctlgeom-types.h nlopt-constants.scm
diff --git a/utils/ctlgeom.h b/utils/ctlgeom.h
index f1f8f85..d963e45 100644
--- a/utils/ctlgeom.h
+++ b/utils/ctlgeom.h
@@ -69,6 +69,9 @@ extern vector3 normal_to_object(vector3 p, GEOMETRIC_OBJECT o);
 extern vector3 normal_to_fixed_object(vector3 p, GEOMETRIC_OBJECT o);
 extern int intersect_line_with_object(vector3 p, vector3 d, GEOMETRIC_OBJECT o,
 				      double s[2]);
+extern double intersect_line_segment_with_object(vector3 p, vector3 d,
+                                                 GEOMETRIC_OBJECT o,
+                                                 double a, double b);
 extern MATERIAL_TYPE material_of_point_inobject(vector3 p, boolean *inobject);
 extern MATERIAL_TYPE material_of_point_inobject0(
      GEOMETRIC_OBJECT_LIST geometry, vector3 p, boolean *inobject);
@@ -144,7 +147,11 @@ GEOMETRIC_OBJECT make_block(MATERIAL_TYPE material, vector3 center,
 GEOMETRIC_OBJECT make_ellipsoid(MATERIAL_TYPE material, vector3 center,
 				vector3 e1, vector3 e2, vector3 e3,
 				vector3 size);
+GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material,
+			    const vector3 *vertices, int num_vertices,
+			    double height, vector3 axis);
 
+int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance);
 
 /**************************************************************************/
 
diff --git a/utils/geom.c b/utils/geom.c
index 411c467..5fd6562 100644
--- a/utils/geom.c
+++ b/utils/geom.c
@@ -21,6 +21,7 @@
 
 #include <stdlib.h>
 #include <stdio.h>
+#include <string.h>
 #include <math.h>
 
 #ifndef LIBCTLGEOM
@@ -59,7 +60,20 @@ using namespace ctlio;
 #endif
 
 #define K_PI 3.14159265358979323846
+#define CHECK(cond, s) if (!(cond)){fprintf(stderr,s "\n");exit(EXIT_FAILURE);}
+
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
 
+// forward declarations of prism-related routines, at the bottom of this file
+static vector3 prism_vector_p2c(prism *prsm, vector3 vp);
+static vector3 prism_vector_c2p(prism *prsm, vector3 vc);
+static void get_prism_bounding_box(prism *prsm, geom_box *box);
+static double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b);
+static boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices);
+static boolean point_in_prism(prism *prsm, vector3 p);
+static vector3 normal_to_prism(prism *prsm, vector3 p);
+static void display_prism_info(int indentby, prism *prsm);
 /**************************************************************************/
 
 /* If v is a vector in the lattice basis, normalize v so that
@@ -123,6 +137,10 @@ void geom_fix_object(geometric_object o)
 	      o.subclass.block_data->projection_matrix = matrix3x3_inverse(m);
 	      break;
 	 }
+	 case GEOM PRISM:
+	 {    
+	 }
+	 break;
 	 case GEOM COMPOUND_GEOMETRIC_OBJECT:
 	 {
 	      int i;
@@ -297,6 +315,10 @@ boolean point_in_fixed_pobjectp(vector3 p, geometric_object *o)
 	}
       }
     }
+  case GEOM PRISM:
+    {
+      return point_in_prism(o->subclass.prism_data, p);
+    }
   case GEOM COMPOUND_GEOMETRIC_OBJECT:
     {
 	 int i;
@@ -348,6 +370,8 @@ vector3 to_geom_object_coords(vector3 p, geometric_object o)
       if (size.z != 0.0) proj.z /= size.z;
       return vector3_plus(half, proj);
     }
+/* case GEOM PRISM:
+    NOT YET IMPLEMENTED */
   }
 }
 
@@ -378,6 +402,8 @@ vector3 from_geom_object_coords(vector3 p, geometric_object o)
 		     vector3_scale(size.z * p.z, o.subclass.block_data->e3))
 		));
     }
+/* case GEOM PRISM:
+    NOT YET IMPLEMENTED */
   }
 }
 
@@ -448,7 +474,11 @@ vector3 normal_to_fixed_object(vector3 p, geometric_object o)
 	  return matrix3x3_transpose_vector3_mult(
 	       o.subclass.block_data->projection_matrix, proj);
 	}
-      }
+      } // switch
+    } // case GEOM BLOCK
+  case GEOM PRISM:
+    { 
+      return normal_to_prism(o.subclass.prism_data, p);
     }
   default:
        return r;
@@ -641,6 +671,9 @@ void CTLIO display_geometric_object_info(int indentby, geometric_object o)
 		       break;
 	      }
 	      break;
+	 case GEOM PRISM:
+	      printf("prism");
+	      break;
 	 case GEOM COMPOUND_GEOMETRIC_OBJECT:
 	      printf("compound object");
 	      break;
@@ -686,6 +719,9 @@ void CTLIO display_geometric_object_info(int indentby, geometric_object o)
                      o.subclass.block_data->e3.y,
                      o.subclass.block_data->e3.z);
 	      break;
+	 case GEOM PRISM:
+	      display_prism_info(indentby, o.subclass.prism_data);
+	      break;
 	 case GEOM COMPOUND_GEOMETRIC_OBJECT:
 	 {
 	      int i;
@@ -711,6 +747,7 @@ void CTLIO display_geometric_object_info(int indentby, geometric_object o)
 int intersect_line_with_object(vector3 p, vector3 d, geometric_object o,
 			       double s[2])
 {
+     vector3 p0=p;
      p = vector3_minus(p, o.center);
      s[0] = s[1] = 0;
      switch (o.which_subclass) {
@@ -734,7 +771,7 @@ int intersect_line_with_object(vector3 p, vector3 d, geometric_object o,
 		   s[1] = (b2 - discrim) / a;
 		   return 2;
 	      }
-	 }
+	 } // case GEOM SPHERE
 	 case GEOM CYLINDER: {
 	      vector3 dm = matrix3x3_vector3_mult(geometry_lattice.metric, d);
 	      vector3 pm = matrix3x3_vector3_mult(geometry_lattice.metric, p);
@@ -796,7 +833,7 @@ int intersect_line_with_object(vector3 p, vector3 d, geometric_object o,
 	      }
 	      if (ret == 2 && s[0] == s[1]) ret = 1;
 	      return ret;
-	 }
+	 } // case GEOM CYLINDER
 	 case GEOM BLOCK:
 	 {
 	      vector3 dproj = matrix3x3_vector3_mult(o.subclass.block_data->projection_matrix, d);
@@ -867,12 +904,34 @@ int intersect_line_with_object(vector3 p, vector3 d, geometric_object o,
 		       }
 		  }
 	      }
-	 }
+	 } // case GEOM BLOCK
+
 	 default:
 	      return 0;
      }
 }
 
+/* Compute the intersections with o of a line along p+s*d in the interval s in [a,b], returning
+    the length of the s intersection in this interval.  (Note: o must not be a compound object.) */
+double intersect_line_segment_with_object(vector3 p, vector3 d, geometric_object o, double a, double b)
+{
+  if (o.which_subclass==GEOM PRISM)
+   { 
+     return intersect_line_segment_with_prism(o.subclass.prism_data, p, d, a, b);
+   }
+  else
+   { double s[2];
+     if (2 == intersect_line_with_object(p, d, o, s)) 
+      { double ds = (s[0] < s[1] ? MIN(s[1],b) - MAX(s[0],a)
+	           	         : MIN(s[0],b) - MAX(s[1],a)
+                    );
+	return (ds > 0 ? ds : 0.0);
+      }
+     else 
+      return 0.0;
+   }
+}
+
 /**************************************************************************/
 
 /* Given a basis (matrix columns are the basis unit vectors) and the
@@ -925,9 +984,6 @@ matrix3x3 CTLIO square_basis(matrix3x3 basis, vector3 size)
 /**************************************************************************/
 
 /* geom_box utilities: */
-
-#define MAX(a,b) ((a) > (b) ? (a) : (b))
-#define MIN(a,b) ((a) < (b) ? (a) : (b))
 static void geom_box_union(geom_box *bu,
 			   const geom_box *b1, const geom_box *b2)
 {
@@ -1140,6 +1196,11 @@ void geom_get_bounding_box(geometric_object o, geom_box *box)
 	        vector3_plus(corner, vector3_plus(s1, vector3_plus(s2, s3))));
 	      break;
 	 }
+	 case GEOM PRISM:
+	 {
+           get_prism_bounding_box(o.subclass.prism_data, box);
+           break;
+	 }
 	 case GEOM COMPOUND_GEOMETRIC_OBJECT:
 	 {
 	      int i;
@@ -1223,13 +1284,7 @@ static double overlap_integrand(integer ndim, number *x, void *data_)
 	  a0 = data->c0 - w; b0 = data->c0 + w;
      }
 
-     if (2 == intersect_line_with_object(p, data->dir, data->o, s)) {
-	  double ds = (s[0] < s[1]
-		       ? MIN(s[1],b0) - MAX(s[0],a0)
-		       : MIN(s[0],b0) - MAX(s[1],a0));
-	  return (ds > 0 ? ds * scale_result : 0.0);
-     }
-     return 0.0;
+    return intersect_line_segment_with_object(p, data->dir, data->o, a0, b0) * scale_result;
 }
 
 number overlap_with_object(geom_box b, int is_ellipsoid, geometric_object o,
@@ -1411,8 +1466,6 @@ static int object_in_box(geometric_object o, vector3 shiftby,
      return geom_boxes_intersect(obj_b, b);
 }
 
-#define CHECK(cond, s) if (!(cond)){fprintf(stderr,s "\n");exit(EXIT_FAILURE);}
-
 static geom_box_tree new_geom_box_tree(void)
 {
      geom_box_tree t;
@@ -2014,3 +2067,451 @@ geometric_object make_ellipsoid(material_type material, vector3 center,
 	  = 2.0 / size.z;
      return o;
 }
+
+/*************************************************************** 
+ * the remainder of this file is the content of `prism.c`
+ * implementing geometric primitives for prisms
+ *
+ * prism.c -- geometry routines for prisms.
+ * a prism is a planar polygon, consisting of 3 or more user-specified
+ * vertices, extruded through a given thickness (the "height") in the
+ * direction of a given axis.
+ * most calculations are done in the "prism coordinate system",
+ * in which the prism floor lies in the XY plane with centroid
+ * at the origin and the prism axis is the positive Z-axis.
+ * homer reid 4/2018 
+ * 
+ ***************************************************************/
+
+/***************************************************************/
+/* given coordinates of a point in the prism coordinate system,*/
+/* return cartesian coordinates of that point                  */
+/***************************************************************/
+vector3 prism_coordinate_p2c(prism *prsm, vector3 vp)
+{ return vector3_plus(prsm->centroid, matrix3x3_vector3_mult(prsm->m_p2c,vp)); }
+
+vector3 prism_vector_p2c(prism *prsm, vector3 vp)
+{ return matrix3x3_vector3_mult(prsm->m_p2c, vp); }
+
+vector3 prism_coordinate_c2p(prism *prsm, vector3 vc)
+{ return matrix3x3_vector3_mult(prsm->m_c2p, vector3_minus(vc,prsm->centroid)); }
+
+vector3 prism_vector_c2p(prism *prsm, vector3 vc)
+{ return matrix3x3_vector3_mult(prsm->m_c2p, vc); }
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+void get_prism_bounding_box(prism *prsm, geom_box *box)
+{
+  vector3 *vertices = prsm->vertices.items;
+  int num_vertices  = prsm->vertices.num_items;
+  double height     = prsm->height;
+
+  // set x,y coordinates of low, high to bounding box 
+  // of prism base polygon (in prism coordinate system)
+  vector3 low=vertices[0], high=vertices[0];
+  int nv;
+  for(nv=1; nv<num_vertices; nv++)
+   { low.x  = fmin(low.x, vertices[nv].x);
+     low.y  = fmin(low.y, vertices[nv].y);
+     high.x = fmax(high.x, vertices[nv].x);
+     high.y = fmax(high.y, vertices[nv].y);
+   };
+
+  // set z coordinates of low, high to upper and lower 
+  // prism surfaces (in prism coordinate system)
+  low.z  = 0.0;
+  high.z = height;
+  
+  // convert from prism coordinates to cartesian coordinates
+  box->low  = prism_coordinate_p2c(prsm, low);
+  box->high = prism_coordinate_p2c(prsm, high);
+}
+
+/***************************************************************/
+/* find the value of s at which the line p+s*d intersects the  */
+/* line segment connecting v1 to v2 (in 2 dimensions)          */
+/* algorithm: solve the 2x2 linear system p+s*d = a+t*b        */
+/* where s,t are scalars and p,d,a,b are 2-vectors with        */
+/* a=v1, b=v2-v1                                               */
+/***************************************************************/
+boolean intersect_line_with_segment(double px, double py, double dx, double dy,
+                                    vector3 v1, vector3 v2, double *s)
+{
+  double ax   = v1.x,       ay  = v1.y; 
+  double bx   = v2.x-v1.x,  by  = v2.y-v1.y;
+  double M00  = dx,         M10 = dy;
+  double M01  = -1.0*bx,    M11 = -1.0*by;
+  double RHSx = ax - px,    RHSy = ay - py;
+  double DetM = M00*M11 - M01*M10;
+  double L2 = bx*bx + by*by; // squared length of edge
+  if ( fabs(DetM) < 1.0e-10*L2 ) // d zero or nearly parallel to edge-->no intersection
+   return 0;
+
+  double t = (M00*RHSy-M10*RHSx)/DetM;
+  if (s) *s = (M11*RHSx-M01*RHSy)/DetM;
+
+  if (t<0.0 || t>1.0) // intersection of lines does not lie between vertices
+   return 0;
+
+  return 1;
+}
+
+// like the previous routine, but only count intersections if s>=0
+boolean intersect_ray_with_segment(double px, double py, double dx, double dy,
+                                   vector3 v1, vector3 v2, double *s)
+{
+  double ss;
+  int status=intersect_line_with_segment(px,py,dx,dy,v1,v2,&ss);
+  if (status==0 || ss<0.0)
+   return 0;
+  if (s) *s=ss;
+  return 1;
+}
+
+/***************************************************************/
+/* 2D point-in-polygon test: return 1 if the point lies within */
+/* the polygon with the given vertices, 0 otherwise.           */
+// method: cast a plumb line in the negative y direction from  */
+/* p to infinity and count the number of edges intersected;    */
+/* point lies in polygon iff this is number is odd.            */
+/***************************************************************/
+boolean point_in_polygon(double px, double py, vector3 *vertices, int num_vertices)
+{
+  double dx=0.0, dy=-1.0;
+  int num_side_intersections=0;
+  int nv;
+  for(nv=0; nv<num_vertices; nv++)
+   num_side_intersections
+    +=intersect_ray_with_segment(px, py, dx, dy,
+                                 vertices[nv], vertices[(nv+1)%num_vertices],
+                                 0);
+  return num_side_intersections%2;
+}
+
+void add_to_list(double s, double *slist, int length)
+{
+  switch(length) 
+   { case 0: 
+       slist[0] = s;
+       break;
+     case 1: 
+       if (s>=slist[0])
+        slist[1]=s;
+       else
+        { slist[1]=slist[0]; slist[0]=s; }
+       break;
+     default:
+       if (s<slist[0])
+        { slist[1]=slist[0]; slist[0]=s; }
+       else if (s<slist[1])
+        slist[1]=s;
+       break;
+   }
+     
+}
+
+/***************************************************************/
+/* return 1 or 0 if xc lies inside or outside the prism        */
+/***************************************************************/
+boolean point_in_prism(prism *prsm, vector3 xc)
+{ 
+  vector3 *vertices = prsm->vertices.items;
+  int num_vertices  = prsm->vertices.num_items;
+  double height     = prsm->height;
+  vector3 xp        = prism_coordinate_c2p(prsm, xc);
+  if ( xp.z<0.0 || xp.z>prsm->height)
+   return 0;
+  return point_in_polygon(xp.x, xp.y, vertices, num_vertices);
+}
+
+/***************************************************************/
+/* insert s into slist, which is sorted in increasing order    */
+/* has maximum length 2                                        */
+/***************************************************************/
+int insert_s_in_list(double s, double slist[2], int length)
+{
+  if (length==0) 
+   { slist[0]=s; 
+     return 1;
+   }
+  if ( s < slist[0] )
+   { slist[1]=slist[0];
+     slist[0]=s;
+   }
+  else if (length==1 || s < slist[1] )
+   slist[1]=s;
+
+  return 2;
+}
+
+/***************************************************************/
+/* compute all values of s at which the line p+s*d intersects  */
+/* a prism face, retaining the first two entries in a list     */
+/* of s values sorted in increasing order m                    */
+/***************************************************************/
+int intersect_line_with_prism(prism *prsm, vector3 p, vector3 d, double slist[2])
+{
+  int num_vertices  = prsm->vertices.num_items;
+  vector3 *vertices = prsm->vertices.items;
+  double height     = prsm->height;
+  
+  // get coords of p (line origin) and components of d (line slope)
+  // in prism coordinate system
+  vector3 p_prsm  = prism_coordinate_c2p(prsm,p);
+  vector3 d_prsm  = prism_vector_c2p(prsm,d);
+
+  // use length of first edge as a length scale for judging
+  // lengths to be small or large
+  double length_scale = vector3_norm(vector3_minus(vertices[1], vertices[0]));
+
+  // identify intersections with prism side faces 
+  int num_intersections=0;
+  int nv;
+  for(nv=0; nv<num_vertices; nv++)
+   { 
+     int nvp1 = nv+1; if (nvp1==num_vertices) nvp1=0;
+
+     // first solve the in-plane portion of the problem: determine the
+     // intersection of the XY-plane projection of the line with the
+     // polygon edge between vertices (nv,nv+1).
+     double s;
+     if (!intersect_line_with_segment(p_prsm.x, p_prsm.y, d_prsm.x, d_prsm.y,
+                                      vertices[nv], vertices[nvp1], &s)
+        ) continue;
+
+     // OK, we know the XY-plane projection of the line intersects the polygon edge;
+     // now go back to 3D, ask for the z-coordinate at which the line intersects
+     // the z-axis extrusion of the edge, and determine whether this point lies 
+     // inside or outside the height of the prism.
+     double z_intersect = p_prsm.z + s*d_prsm.z;
+     double AbsTol = 1.0e-7*length_scale;
+     double z_min  = 0.0    - AbsTol;
+     double z_max  = height + AbsTol;
+     if ( (z_intersect<z_min) || (z_intersect>z_max) )
+      continue;
+
+     num_intersections=insert_s_in_list(s, slist, num_intersections);
+   }
+
+  // identify intersections with prism ceiling and floor faces
+  double dz = d_prsm.z;
+  int LowerUpper;
+  if ( fabs(dz) > 1.0e-7*vector3_norm(d_prsm))
+   for(LowerUpper=0; LowerUpper<2; LowerUpper++)
+    { double z0    = LowerUpper ? height : 0.0;
+      double s     = (z0 - p_prsm.z)/dz;
+      if (!point_in_polygon(p_prsm.x + s*d_prsm.x, p_prsm.y+s*d_prsm.y, vertices, num_vertices))
+       continue;
+      num_intersections=insert_s_in_list(s,slist,num_intersections);
+    }
+  return num_intersections;
+}
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+double intersect_line_segment_with_prism(prism *prsm, vector3 p, vector3 d, double a, double b)
+{
+  double slist[2]={0.0, 0.0};
+  if (2!=intersect_line_with_prism(prsm, p, d, slist))
+   return 0.0;
+  double ds = fmin(slist[1],b) - fmax(slist[0],a);
+  return ds > 0.0 ? ds : 0.0;
+}
+
+/***************************************************************/
+/* compute the minimum-length vector from p to the plane       */
+/* containing o (origin) and spanned by basis vectors v1,v2    */
+/* algorithm: solve the 3x3 system p-s*v3 = o + t*v1 + u*v2    */
+/* where v3 = v1 x v2 and s,t,u are unknowns                   */
+/***************************************************************/
+vector3 normal_to_plane(vector3 o, vector3 v1, vector3 v2, vector3 p)
+{
+  vector3 RHS = vector3_minus(p,o);
+
+  // handle the degenerate-to-2D case
+  if ( (fabs(v1.z) + fabs(v2.z)) < 2.0e-7 && fabs(RHS.z)<1.0e-7 )
+   { vector3 zhat={0,0,1};
+     vector3 v3 = vector3_cross(zhat, v1);
+     double M00  = v1.x;     double M10  = v3.x;
+     double M01  = v1.y;     double M11  = v3.y;
+     double DetM = M00*M11 - M01*M10;
+     if ( fabs(DetM) < 1.0e-10 )
+      return vector3_scale(0.0, v3);
+     // double t = (M00*RHSy-M10*RHSx)/DetM;
+     double s= (M11*RHS.x-M01*RHS.y)/DetM;
+     return vector3_scale(-1.0*s,v3);
+   }
+
+  vector3 v3 = vector3_cross(v1, v2);
+  matrix3x3 M;
+  M.c0 = v1;
+  M.c1 = v2;
+  M.c2 = vector3_scale(-1.0, v3);
+  vector3 tus = matrix3x3_vector3_mult(matrix3x3_inverse(M),RHS);
+  return vector3_scale(-1.0*tus.z, v3);
+}
+
+/***************************************************************/
+/* find the face of the prism for which the normal distance    */
+/* from x to the plane of that face is the shortest, then      */
+/* return the normal vector to that plane.                     */
+/***************************************************************/
+vector3 normal_to_prism(prism *prsm, vector3 xc)
+{
+  vector3 centroid  = prsm->centroid;
+  double height     = prsm->height;
+  vector3 *vertices = prsm->vertices.items;
+  int num_vertices  = prsm->vertices.num_items;
+
+  vector3 xp   = prism_coordinate_c2p(prsm, xc);
+  vector3 axis = {0,0,0}; axis.z=height;
+ 
+  vector3 retval;
+  double min_distance;
+
+  // side walls
+  int nv;
+  for(nv=0; nv<num_vertices; nv++)
+   { int nvp1 = ( nv==(num_vertices-1) ? 0 : nv+1 );
+     vector3 v1 = vector3_minus(vertices[nvp1],vertices[nv]);
+     vector3 v2 = axis;
+     vector3 v3 = normal_to_plane(vertices[nv], v1, v2, xp);
+     double distance = vector3_norm(v3);
+     if (nv==0 || distance < min_distance)
+      { min_distance = distance;
+        retval = v3;
+      }
+   }
+
+  // floor and ceiling
+  int UpperLower;
+  for(UpperLower=0; UpperLower<2; UpperLower++)
+   { vector3 zhat={0,0,1.0};
+     vector3 v1 = vector3_minus(vertices[1],vertices[0]);
+     vector3 v2 = vector3_cross(zhat,v1);
+     vector3 o  = {0,0,0};
+     if (UpperLower) o.z = height;
+     vector3 v3 = normal_to_plane(o, v1, v2, xp);
+     double distance = vector3_norm(v3);
+     if (distance < min_distance)
+      { min_distance = distance;
+        retval = v3;
+      }
+   }
+
+  return prism_vector_p2c(prsm, retval);
+}
+
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+void display_prism_info(int indentby, prism *prsm)
+{
+  vector3 *vertices = prsm->vertices.items;
+  int num_vertices  = prsm->vertices.num_items;
+  double height     = prsm->height;
+  vector3 z0        = {0.0, 0.0, 1.0};
+  vector3 axis      = matrix3x3_vector3_mult(prsm->m_c2p,z0);
+
+  printf("%*s     height %g, axis (%g,%g,%g), %i vertices:\n", indentby, "",
+          height, axis.x, axis.y, axis.z, num_vertices);
+  matrix3x3 m_p2c = matrix3x3_inverse(prsm->m_c2p);
+  int nv;
+  for(nv=0; nv<num_vertices; nv++)
+   { vector3 v = matrix3x3_vector3_mult(m_p2c, vertices[nv]);
+     printf("%*s     {%e,%e,%e}\n",indentby,"",v.x,v.y,v.z);
+   };
+
+}
+
+/***************************************************************/
+// like vector3_equal but tolerant of small floating-point discrepancies
+/***************************************************************/
+int vector3_nearly_equal(vector3 v1, vector3 v2, double tolerance)
+{
+  return (vector3_norm( vector3_minus(v1,v2) ) <= tolerance*vector3_norm(v1));
+}
+
+/***************************************************************/
+/* return the unit normal vector to the triangle with the given*/
+/* vertices                                                    */
+/***************************************************************/
+vector3 triangle_normal(vector3 v1, vector3 v2, vector3 v3)
+{
+  vector3 nv=vector3_cross( vector3_minus(v2,v1), vector3_minus(v3,v1) );
+  double nvnorm=vector3_norm(nv);
+ // if (area) *area += 0.5*nvnorm;
+  return unit_vector3(nv);
+}
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+GEOMETRIC_OBJECT make_prism(MATERIAL_TYPE material,
+			    const vector3 *vertices, int num_vertices,
+			    double height, vector3 axis)
+{
+  CHECK(num_vertices>=3, "fewer than 3 vertices in make_prism");
+   
+  // compute centroid of vertices
+  vector3 centroid = {0.0, 0.0, 0.0};
+  int nv;
+  for(nv=0; nv<num_vertices; nv++)
+   centroid = vector3_plus(centroid, vertices[nv]);
+  centroid = vector3_scale(1.0/((double)num_vertices), centroid);
+
+  // make sure all vertices lie in a plane normal to the given axis
+  vector3 zhat = unit_vector3(axis);
+  double tolerance=1.0e-6;
+  for(nv=0; nv<num_vertices; nv++)
+   { int nvp1 = (nv+1) % num_vertices;
+     vector3 zhatp = triangle_normal(centroid,vertices[nv],vertices[nvp1]);
+     boolean axis_normal 
+      = (    vector3_nearly_equal(zhat, zhatp, tolerance)
+          || vector3_nearly_equal(zhat, vector3_scale(-1.0,zhatp), tolerance)
+        );
+     CHECK(axis_normal, "axis not normal to vertex plane in make_prism");
+   }
+
+  // compute rotation matrix that operates on a vector of cartesian coordinates
+  // to yield the coordinates of the same point in the prism coordinate system.
+  // the prism coordinate system is a right-handed coordinate system
+  // in which the prism lies in the xy plane (extrusion axis is the positive z-axis)
+  // with centroid at the origin and the first edge lying on the positive x-axis.
+  // note: the prism *centroid* is the center of mass of the planar vertex polygon,
+  //       i.e. it is a point lying on the bottom surface of the prism.
+  //       This is the origin of coordinates in the prism system.
+  //       The *center* of the geometric object is the center of mass of the 
+  //       3D prism. So center = centroid + 0.5*zHat.
+  vector3 xhat    = unit_vector3(vector3_minus(vertices[1],vertices[0]));
+  vector3 yhat    = unit_vector3(vector3_cross(zhat,xhat));
+  matrix3x3 m_p2c = {xhat, yhat, zhat };
+  matrix3x3 m_c2p = matrix3x3_inverse(m_p2c);
+
+  prism *prsm = MALLOC1(prism);
+  CHECK(prsm, "out of memory");
+  prsm->centroid  = centroid;
+  prsm->height    = height;
+  prsm->m_p2c     = m_p2c;
+  prsm->m_c2p     = m_c2p;
+
+  // note that the vertices stored in the prism_data structure
+  // are the vertices in the *prism* coordinate system, which means
+  // their z-coordinates are all zero and in principle need not be stored
+  prsm->vertices.num_items = num_vertices;
+  prsm->vertices.items     = (vector3 *)malloc(num_vertices*sizeof(vector3));
+  for(nv=0; nv<num_vertices; nv++)
+   prsm->vertices.items[nv] = prism_coordinate_c2p(prsm,vertices[nv]);
+ 
+  // note the center and centroid are different!
+  vector3 center = vector3_plus(centroid, vector3_scale(0.5*height,zhat) );
+  geometric_object o=make_geometric_object(material, center);
+  o.which_subclass=GEOM PRISM;
+  o.subclass.prism_data = prsm;
+
+  return o;
+}
diff --git a/utils/geom.scm b/utils/geom.scm
index ee0325a..42aa40e 100644
--- a/utils/geom.scm
+++ b/utils/geom.scm
@@ -85,6 +85,17 @@
 	(object-property-value object 'e2)
 	(object-property-value object 'e3))))))
 
+(define identity_matrix (matrix3x3 (vector3 1 0 0) 
+                                   (vector3 0 1 0) 
+                                   (vector3 0 0 1)))
+
+(define-class prism geometric-object
+  (define-property vertices '() (make-list-type 'vector3))
+  (define-property centroid (vector3 0 0 0) 'vector3)
+  (define-property height 0 'number)
+  (define-property m_c2p identity_matrix 'matrix3x3)
+  (define-property m_p2c identity_matrix 'matrix3x3))
+
 (define-class ellipsoid block
   (define-derived-property inverse-semi-axes 'vector3
     (lambda (object)
diff --git a/utils/test-prism.c b/utils/test-prism.c
new file mode 100644
index 0000000..97c6896
--- /dev/null
+++ b/utils/test-prism.c
@@ -0,0 +1,415 @@
+/* libctl: flexible Guile-based control files for scientific software
+ * Copyright (C) 1998-2014 Massachusetts Institute of Technology and Steven G. Johnson
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the
+ * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ * Boston, MA  02111-1307, USA.
+ *
+ * Steven G. Johnson can be contacted at stevenj@alum.mit.edu.
+ */
+
+/************************************************************************/
+/* test-prism.c: unit test for prisms in libctlgeom                     */
+/* homer reid 5/2018                                                    */
+/************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include <string.h>
+
+#include "ctlgeom.h"
+
+#define K_PI 3.141592653589793238462643383279502884197
+
+#define NUMPTS   10000
+#define NUMLINES 1000
+
+#define LX 0.5
+#define LY 1.0
+#define LZ 1.5
+
+// routine from geom.c that rotates the coordinate of a point
+// from the prism coordinate system to the cartesian coordinate system
+vector3 prism_coordinate_p2c(prism *prsm, vector3 vp);
+
+static vector3 make_vector3(double x, double y, double z)
+{
+  vector3 v;
+  v.x = x; v.y = y; v.z = z;
+  return v;
+}
+
+/************************************************************************/
+/* return a uniform random number in [a,b] */
+/************************************************************************/
+static double urand(double a, double b)
+{ return a + (b-a)*(rand()/((double)RAND_MAX)); }
+
+static double drand() 
+{ return urand(0.0,1.0); }
+
+/************************************************************************/
+/* random point uniformly distributed over a parallelepiped             */
+/************************************************************************/
+vector3 random_point(vector3 min_corner, vector3 max_corner)
+{ return make_vector3( urand(min_corner.x, max_corner.x),
+                       urand(min_corner.y, max_corner.y),
+                       urand(min_corner.z, max_corner.z)
+                     );
+}
+
+/************************************************************************/
+/* random unit vector uniformly distributed over a sphere               */
+/************************************************************************/
+vector3 random_unit_vector3(void)
+{
+  double z   = urand(-1.0,1.0);
+  double rho = sqrt(1 - z*z);
+  double phi = urand(0.0, 2.0*K_PI);
+  return make_vector3( rho*cos(phi), rho*sin(phi), z );
+}
+
+/***************************************************************/
+/* write prism vertices and edges to text file.                */
+/* after running this routine to produce a file named MyFile,  */
+/* the prism may be plotted in gnuplot like this:              */
+/* gnuplot> splot 'MyFile' u 1:2:3 w lp pt 7 ps 1              */
+/***************************************************************/
+void prism2gnuplot(prism *prsm, char *filename)
+{ 
+  vector3 *vertices = prsm->vertices.items;
+  int num_vertices  = prsm->vertices.num_items;
+  double height     = prsm->height;
+
+  FILE *f=fopen(filename,"w");
+  int nv;
+  for(nv=0; nv<num_vertices; nv++)
+   {  
+     vector3 vap = vertices[nv]; vap.z = 0.0;
+     vector3 vbp = vertices[nv]; vbp.z = height;
+     vector3 vcp = vertices[(nv+1)%num_vertices]; vcp.z = height;
+     vector3 vdp = vertices[(nv+1)%num_vertices]; vdp.z = 0.0;
+     vector3 vac = prism_coordinate_p2c(prsm,vap);
+     vector3 vbc = prism_coordinate_p2c(prsm,vbp);
+     vector3 vcc = prism_coordinate_p2c(prsm,vcp);
+     vector3 vdc = prism_coordinate_p2c(prsm,vdp);
+
+     fprintf(f,"%e %e %e %e %e %e \n",vac.x,vac.y,vac.z,vap.x,vap.y,vap.z);
+     fprintf(f,"%e %e %e %e %e %e \n",vbc.x,vbc.y,vbc.z,vbp.x,vbp.y,vbp.z);
+     fprintf(f,"%e %e %e %e %e %e \n",vcc.x,vcc.y,vcc.z,vcp.x,vcp.y,vcp.z);
+     fprintf(f,"%e %e %e %e %e %e \n",vdc.x,vdc.y,vdc.z,vdp.x,vdp.y,vdp.z);
+     fprintf(f,"%e %e %e %e %e %e \n",vac.x,vac.y,vac.z,vap.x,vap.y,vap.z);
+     fprintf(f,"\n\n");
+   }
+  fclose(f);
+}
+
+/***************************************************************/
+/* write prism vertices and edges to GMSH geometry (.geo) file */
+/***************************************************************/
+void prism2gmsh(prism *prsm, char *filename)
+{ 
+  vector3 *vertices = prsm->vertices.items;
+  int num_vertices  = prsm->vertices.num_items;
+  double height     = prsm->height;
+  vector3 zhat      = prsm->m_p2c.c2;
+  vector3 axis      = vector3_scale(height, zhat);
+
+  FILE *f=fopen(filename,"w");
+  int nv;
+  for(nv=0; nv<num_vertices; nv++)
+   { vector3 vp=vertices[nv];
+     vector3 vc=prism_coordinate_p2c(prsm, vp);
+     fprintf(f,"Point(%i)={%e, %e, %e};\n",nv,vc.x,vc.y,vc.z);
+   }
+  for(nv=0; nv<num_vertices; nv++)
+   fprintf(f,"Line(%i)={%i, %i};\n",nv,nv,(nv+1)%num_vertices);
+  fprintf(f,"Line Loop(0)={0");
+  for(nv=1; nv<num_vertices; nv++)
+   fprintf(f,",%i",nv);
+  fprintf(f,"};\n");
+  fprintf(f,"Plane Surface(0)={0};\n");
+  fprintf(f,"Extrude { %e,%e,%e } { Surface{0}; }\n",height*zhat.x,height*zhat.y,height*zhat.z);
+  fclose(f); 
+  
+}
+
+/* "standardize" a vector for vector comparisons up to normalization and sign flip */
+double sgn(double x) { return x>=0.0 ? 1.0 : -1.0; }
+
+vector3 standardize(vector3 v)
+{ vector3 sv=unit_vector3(v);
+  double sign = (sv.z!=0.0 ? sgn(sv.z) : sv.y!=0.0 ? sgn(sv.y) : sgn(sv.x));
+  return vector3_scale(sign,sv);
+}
+
+/************************************************************************/
+/* first unit test: check inclusion of randomly-generated points        */
+/************************************************************************/
+int test_point_inclusion(geometric_object the_block, geometric_object the_prism,
+                         int num_tests, int write_log)
+{
+  vector3 size = the_block.subclass.block_data->size;
+  vector3 min_corner = vector3_scale(-1.0, size);
+  vector3 max_corner = vector3_scale(+1.0, size);
+  FILE *f = write_log ? fopen("/tmp/test-prism.points","w") : 0;
+  int num_failed=0;
+  int n;
+  for(n=0; n<num_tests; n++)
+   { 
+     vector3 p = random_point(min_corner, max_corner);
+     boolean in_block=point_in_objectp(p,the_block);
+     boolean in_prism=point_in_objectp(p,the_prism);
+     if (in_block!=in_prism)
+      num_failed++;
+     if (f) fprintf(f,"%i %i %e %e %e \n",in_block,in_prism,p.x,p.y,p.z);
+   }
+  if (f) fclose(f);
+  
+  printf("point inclusion: %i/%i points failed\n",num_failed,num_tests);
+  return num_failed;
+}
+
+/************************************************************************/
+/* second unit test: check calculation of normals to objects            */
+/************************************************************************/
+int test_normal_to_object(geometric_object the_block, geometric_object the_prism,
+                          int num_tests, int write_log)
+{
+  vector3 size = the_block.subclass.block_data->size;
+  vector3 min_corner = vector3_scale(-1.0, size);
+  vector3 max_corner = vector3_scale(+1.0, size);
+  FILE *f = write_log ? fopen("/tmp/test-prism.normals","w") : 0;
+
+  int num_failed=0;
+  int n;
+  double tolerance=1.0e-6;
+  for(n=0; n<num_tests; n++)
+   { 
+     // randomly generated base point within enlarged bounding box
+     vector3 p = random_point(min_corner, max_corner);
+
+     vector3 nhat_block=standardize(normal_to_object(p, the_block));
+     vector3 nhat_prism=standardize(normal_to_object(p, the_prism));
+     if (!vector3_nearly_equal(nhat_block, nhat_prism, tolerance))
+      num_failed++;
+
+     if (f)
+      fprintf(f,"%e %e %e %e %e %e %e %e %e %i\n\n\n",p.x,p.y,p.z,
+                 nhat_block.x,nhat_block.y,nhat_block.z,
+                 nhat_prism.x,nhat_prism.y,nhat_prism.z,
+                 vector3_nearly_equal(nhat_block,nhat_prism,tolerance));
+   }
+  if (f) fclose(f);
+  
+  printf("%i/%i normals failed\n",num_failed,num_tests);
+  return num_failed;
+}
+
+/************************************************************************/
+/* third unit test: check-line segment intersections                   */
+/************************************************************************/
+int test_line_segment_intersection(geometric_object the_block, geometric_object the_prism,
+                                   int num_tests, int write_log)
+{
+  vector3 size = the_block.subclass.block_data->size;
+  vector3 min_corner = vector3_scale(-1.0, size);
+  vector3 max_corner = vector3_scale(+1.0, size);
+  FILE *f = write_log ? fopen("/tmp/test-prism.segments","w") : 0;
+
+  int num_failed=0;
+  int n;
+  for(n=0; n<num_tests; n++)
+   { 
+     // randomly generated base point within enlarged bounding box
+     vector3 p = random_point(min_corner, max_corner);
+     vector3 d = random_unit_vector3();
+     double a  = urand(0.0,1.0);
+     double b  = urand(0.0,1.0);
+
+     double sblock = intersect_line_segment_with_object(p, d, the_block, a, b);
+     double sprism = intersect_line_segment_with_object(p, d, the_prism, a, b);
+     if ( fabs(sblock-sprism) > 1.0e-6*fmax(fabs(sblock),fabs(sprism)) )
+      num_failed++;
+
+     if (f)
+      fprintf(f," %e %e %s\n",sblock,sprism,fabs(sblock-sprism) > 1.0e-6*fmax(fabs(sblock),fabs(sprism)) ? "fail" : "success");
+   }
+  if (f) fclose(f);
+  
+  printf("%i/%i segments failed\n",num_failed,num_tests);
+  return num_failed;
+}
+
+/***************************************************************/
+/* unit tests: create the same parallelepiped two ways (as a   */
+/* block and as a prism) and verify that geometric primitives  */
+/* give identical results                                      */
+/***************************************************************/
+int run_unit_tests()
+{
+  void* m = NULL;
+  vector3 c = { 0, 0, 0 };
+  vector3 xhat = make_vector3(1,0,0);
+  vector3 yhat = make_vector3(0,1,0);
+  vector3 zhat = make_vector3(0,0,1);
+  vector3 size = make_vector3(LX,LY,LZ);
+  geometric_object the_block = make_block(m, c, xhat, yhat, zhat, size);
+
+  vector3 v[4]; 
+  v[0].x=-0.5*LX; v[0].y=-0.5*LY; v[0].z=-0.5*LZ;
+  v[1].x=+0.5*LX; v[1].y=-0.5*LY; v[1].z=-0.5*LZ;
+  v[2].x=+0.5*LX; v[2].y=+0.5*LY; v[2].z=-0.5*LZ;
+  v[3].x=-0.5*LX; v[3].y=+0.5*LY; v[3].z=-0.5*LZ;
+  geometric_object the_prism=make_prism(m, v, 4, LZ, zhat);
+
+  int write_log=0;
+
+  if (write_log)
+   prism2gnuplot(the_prism.subclass.prism_data, "/tmp/test-prism.prism");
+
+  int num_failed_1 = test_point_inclusion(the_block, the_prism, NUMPTS, write_log);
+  int num_failed_2 = test_normal_to_object(the_block, the_prism, NUMLINES, write_log);
+  int num_failed_3 = test_line_segment_intersection(the_block, the_prism, NUMLINES, write_log);
+
+  return num_failed_1 + num_failed_2 + num_failed_3;
+}
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+void print_usage(char *msg, int print_usage)
+{ 
+  if (!msg)
+   fprintf(stderr,"%s\n",msg);
+  if (print_usage)
+   { printf("usage: \n");
+     printf(" --vertexfile MyVertices\n");
+     printf(" --height     height\n");
+     printf(" --axis       x y z\n");
+     printf("\n");
+     printf(" --point      x y z\n");
+     printf(" --dir        x y z\n");
+     printf(" --a          a\n");
+     printf(" --b          b\n");
+   }
+  exit(1);
+}
+
+void quit(char *msg)
+{ print_usage(msg, 0); }
+
+void usage(char *msg)
+{ print_usage(msg, 1); }
+
+/************************************************************************/
+/************************************************************************/
+/************************************************************************/
+int main(int argc, char *argv[])
+{
+  srand(time(NULL));
+  geom_initialize();
+
+  if (argc<=1) // if no arguments, run unit tests
+   return run_unit_tests();
+
+  /***************************************************************/
+  /* process arguments *******************************************/
+  /***************************************************************/
+  char *vertexfile=0;
+  vector3 axis={0,0,1};
+  double height=1.5;
+  vector3 test_point={0,0,0};
+  vector3 test_dir={0,0,1};
+  double a = 0.2, b=0.3; 
+  int narg;
+  for(narg=1; narg<argc-1; narg++)
+   { if (!strcmp(argv[narg],"--vertexfile"))
+      vertexfile=argv[++narg];
+     else if (!strcmp(argv[narg],"--axis"))
+      { if (narg+3>=argc) usage("too few arguments to --axis");
+        sscanf(argv[narg+1],"%le",&(axis.x));
+        sscanf(argv[narg+2],"%le",&(axis.y));
+        sscanf(argv[narg+3],"%le",&(axis.z));
+        narg+=3;
+      }
+     else if (!strcmp(argv[narg],"--point"))
+      { if (narg+3>=argc) usage("too few arguments to --point");
+        sscanf(argv[narg+1],"%le",&(test_point.x));
+        sscanf(argv[narg+2],"%le",&(test_point.y));
+        sscanf(argv[narg+3],"%le",&(test_point.z));
+        narg+=3;
+      }
+     else if (!strcmp(argv[narg],"--dir"))
+      { if (narg+5>=argc) usage("too few arguments to --lineseg");
+        sscanf(argv[narg+1],"%le",&(test_dir.x));
+        sscanf(argv[narg+2],"%le",&(test_dir.y));
+        sscanf(argv[narg+3],"%le",&(test_dir.z));
+        narg+=3;
+      }
+     else if (!strcmp(argv[narg],"--height"))
+      sscanf(argv[++narg],"%le",&height);
+     else if (!strcmp(argv[narg],"--a"))
+      sscanf(argv[++narg],"%le",&a);
+     else if (!strcmp(argv[narg],"--b"))
+      sscanf(argv[++narg],"%le",&b);
+     else 
+      usage("unknown argument");
+   }
+  if (!vertexfile) usage("no --vertexfile specified");
+
+  /***************************************************************/
+  /* read vertices from vertex file and create prism *************/
+  /***************************************************************/
+  vector3 *vertices=0;
+  int num_vertices=0;
+  FILE *f=fopen(vertexfile,"r");
+  if (!f) usage("could not open vertexfile");
+  char Line[100];
+  int LineNum=0;
+  while( fgets(Line,100,f) )
+   { if (Line[0]=='\n' || Line[0]=='#') continue;
+      num_vertices++;
+     vector3 v;
+     if ( 3!=sscanf(Line,"%le %le %le\n",&(v.x),&(v.y),&(v.z)) )
+      { fprintf(stderr,"bad vertex on line %i of %s",num_vertices,vertexfile);
+        exit(1);
+      }
+     vertices = (vector3 *)realloc(vertices, num_vertices*sizeof(vector3));
+     vertices[num_vertices-1]=v;
+   }
+  fclose(f);
+
+  geometric_object the_prism=make_prism(NULL, vertices, num_vertices, height, axis);
+  prism *prsm=the_prism.subclass.prism_data;
+  prism2gmsh(prsm, "test-prism.pp");
+  prism2gnuplot(prsm, "test-prism.gp");
+  printf("Wrote prism description to GNUPLOT file test-prism.gp.\n");
+  printf("Wrote prism description to GMSH file test-prism.geo.\n");
+
+  /***************************************************************/
+  /* test point inclusion, normal to object, and line-segment    */
+  /* intersection with specified data                            */
+  /***************************************************************/
+  boolean in_prism=point_in_objectp(test_point,the_prism);
+  vector3 nhat=normal_to_object(test_point, the_prism);
+  double s= intersect_line_segment_with_object(test_point, test_dir, the_prism, a, b);
+  printf("point {%e,%e,%e}: \n",test_point.x,test_point.y,test_point.z);
+  printf(" %s prism\n", in_prism ? "in" : "not in");
+  printf(" normal to prism: {%e,%e,%e}\n",nhat.x,nhat.y,nhat.z);
+  printf(" intersection with line segment {%e,%e,%e} + (%e,%e)*{%e,%e,%e}: %e\n",
+           test_point.x, test_point.y, test_point.z, 
+           a,b,test_dir.x, test_dir.y, test_dir.z,s);
+}