Introduce Position/Angle types

CMakeLists.txt

@@ -386,6 +386,7 @@ add_library(libopenmc SHARED
   src/cell.cpp
   src/initialize.cpp
   src/finalize.cpp
+  src/geometry.cpp
   src/geometry_aux.cpp
   src/hdf5_interface.cpp
   src/lattice.cpp

src/cell.cpp

@@ -279,33 +279,31 @@ Cell::Cell(pugi::xml_node cell_node)
 //==============================================================================
 
 bool
-Cell::contains(const double xyz[3], const double uvw[3],
-               int32_t on_surface) const
+Cell::contains(Position r, Angle a, int32_t on_surface) const
 {
   if (simple) {
-    return contains_simple(xyz, uvw, on_surface);
+    return contains_simple(r, a, on_surface);
   } else {
-    return contains_complex(xyz, uvw, on_surface);
+    return contains_complex(r, a, on_surface);
   }
 }
 
 //==============================================================================
 
 std::pair<double, int32_t>
-Cell::distance(const double xyz[3], const double uvw[3],
-               int32_t on_surface) const
+Cell::distance(Position r, Angle a, int32_t on_surface) const
 {
   double min_dist {INFTY};
   int32_t i_surf {std::numeric_limits<int32_t>::max()};
 
   for (int32_t token : rpn) {
     // Ignore this token if it corresponds to an operator rather than a region.
-    if (token >= OP_UNION) {continue;}
+    if (token >= OP_UNION) continue;
 
     // Calculate the distance to this surface.
     // Note the off-by-one indexing
     bool coincident {token == on_surface};
-    double d {surfaces_c[abs(token)-1]->distance(xyz, uvw, coincident)};
+    double d {surfaces_c[abs(token)-1]->distance(r, a, coincident)};
 
     // Check if this distance is the new minimum.
     if (d < min_dist) {
@@ -356,8 +354,7 @@ Cell::to_hdf5(hid_t cell_group) const
 //==============================================================================
 
 bool
-Cell::contains_simple(const double xyz[3], const double uvw[3],
-                      int32_t on_surface) const
+Cell::contains_simple(Position r, Angle a, int32_t on_surface) const
 {
   for (int32_t token : rpn) {
     if (token < OP_UNION) {
@@ -370,7 +367,7 @@ Cell::contains_simple(const double xyz[3], const double uvw[3],
         return false;
       } else {
         // Note the off-by-one indexing
-        bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);
+        bool sense = surfaces_c[abs(token)-1]->sense(r, a);
         if (sense != (token > 0)) {return false;}
       }
     }
@@ -381,8 +378,7 @@ Cell::contains_simple(const double xyz[3], const double uvw[3],
 //==============================================================================
 
 bool
-Cell::contains_complex(const double xyz[3], const double uvw[3],
-                       int32_t on_surface) const
+Cell::contains_complex(Position r, Angle a, int32_t on_surface) const
 {
   // Make a stack of booleans.  We don't know how big it needs to be, but we do
   // know that rpn.size() is an upper-bound.
@@ -413,7 +409,7 @@ Cell::contains_complex(const double xyz[3], const double uvw[3],
         stack[i_stack] = false;
       } else {
         // Note the off-by-one indexing
-        bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);;
+        bool sense = surfaces_c[abs(token)-1]->sense(r, a);;
         stack[i_stack] = (sense == (token > 0));
       }
     }
@@ -494,12 +490,18 @@ extern "C" {
   bool cell_simple(Cell *c) {return c->simple;}
 
   bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface)
-  {return c->contains(xyz, uvw, on_surface);}
+  {
+    Position r {xyz};
+    Angle a {uvw};
+    return c->contains(r, a, on_surface);
+  }
 
   void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface,
                      double *min_dist, int32_t *i_surf)
   {
-    std::pair<double, int32_t> out = c->distance(xyz, uvw, on_surface);
+    Position r {xyz};
+    Angle a {uvw};
+    std::pair<double, int32_t> out = c->distance(r, a, on_surface);
     *min_dist = out.first;
     *i_surf = out.second;
   }

src/cell.h

@@ -8,6 +8,7 @@
 
 #include "hdf5.h"
 #include "pugixml.hpp"
+#include "geometry.h"
 
 
 namespace openmc {
@@ -97,21 +98,19 @@ public:
   //!   known to be on.  This index takes precedence over surface sense
   //!   calculations.
   bool
-  contains(const double xyz[3], const double uvw[3], int32_t on_surface) const;
+  contains(Position r, Angle a, int32_t on_surface) const;
 
   //! Find the oncoming boundary of this cell.
   std::pair<double, int32_t>
-  distance(const double xyz[3], const double uvw[3], int32_t on_surface) const;
+  distance(Position r, Angle a, int32_t on_surface) const;
 
   //! \brief Write cell information to an HDF5 group.
   //! @param group_id An HDF5 group id.
   void to_hdf5(hid_t group_id) const;
 
 protected:
-  bool contains_simple(const double xyz[3], const double uvw[3],
-                       int32_t on_surface) const;
-  bool contains_complex(const double xyz[3], const double uvw[3],
-                        int32_t on_surface) const;
+  bool contains_simple(Position r, Angle a, int32_t on_surface) const;
+  bool contains_complex(Position r, Angle a, int32_t on_surface) const;
 };
 
 } // namespace openmc

src/geometry.cpp

@@ -0,0 +1,59 @@
+#include "geometry.h"
+
+namespace openmc {
+
+Position&
+Position::operator+=(Position other)
+{
+  x += other.x;
+  y += other.y;
+  z += other.z;
+  return *this;
+}
+
+Position&
+Position::operator+=(double v)
+{
+  x += v;
+  y += v;
+  z += v;
+  return *this;
+}
+
+Position&
+Position::operator-=(Position other)
+{
+  x -= other.x;
+  y -= other.y;
+  z -= other.z;
+  return *this;
+}
+
+Position&
+Position::operator-=(double v)
+{
+  x -= v;
+  y -= v;
+  z -= v;
+  return *this;
+}
+
+Position&
+Position::operator*=(Position other)
+{
+  x *= other.x;
+  y *= other.y;
+  z *= other.z;
+  return *this;
+}
+
+Position&
+Position::operator*=(double v)
+{
+  x *= v;
+  y *= v;
+  z *= v;
+  return *this;
+}
+
+} // namespace openmc

src/geometry.h

@@ -1,6 +1,68 @@
 #ifndef GEOMETRY_H
 #define GEOMETRY_H
 
+namespace openmc {
+
 extern "C" int openmc_root_universe;
 
+struct Position {
+  double x = 0.;
+  double y = 0.;
+  double z = 0.;
+
+  Position() = default;
+  Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { };
+  Position(double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
+
+  Position& operator+=(Position);
+  Position& operator+=(double);
+  Position& operator-=(Position);
+  Position& operator-=(double);
+  Position& operator*=(Position);
+  Position& operator*=(double);
+  const double& operator[](int i) const {
+    switch (i) {
+      case 0: return x;
+      case 1: return y;
+      case 2: return z;
+    }
+  }
+  double& operator[](int i) {
+    switch (i) {
+      case 0: return x;
+      case 1: return y;
+      case 2: return z;
+    }
+  }
+
+  inline double dot(Position other) {
+    return x*other.x + y*other.y + z*other.z;
+  }
+};
+
+inline Position operator+(Position a, Position b) { return a += b; }
+inline Position operator+(Position a, double b)   { return a += b; }
+inline Position operator+(double a, Position b)   { return b += a; }
+
+inline Position operator-(Position a, Position b) { return a -= b; }
+inline Position operator-(Position a, double b)   { return a -= b; }
+inline Position operator-(double a, Position b)   { return b -= a; }
+
+inline Position operator*(Position a, Position b) { return a *= b; }
+inline Position operator*(Position a, double b)   { return a *= b; }
+inline Position operator*(double a, Position b)   { return b *= a; }
+
+
+struct Angle : Position {
+  double& u() { return x; }
+  double& v() { return y; }
+  double& w() { return z; }
+  Angle() = default;
+  Angle(double u, double v, double w) : Position{u, v, w} { };
+  Angle(double uvw[]) : Position{uvw} { };
+  Angle(Position r) : Position{r} { };
+};
+
+} // namespace openmc
+
 #endif // GEOMETRY_H

src/surface.cpp

@@ -3,6 +3,7 @@
 #include <array>
 #include <cmath>
 #include <sstream>
+#include <utility>
 
 #include "error.h"
 #include "hdf5_interface.h"
@@ -177,38 +178,33 @@ Surface::Surface(pugi::xml_node surf_node)
 }
 
 bool
-Surface::sense(const double xyz[3], const double uvw[3]) const
+Surface::sense(Position r, Angle a) const
 {
   // Evaluate the surface equation at the particle's coordinates to determine
   // which side the particle is on.
-  const double f = evaluate(xyz);
+  const double f = evaluate(r);
 
   // Check which side of surface the point is on.
   if (std::abs(f) < FP_COINCIDENT) {
     // Particle may be coincident with this surface. To determine the sense, we
     // look at the direction of the particle relative to the surface normal (by
     // default in the positive direction) via their dot product.
-    double norm[3];
-    normal(xyz, norm);
-    return uvw[0] * norm[0] + uvw[1] * norm[1] + uvw[2] * norm[2] > 0.0;
+    return a.dot(normal(r)) > 0.0;
   }
   return f > 0.0;
 }
 
-void
-Surface::reflect(const double xyz[3], double uvw[3]) const
+Angle
+Surface::reflect(Position r, Angle a) const
 {
   // Determine projection of direction onto normal and squared magnitude of
   // normal.
-  double norm[3];
-  normal(xyz, norm);
-  const double projection = norm[0]*uvw[0] + norm[1]*uvw[1] + norm[2]*uvw[2];
-  const double magnitude = norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2];
+  Angle n = normal(r);
+  const double projection = n.dot(a);
+  const double magnitude = n.dot(n);
 
   // Reflect direction according to normal.
-  uvw[0] -= 2.0 * projection / magnitude * norm[0];
-  uvw[1] -= 2.0 * projection / magnitude * norm[1];
-  uvw[2] -= 2.0 * projection / magnitude * norm[2];
+  return a -= (2.0 * projection / magnitude) * n;
 }
 
 void
@@ -261,33 +257,15 @@ PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node)
 
 // The template parameter indicates the axis normal to the plane.
 template<int i> double
-axis_aligned_plane_evaluate(const double xyz[3], double offset)
+axis_aligned_plane_distance(Position r, Angle a, bool coincident, double offset)
 {
-  return xyz[i] - offset;
-}
-
-// The template parameter indicates the axis normal to the plane.
-template<int i> double
-axis_aligned_plane_distance(const double xyz[3], const double uvw[3],
-                            bool coincident, double offset)
-{
-  const double f = offset - xyz[i];
-  if (coincident or std::abs(f) < FP_COINCIDENT or uvw[i] == 0.0) return INFTY;
-  const double d = f / uvw[i];
+  const double f = offset - r[i];
+  if (coincident or std::abs(f) < FP_COINCIDENT or a[i] == 0.0) return INFTY;
+  const double d = f / a[i];
   if (d < 0.0) return INFTY;
   return d;
 }
 
-// The first template parameter indicates the axis normal to the plane.  The
-// other two parameters indicate the other two axes.
-template<int i1, int i2, int i3> void
-axis_aligned_plane_normal(const double xyz[3], double uvw[3])
-{
-  uvw[i1] = 1.0;
-  uvw[i2] = 0.0;
-  uvw[i3] = 0.0;
-}
-
 //==============================================================================
 // SurfaceXPlane implementation
 //==============================================================================
@@ -298,20 +276,19 @@ SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, x0);
 }
 
-inline double SurfaceXPlane::evaluate(const double xyz[3]) const
+double SurfaceXPlane::evaluate(Position r) const
 {
-  return axis_aligned_plane_evaluate<0>(xyz, x0);
+  return r.x - x0;
 }
 
-inline double SurfaceXPlane::distance(const double xyz[3], const double uvw[3],
-                                      bool coincident) const
+double SurfaceXPlane::distance(Position r, Angle a, bool coincident) const
 {
-  return axis_aligned_plane_distance<0>(xyz, uvw, coincident, x0);
+  return axis_aligned_plane_distance<0>(r, a, coincident, x0);
 }
 
-inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceXPlane::normal(Position r) const
 {
-  axis_aligned_plane_normal<0, 1, 2>(xyz, uvw);
+  return {1., 0., 0.};
 }
 
 void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
@@ -321,26 +298,24 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
-     double uvw[3]) const
+bool SurfaceXPlane::periodic_translate(const PeriodicSurface *other, Position& r,
+                                       Angle& a) const
 {
-  double other_norm[3];
-  other->normal(xyz, other_norm);
-  if (other_norm[0] == 1 and other_norm[1] == 0 and other_norm[2] == 0) {
-    xyz[0] = x0;
+  Angle other_n = other->normal(r);
+  if (other_n.u() == 1 and other_n.v() == 0 and other_n.w() == 0) {
+    r.x = x0;
     return false;
   } else {
     // Assume the partner is an YPlane (the only supported partner).  Use the
-    // evaluate function to find y0, then adjust xyz and uvw for rotational
+    // evaluate function to find y0, then adjust position/angle for rotational
     // symmetry.
-    double xyz_test[3] {0, 0, 0};
-    double y0 = -other->evaluate(xyz_test);
-    xyz[1] = xyz[0] - x0 + y0;
-    xyz[0] = x0;
+    double y0 = -other->evaluate({0., 0., 0.});
+    r.y = r.x - x0 + y0;
+    r.x = x0;
 
-    double u = uvw[0];
-    uvw[0] = -uvw[1];
-    uvw[1] = u;
+    double u = a.u();
+    a.u() = -a.v();
+    a.v() = u;
 
     return true;
   }
@@ -349,8 +324,7 @@ bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
 BoundingBox
 SurfaceXPlane::bounding_box() const
 {
-  BoundingBox out {x0, x0, -INFTY, INFTY, -INFTY, INFTY};
-  return out;
+  return {x0, x0, -INFTY, INFTY, -INFTY, INFTY};
 }
 
 //==============================================================================
@@ -363,20 +337,19 @@ SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, y0);
 }
 
-inline double SurfaceYPlane::evaluate(const double xyz[3]) const
+double SurfaceYPlane::evaluate(Position r) const
 {
-  return axis_aligned_plane_evaluate<1>(xyz, y0);
+  return r.y - y0;
 }
 
-inline double SurfaceYPlane::distance(const double xyz[3], const double uvw[3],
-                                      bool coincident) const
+double SurfaceYPlane::distance(Position r, Angle a, bool coincident) const
 {
-  return axis_aligned_plane_distance<1>(xyz, uvw, coincident, y0);
+  return axis_aligned_plane_distance<1>(r, a, coincident, y0);
 }
 
-inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceYPlane::normal(Position r) const
 {
-  axis_aligned_plane_normal<1, 0, 2>(xyz, uvw);
+  return {0., 1., 0.};
 }
 
 void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
@@ -386,27 +359,25 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
-                                       double uvw[3]) const
+bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, Position& r,
+                                       Angle& a) const
 {
-  double other_norm[3];
-  other->normal(xyz, other_norm);
-  if (other_norm[0] == 0 and other_norm[1] == 1 and other_norm[2] == 0) {
+  Angle other_n = other->normal(r);
+  if (other_n.u() == 0 and other_n.v() == 1 and other_n.w() == 0) {
     // The periodic partner is also aligned along y.  Just change the y coord.
-    xyz[1] = y0;
+    r.y = y0;
     return false;
   } else {
     // Assume the partner is an XPlane (the only supported partner).  Use the
-    // evaluate function to find x0, then adjust xyz and uvw for rotational
+    // evaluate function to find x0, then adjust position/angle for rotational
     // symmetry.
-    double xyz_test[3] {0, 0, 0};
-    double x0 = -other->evaluate(xyz_test);
-    xyz[0] = xyz[1] - y0 + x0;
-    xyz[1] = y0;
+    double x0 = -other->evaluate({0., 0., 0.});
+    r.x = r.y - y0 + x0;
+    r.y = y0;
 
-    double u = uvw[0];
-    uvw[0] = uvw[1];
-    uvw[1] = -u;
+    double u = a.u();
+    a.u() = a.v();
+    a.v() = -u;
 
     return true;
   }
@@ -415,8 +386,7 @@ bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
 BoundingBox
 SurfaceYPlane::bounding_box() const
 {
-  BoundingBox out {-INFTY, INFTY, y0, y0, -INFTY, INFTY};
-  return out;
+  return {-INFTY, INFTY, y0, y0, -INFTY, INFTY};
 }
 
 //==============================================================================
@@ -429,20 +399,19 @@ SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, z0);
 }
 
-inline double SurfaceZPlane::evaluate(const double xyz[3]) const
+double SurfaceZPlane::evaluate(Position r) const
 {
-  return axis_aligned_plane_evaluate<2>(xyz, z0);
+  return r.z - z0;
 }
 
-inline double SurfaceZPlane::distance(const double xyz[3], const double uvw[3],
-                                      bool coincident) const
+double SurfaceZPlane::distance(Position r, Angle a, bool coincident) const
 {
-  return axis_aligned_plane_distance<2>(xyz, uvw, coincident, z0);
+  return axis_aligned_plane_distance<2>(r, a, coincident, z0);
 }
 
-inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceZPlane::normal(Position r) const
 {
-  axis_aligned_plane_normal<2, 0, 1>(xyz, uvw);
+  return {0., 0., 1.};
 }
 
 void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
@@ -452,19 +421,18 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
-                                       double uvw[3]) const
+bool SurfaceZPlane::periodic_translate(const PeriodicSurface *other, Position& r,
+                                       Angle& a) const
 {
   // Assume the other plane is aligned along z.  Just change the z coord.
-  xyz[2] = z0;
+  r.z = z0;
   return false;
 }
 
 BoundingBox
 SurfaceZPlane::bounding_box() const
 {
-  BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, z0, z0};
-  return out;
+  return {-INFTY, INFTY, -INFTY, INFTY, z0, z0};
 }
 
 //==============================================================================
@@ -478,17 +446,16 @@ SurfacePlane::SurfacePlane(pugi::xml_node surf_node)
 }
 
 double
-SurfacePlane::evaluate(const double xyz[3]) const
+SurfacePlane::evaluate(Position r) const
 {
-  return A*xyz[0] + B*xyz[1] + C*xyz[2] - D;
+  return A*r.x + B*r.y + C*r.z - D;
 }
 
 double
-SurfacePlane::distance(const double xyz[3], const double uvw[3],
-                       bool coincident) const
+SurfacePlane::distance(Position r, Angle a, bool coincident) const
 {
-  const double f = A*xyz[0] + B*xyz[1] + C*xyz[2] - D;
-  const double projection = A*uvw[0] + B*uvw[1] + C*uvw[2];
+  const double f = A*r.x + B*r.y + C*r.z - D;
+  const double projection = A*a.u() + B*a.v() + C*a.w();
   if (coincident or std::abs(f) < FP_COINCIDENT or projection == 0.0) {
     return INFTY;
   } else {
@@ -498,12 +465,10 @@ SurfacePlane::distance(const double xyz[3], const double uvw[3],
   }
 }
 
-void
-SurfacePlane::normal(const double xyz[3], double uvw[3]) const
+Angle
+SurfacePlane::normal(Position r) const
 {
-  uvw[0] = A;
-  uvw[1] = B;
-  uvw[2] = C;
+  return {A, B, C};
 }
 
 void SurfacePlane::to_hdf5_inner(hid_t group_id) const
@@ -513,18 +478,18 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const
   write_dataset(group_id, "coefficients", coeffs);
 }
 
-bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3],
-                                      double uvw[3]) const
+bool SurfacePlane::periodic_translate(const PeriodicSurface *other, Position& r,
+                                      Angle& a) const
 {
   // This function assumes the other plane shares this plane's normal direction.
 
   // Determine the distance to intersection.
-  double d = evaluate(xyz) / (A*A + B*B + C*C);
+  double d = evaluate(r) / (A*A + B*B + C*C);
 
   // Move the particle that distance along the normal vector.
-  xyz[0] -= d * A;
-  xyz[1] -= d * B;
-  xyz[2] -= d * C;
+  r.x -= d * A;
+  r.y -= d * B;
+  r.z -= d * C;
 
   return false;
 }
@@ -532,8 +497,7 @@ bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3],
 BoundingBox
 SurfacePlane::bounding_box() const
 {
-  BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
-  return out;
+  return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY};
 }
 
 //==============================================================================
@@ -544,28 +508,28 @@ SurfacePlane::bounding_box() const
 // cylinder.  offset1 and offset2 should correspond with i1 and i2,
 // respectively.
 template<int i1, int i2> double
-axis_aligned_cylinder_evaluate(const double xyz[3], double offset1,
+axis_aligned_cylinder_evaluate(Position r, double offset1,
                                double offset2, double radius)
 {
-  const double xyz1 = xyz[i1] - offset1;
-  const double xyz2 = xyz[i2] - offset2;
-  return xyz1*xyz1 + xyz2*xyz2 - radius*radius;
+  const double r1 = r[i1] - offset1;
+  const double r2 = r[i2] - offset2;
+  return r1*r1 + r2*r2 - radius*radius;
 }
 
 // The first template parameter indicates which axis the cylinder is aligned to.
 // The other two parameters indicate the other two axes.  offset1 and offset2
 // should correspond with i2 and i3, respectively.
 template<int i1, int i2, int i3> double
-axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3],
+axis_aligned_cylinder_distance(Position r, Angle u,
      bool coincident, double offset1, double offset2, double radius)
 {
-  const double a = 1.0 - uvw[i1]*uvw[i1];  // u^2 + v^2
+  const double a = 1.0 - u[i1]*u[i1];  // u^2 + v^2
   if (a == 0.0) return INFTY;
 
-  const double xyz2 = xyz[i2] - offset1;
-  const double xyz3 = xyz[i3] - offset2;
-  const double k = xyz2 * uvw[i2] + xyz3 * uvw[i3];
-  const double c = xyz2*xyz2 + xyz3*xyz3 - radius*radius;
+  const double r2 = r[i2] - offset1;
+  const double r3 = r[i3] - offset2;
+  const double k = r2 * u[i2] + r3 * u[i3];
+  const double c = r2*r2 + r3*r3 - radius*radius;
   const double quad = k*k - a*c;
 
   if (quad < 0.0) {
@@ -601,13 +565,14 @@ axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3],
 // The first template parameter indicates which axis the cylinder is aligned to.
 // The other two parameters indicate the other two axes.  offset1 and offset2
 // should correspond with i2 and i3, respectively.
-template<int i1, int i2, int i3> void
-axis_aligned_cylinder_normal(const double xyz[3], double uvw[3], double offset1,
-                             double offset2)
+template<int i1, int i2, int i3> Angle
+axis_aligned_cylinder_normal(Position r, double offset1, double offset2)
 {
-  uvw[i2] = 2.0 * (xyz[i2] - offset1);
-  uvw[i3] = 2.0 * (xyz[i3] - offset2);
-  uvw[i1] = 0.0;
+  Angle a;
+  a[i2] = 2.0 * (r[i2] - offset1);
+  a[i3] = 2.0 * (r[i3] - offset2);
+  a[i1] = 0.0;
+  return a;
 }
 
 //==============================================================================
@@ -620,21 +585,20 @@ SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, y0, z0, r);
 }
 
-inline double SurfaceXCylinder::evaluate(const double xyz[3]) const
+double SurfaceXCylinder::evaluate(Position r) const
 {
-  return axis_aligned_cylinder_evaluate<1, 2>(xyz, y0, z0, r);
+  return axis_aligned_cylinder_evaluate<1, 2>(r, y0, z0, this->r);
 }
 
-inline double SurfaceXCylinder::distance(const double xyz[3],
-     const double uvw[3], bool coincident) const
+double SurfaceXCylinder::distance(Position r, Angle a, bool coincident) const
 {
-  return axis_aligned_cylinder_distance<0, 1, 2>(xyz, uvw, coincident, y0, z0,
-                                                 r);
+  return axis_aligned_cylinder_distance<0, 1, 2>(r, a, coincident, y0, z0,
+                                                 this->r);
 }
 
-inline void SurfaceXCylinder::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceXCylinder::normal(Position r) const
 {
-  axis_aligned_cylinder_normal<0, 1, 2>(xyz, uvw, y0, z0);
+  return axis_aligned_cylinder_normal<0, 1, 2>(r, y0, z0);
 }
 
 
@@ -655,21 +619,20 @@ SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, x0, z0, r);
 }
 
-inline double SurfaceYCylinder::evaluate(const double xyz[3]) const
+double SurfaceYCylinder::evaluate(Position r) const
 {
-  return axis_aligned_cylinder_evaluate<0, 2>(xyz, x0, z0, r);
+  return axis_aligned_cylinder_evaluate<0, 2>(r, x0, z0, this->r);
 }
 
-inline double SurfaceYCylinder::distance(const double xyz[3],
-     const double uvw[3], bool coincident) const
+double SurfaceYCylinder::distance(Position r, Angle a, bool coincident) const
 {
-  return axis_aligned_cylinder_distance<1, 0, 2>(xyz, uvw, coincident, x0, z0,
-                                                 r);
+  return axis_aligned_cylinder_distance<1, 0, 2>(r, a, coincident, x0, z0,
+                                                 this->r);
 }
 
-inline void SurfaceYCylinder::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceYCylinder::normal(Position r) const
 {
-  axis_aligned_cylinder_normal<1, 0, 2>(xyz, uvw, x0, z0);
+  return axis_aligned_cylinder_normal<1, 0, 2>(r, x0, z0);
 }
 
 void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
@@ -689,21 +652,20 @@ SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, x0, y0, r);
 }
 
-inline double SurfaceZCylinder::evaluate(const double xyz[3]) const
+double SurfaceZCylinder::evaluate(Position r) const
 {
-  return axis_aligned_cylinder_evaluate<0, 1>(xyz, x0, y0, r);
+  return axis_aligned_cylinder_evaluate<0, 1>(r, x0, y0, this->r);
 }
 
-inline double SurfaceZCylinder::distance(const double xyz[3],
-     const double uvw[3], bool coincident) const
+double SurfaceZCylinder::distance(Position r, Angle a, bool coincident) const
 {
-  return axis_aligned_cylinder_distance<2, 0, 1>(xyz, uvw, coincident, x0, y0,
-                                                 r);
+  return axis_aligned_cylinder_distance<2, 0, 1>(r, a, coincident, x0, y0,
+                                                 this->r);
 }
 
-inline void SurfaceZCylinder::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceZCylinder::normal(Position r) const
 {
-  axis_aligned_cylinder_normal<2, 0, 1>(xyz, uvw, x0, y0);
+  return axis_aligned_cylinder_normal<2, 0, 1>(r, x0, y0);
 }
 
 void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
@@ -723,22 +685,21 @@ SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, x0, y0, z0, r);
 }
 
-double SurfaceSphere::evaluate(const double xyz[3]) const
+double SurfaceSphere::evaluate(Position r) const
 {
-  const double x = xyz[0] - x0;
-  const double y = xyz[1] - y0;
-  const double z = xyz[2] - z0;
-  return x*x + y*y + z*z - r*r;
+  const double x = r.x - x0;
+  const double y = r.y - y0;
+  const double z = r.z - z0;
+  return x*x + y*y + z*z - this->r*this->r;
 }
 
-double SurfaceSphere::distance(const double xyz[3], const double uvw[3],
-                               bool coincident) const
+double SurfaceSphere::distance(Position r, Angle a, bool coincident) const
 {
-  const double x = xyz[0] - x0;
-  const double y = xyz[1] - y0;
-  const double z = xyz[2] - z0;
-  const double k = x*uvw[0] + y*uvw[1] + z*uvw[2];
-  const double c = x*x + y*y + z*z - r*r;
+  const double x = r.x - x0;
+  const double y = r.y - y0;
+  const double z = r.z - z0;
+  const double k = x*a.u() + y*a.v() + z*a.w();
+  const double c = x*x + y*y + z*z - this->r*this->r;
   const double quad = k*k - c;
 
   if (quad < 0.0) {
@@ -770,11 +731,9 @@ double SurfaceSphere::distance(const double xyz[3], const double uvw[3],
   }
 }
 
-inline void SurfaceSphere::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceSphere::normal(Position r) const
 {
-  uvw[0] = 2.0 * (xyz[0] - x0);
-  uvw[1] = 2.0 * (xyz[1] - y0);
-  uvw[2] = 2.0 * (xyz[2] - z0);
+  return {2.0*(r.x - x0), 2.0*(r.y - y0), 2.0*(r.z - z0)};
 }
 
 void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
@@ -792,30 +751,30 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
 // The other two parameters indicate the other two axes.  offset1, offset2,
 // and offset3 should correspond with i1, i2, and i3, respectively.
 template<int i1, int i2, int i3> double
-axis_aligned_cone_evaluate(const double xyz[3], double offset1,
+axis_aligned_cone_evaluate(Position r, double offset1,
                            double offset2, double offset3, double radius_sq)
 {
-  const double xyz1 = xyz[i1] - offset1;
-  const double xyz2 = xyz[i2] - offset2;
-  const double xyz3 = xyz[i3] - offset3;
-  return xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1;
+  const double r1 = r[i1] - offset1;
+  const double r2 = r[i2] - offset2;
+  const double r3 = r[i3] - offset3;
+  return r2*r2 + r3*r3 - radius_sq*r1*r1;
 }
 
 // The first template parameter indicates which axis the cone is aligned to.
 // The other two parameters indicate the other two axes.  offset1, offset2,
 // and offset3 should correspond with i1, i2, and i3, respectively.
 template<int i1, int i2, int i3> double
-axis_aligned_cone_distance(const double xyz[3], const double uvw[3],
+axis_aligned_cone_distance(Position r, Angle u,
      bool coincident, double offset1, double offset2, double offset3,
      double radius_sq)
 {
-  const double xyz1 = xyz[i1] - offset1;
-  const double xyz2 = xyz[i2] - offset2;
-  const double xyz3 = xyz[i3] - offset3;
-  const double a = uvw[i2]*uvw[i2] + uvw[i3]*uvw[i3]
-                   - radius_sq*uvw[i1]*uvw[i1];
-  const double k = xyz2*uvw[i2] + xyz3*uvw[i3] - radius_sq*xyz1*uvw[i1];
-  const double c = xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1;
+  const double r1 = r[i1] - offset1;
+  const double r2 = r[i2] - offset2;
+  const double r3 = r[i3] - offset3;
+  const double a = u[i2]*u[i2] + u[i3]*u[i3]
+                   - radius_sq*u[i1]*u[i1];
+  const double k = r2*u[i2] + r3*u[i3] - radius_sq*r1*u[i1];
+  const double c = r2*r2 + r3*r3 - radius_sq*r1*r1;
   double quad = k*k - a*c;
 
   double d;
@@ -858,13 +817,15 @@ axis_aligned_cone_distance(const double xyz[3], const double uvw[3],
 // The first template parameter indicates which axis the cone is aligned to.
 // The other two parameters indicate the other two axes.  offset1, offset2,
 // and offset3 should correspond with i1, i2, and i3, respectively.
-template<int i1, int i2, int i3> void
-axis_aligned_cone_normal(const double xyz[3], double uvw[3], double offset1,
-                         double offset2, double offset3, double radius_sq)
+template<int i1, int i2, int i3> Angle
+axis_aligned_cone_normal(Position r, double offset1, double offset2,
+                         double offset3, double radius_sq)
 {
-  uvw[i1] = -2.0 * radius_sq * (xyz[i1] - offset1);
-  uvw[i2] = 2.0 * (xyz[i2] - offset2);
-  uvw[i3] = 2.0 * (xyz[i3] - offset3);
+  Angle a;
+  a[i1] = -2.0 * radius_sq * (r[i1] - offset1);
+  a[i2] = 2.0 * (r[i2] - offset2);
+  a[i3] = 2.0 * (r[i3] - offset3);
+  return a;
 }
 
 //==============================================================================
@@ -877,21 +838,20 @@ SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, x0, y0, z0, r_sq);
 }
 
-inline double SurfaceXCone::evaluate(const double xyz[3]) const
+double SurfaceXCone::evaluate(Position r) const
 {
-  return axis_aligned_cone_evaluate<0, 1, 2>(xyz, x0, y0, z0, r_sq);
+  return axis_aligned_cone_evaluate<0, 1, 2>(r, x0, y0, z0, r_sq);
 }
 
-inline double SurfaceXCone::distance(const double xyz[3],
-     const double uvw[3], bool coincident) const
+double SurfaceXCone::distance(Position r, Angle a, bool coincident) const
 {
-  return axis_aligned_cone_distance<0, 1, 2>(xyz, uvw, coincident, x0, y0, z0,
+  return axis_aligned_cone_distance<0, 1, 2>(r, a, coincident, x0, y0, z0,
                                              r_sq);
 }
 
-inline void SurfaceXCone::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceXCone::normal(Position r) const
 {
-  axis_aligned_cone_normal<0, 1, 2>(xyz, uvw, x0, y0, z0, r_sq);
+  return axis_aligned_cone_normal<0, 1, 2>(r, x0, y0, z0, r_sq);
 }
 
 void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
@@ -911,21 +871,20 @@ SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, x0, y0, z0, r_sq);
 }
 
-inline double SurfaceYCone::evaluate(const double xyz[3]) const
+double SurfaceYCone::evaluate(Position r) const
 {
-  return axis_aligned_cone_evaluate<1, 0, 2>(xyz, y0, x0, z0, r_sq);
+  return axis_aligned_cone_evaluate<1, 0, 2>(r, y0, x0, z0, r_sq);
 }
 
-inline double SurfaceYCone::distance(const double xyz[3],
-     const double uvw[3], bool coincident) const
+double SurfaceYCone::distance(Position r, Angle a, bool coincident) const
 {
-  return axis_aligned_cone_distance<1, 0, 2>(xyz, uvw, coincident, y0, x0, z0,
+  return axis_aligned_cone_distance<1, 0, 2>(r, a, coincident, y0, x0, z0,
                                              r_sq);
 }
 
-inline void SurfaceYCone::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceYCone::normal(Position r) const
 {
-  axis_aligned_cone_normal<1, 0, 2>(xyz, uvw, y0, x0, z0, r_sq);
+  return axis_aligned_cone_normal<1, 0, 2>(r, y0, x0, z0, r_sq);
 }
 
 void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
@@ -945,21 +904,20 @@ SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node)
   read_coeffs(surf_node, id, x0, y0, z0, r_sq);
 }
 
-inline double SurfaceZCone::evaluate(const double xyz[3]) const
+double SurfaceZCone::evaluate(Position r) const
 {
-  return axis_aligned_cone_evaluate<2, 0, 1>(xyz, z0, x0, y0, r_sq);
+  return axis_aligned_cone_evaluate<2, 0, 1>(r, z0, x0, y0, r_sq);
 }
 
-inline double SurfaceZCone::distance(const double xyz[3],
-     const double uvw[3], bool coincident) const
+double SurfaceZCone::distance(Position r, Angle a, bool coincident) const
 {
-  return axis_aligned_cone_distance<2, 0, 1>(xyz, uvw, coincident, z0, x0, y0,
+  return axis_aligned_cone_distance<2, 0, 1>(r, a, coincident, z0, x0, y0,
                                              r_sq);
 }
 
-inline void SurfaceZCone::normal(const double xyz[3], double uvw[3]) const
+Angle SurfaceZCone::normal(Position r) const
 {
-  axis_aligned_cone_normal<2, 0, 1>(xyz, uvw, z0, x0, y0, r_sq);
+  return axis_aligned_cone_normal<2, 0, 1>(r, z0, x0, y0, r_sq);
 }
 
 void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
@@ -980,26 +938,25 @@ SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node)
 }
 
 double
-SurfaceQuadric::evaluate(const double xyz[3]) const
+SurfaceQuadric::evaluate(Position r) const
 {
-  const double &x = xyz[0];
-  const double &y = xyz[1];
-  const double &z = xyz[2];
+  const double x = r.x;
+  const double y = r.y;
+  const double z = r.z;
   return x*(A*x + D*y + G) +
          y*(B*y + E*z + H) +
          z*(C*z + F*x + J) + K;
 }
 
 double
-SurfaceQuadric::distance(const double xyz[3],
-     const double uvw[3], bool coincident) const
+SurfaceQuadric::distance(Position r, Angle ang, bool coincident) const
 {
-  const double &x = xyz[0];
-  const double &y = xyz[1];
-  const double &z = xyz[2];
-  const double &u = uvw[0];
-  const double &v = uvw[1];
-  const double &w = uvw[2];
+  const double &x = r.x;
+  const double &y = r.y;
+  const double &z = r.z;
+  const double &u = ang.u();
+  const double &v = ang.v();
+  const double &w = ang.w();
 
   const double a = A*u*u + B*v*v + C*w*w + D*u*v + E*v*w + F*u*w;
   const double k = (A*u*x + B*v*y + C*w*z + 0.5*(D*(u*y + v*x) +
@@ -1045,15 +1002,15 @@ SurfaceQuadric::distance(const double xyz[3],
   return d;
 }
 
-void
-SurfaceQuadric::normal(const double xyz[3], double uvw[3]) const
+Angle
+SurfaceQuadric::normal(Position r) const
 {
-  const double &x = xyz[0];
-  const double &y = xyz[1];
-  const double &z = xyz[2];
-  uvw[0] = 2.0*A*x + D*y + F*z + G;
-  uvw[1] = 2.0*B*y + D*x + E*z + H;
-  uvw[2] = 2.0*C*z + E*y + F*x + J;
+  const double &x = r.x;
+  const double &y = r.y;
+  const double &z = r.z;
+  return {2.0*A*x + D*y + F*z + G,
+      2.0*B*y + D*x + E*z + H,
+      2.0*C*z + E*y + F*x + J};
 }
 
 void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
@@ -1263,10 +1220,24 @@ extern "C" {
   int surface_bc(Surface *surf) {return surf->bc;}
 
   void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
-  {surf->reflect(xyz, uvw);}
+  {
+    Position r {xyz};
+    Angle a {uvw};
+    a = surf->reflect(r, a);
+
+    uvw[0] = a.u();
+    uvw[1] = a.v();
+    uvw[2] = a.w();
+  }
 
   void surface_normal(Surface *surf, double xyz[3], double uvw[3])
-  {return surf->normal(xyz, uvw);}
+  {
+    Position r {xyz};
+    Angle a = surf->normal(r);
+    uvw[0] = a.u();
+    uvw[1] = a.v();
+    uvw[2] = a.w();
+  }
 
   void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);}
 
@@ -1275,7 +1246,21 @@ extern "C" {
   bool
   surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
                    double uvw[3])
-  {return surf->periodic_translate(other, xyz, uvw);}
+  {
+    Position r {xyz};
+    Angle a {uvw};
+    bool rotational = surf->periodic_translate(other, r, a);
+
+    // Copy back to arrays
+    xyz[0] = r.x;
+    xyz[1] = r.y;
+    xyz[2] = r.z;
+    uvw[0] = a.u();
+    uvw[1] = a.v();
+    uvw[2] = a.w();
+
+    return rotational;
+  }
 
   void free_memory_surfaces_c()
   {

src/surface.h

@@ -9,6 +9,7 @@
 #include "pugixml.hpp"
 
 #include "constants.h"
+#include "geometry.h"
 
 
 namespace openmc {
@@ -65,39 +66,38 @@ public:
   virtual ~Surface() {}
 
   //! Determine which side of a surface a point lies on.
-  //! @param xyz[3] The 3D Cartesian coordinate of a point.
-  //! @param uvw[3] A direction used to "break ties" and pick a sense when the
+  //! @param r The 3D Cartesian coordinate of a point.
+  //! @param o A direction used to "break ties" and pick a sense when the
   //!   point is very close to the surface.
   //! @return true if the point is on the "positive" side of the surface and
   //!   false otherwise.
-  bool sense(const double xyz[3], const double uvw[3]) const;
+  bool sense(Position r, Angle a) const;
 
   //! Determine the direction of a ray reflected from the surface.
-  //! @param xyz[3] The point at which the ray is incident.
-  //! @param uvw[3] A direction.  This is both an input and an output parameter.
+  //! @param r The point at which the ray is incident.
+  //! @param o A direction.  This is both an input and an output parameter.
   //!   It specifies the icident direction on input and the reflected direction
   //!   on output.
-  void reflect(const double xyz[3], double uvw[3]) const;
+  Angle reflect(Position r, Angle a) const;
 
   //! Evaluate the equation describing the surface.
   //!
   //! Surfaces can be described by some function f(x, y, z) = 0.  This member
   //! function evaluates that mathematical function.
-  //! @param xyz[3] A 3D Cartesian coordinate.
-  virtual double evaluate(const double xyz[3]) const = 0;
+  //! @param r A 3D Cartesian coordinate.
+  virtual double evaluate(Position r) const = 0;
 
   //! Compute the distance between a point and the surface along a ray.
-  //! @param xyz[3] A 3D Cartesian coordinate.
-  //! @param uvw[3] The direction of the ray.
+  //! @param r A 3D Cartesian coordinate.
+  //! @param o The direction of the ray.
   //! @param coincident A hint to the code that the given point should lie
   //!   exactly on the surface.
-  virtual double distance(const double xyz[3], const double uvw[3],
-                          bool coincident) const = 0;
+  virtual double distance(Position r, Angle a, bool coincident) const = 0;
 
   //! Compute the local outward normal direction of the surface.
-  //! @param xyz[3] A 3D Cartesian coordinate.
-  //! @param uvw[3] This output argument provides the normal.
-  virtual void normal(const double xyz[3], double uvw[3]) const = 0;
+  //! @param r A 3D Cartesian coordinate.
+  //! @return Normal direction
+  virtual Angle normal(Position r) const = 0;
 
   //! Write all information needed to reconstruct the surface to an HDF5 group.
   //! @param group_id An HDF5 group id.
@@ -125,14 +125,14 @@ public:
 
   //! Translate a particle onto this surface from a periodic partner surface.
   //! @param other A pointer to the partner surface in this periodic BC.
-  //! @param xyz[3] A point on the partner surface that will be translated onto
+  //! @param r A point on the partner surface that will be translated onto
   //!   this surface.
-  //! @param uvw[3] A direction that will be rotated for systems with rotational
+  //! @param a A direction that will be rotated for systems with rotational
   //!   periodicity.
   //! @return true if this surface and its partner make a rotationally-periodic
   //!   boundary condition.
-  virtual bool periodic_translate(PeriodicSurface *other, double xyz[3],
-                                  double uvw[3]) const = 0;
+  virtual bool periodic_translate(const PeriodicSurface *other, Position& r,
+                                  Angle& a) const = 0;
 
   //! Get the bounding box for this surface.
   virtual BoundingBox bounding_box() const = 0;
@@ -149,12 +149,11 @@ class SurfaceXPlane : public PeriodicSurface
   double x0;
 public:
   explicit SurfaceXPlane(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3], bool coincident)
-         const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
+  bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
        const;
   BoundingBox bounding_box() const;
 };
@@ -170,12 +169,11 @@ class SurfaceYPlane : public PeriodicSurface
   double y0;
 public:
   explicit SurfaceYPlane(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
+  bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
        const;
   BoundingBox bounding_box() const;
 };
@@ -191,12 +189,11 @@ class SurfaceZPlane : public PeriodicSurface
   double z0;
 public:
   explicit SurfaceZPlane(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
+  bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
        const;
   BoundingBox bounding_box() const;
 };
@@ -212,12 +209,11 @@ class SurfacePlane : public PeriodicSurface
   double A, B, C, D;
 public:
   explicit SurfacePlane(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3], bool coincident)
-         const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3])
+  bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
        const;
   BoundingBox bounding_box() const;
 };
@@ -234,10 +230,9 @@ class SurfaceXCylinder : public Surface
   double y0, z0, r;
 public:
   explicit SurfaceXCylinder(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -253,10 +248,9 @@ class SurfaceYCylinder : public Surface
   double x0, z0, r;
 public:
   explicit SurfaceYCylinder(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -272,10 +266,9 @@ class SurfaceZCylinder : public Surface
   double x0, y0, r;
 public:
   explicit SurfaceZCylinder(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -291,10 +284,9 @@ class SurfaceSphere : public Surface
   double x0, y0, z0, r;
 public:
   explicit SurfaceSphere(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -310,10 +302,9 @@ class SurfaceXCone : public Surface
   double x0, y0, z0, r_sq;
 public:
   explicit SurfaceXCone(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -329,10 +320,9 @@ class SurfaceYCone : public Surface
   double x0, y0, z0, r_sq;
 public:
   explicit SurfaceYCone(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -348,10 +338,9 @@ class SurfaceZCone : public Surface
   double x0, y0, z0, r_sq;
 public:
   explicit SurfaceZCone(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -367,10 +356,9 @@ class SurfaceQuadric : public Surface
   double A, B, C, D, E, F, G, H, J, K;
 public:
   explicit SurfaceQuadric(pugi::xml_node surf_node);
-  double evaluate(const double xyz[3]) const;
-  double distance(const double xyz[3], const double uvw[3],
-                  bool coincident) const;
-  void normal(const double xyz[3], double uvw[3]) const;
+  double evaluate(Position r) const;
+  double distance(Position r, Angle a, bool coincident) const;
+  Angle normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };