diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0be9ccb32..c94f4df5f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -386,7 +386,6 @@ 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
@@ -396,6 +395,7 @@ add_library(libopenmc SHARED
   src/mgxs_interface.cpp
   src/particle.cpp
   src/plot.cpp
+  src/position.cpp
   src/pugixml/pugixml_c.cpp
   src/random_lcg.cpp
   src/scattdata.cpp
diff --git a/src/cell.cpp b/src/cell.cpp
index 972540fbd..acf1473bb 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -279,19 +279,19 @@ Cell::Cell(pugi::xml_node cell_node)
 //==============================================================================
 
 bool
-Cell::contains(Position r, Angle a, int32_t on_surface) const
+Cell::contains(Position r, Direction u, int32_t on_surface) const
 {
   if (simple) {
-    return contains_simple(r, a, on_surface);
+    return contains_simple(r, u, on_surface);
   } else {
-    return contains_complex(r, a, on_surface);
+    return contains_complex(r, u, on_surface);
   }
 }
 
 //==============================================================================
 
 std::pair<double, int32_t>
-Cell::distance(Position r, Angle a, int32_t on_surface) const
+Cell::distance(Position r, Direction u, int32_t on_surface) const
 {
   double min_dist {INFTY};
   int32_t i_surf {std::numeric_limits<int32_t>::max()};
@@ -303,7 +303,7 @@ Cell::distance(Position r, Angle a, int32_t on_surface) const
     // 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(r, a, coincident)};
+    double d {surfaces_c[abs(token)-1]->distance(r, u, coincident)};
 
     // Check if this distance is the new minimum.
     if (d < min_dist) {
@@ -354,7 +354,7 @@ Cell::to_hdf5(hid_t cell_group) const
 //==============================================================================
 
 bool
-Cell::contains_simple(Position r, Angle a, int32_t on_surface) const
+Cell::contains_simple(Position r, Direction u, int32_t on_surface) const
 {
   for (int32_t token : rpn) {
     if (token < OP_UNION) {
@@ -367,7 +367,7 @@ Cell::contains_simple(Position r, Angle a, int32_t on_surface) const
         return false;
       } else {
         // Note the off-by-one indexing
-        bool sense = surfaces_c[abs(token)-1]->sense(r, a);
+        bool sense = surfaces_c[abs(token)-1]->sense(r, u);
         if (sense != (token > 0)) {return false;}
       }
     }
@@ -378,7 +378,7 @@ Cell::contains_simple(Position r, Angle a, int32_t on_surface) const
 //==============================================================================
 
 bool
-Cell::contains_complex(Position r, Angle a, int32_t on_surface) const
+Cell::contains_complex(Position r, Direction u, 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.
@@ -409,7 +409,7 @@ Cell::contains_complex(Position r, Angle a, int32_t on_surface) const
         stack[i_stack] = false;
       } else {
         // Note the off-by-one indexing
-        bool sense = surfaces_c[abs(token)-1]->sense(r, a);;
+        bool sense = surfaces_c[abs(token)-1]->sense(r, u);;
         stack[i_stack] = (sense == (token > 0));
       }
     }
@@ -492,16 +492,16 @@ extern "C" {
   bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface)
   {
     Position r {xyz};
-    Angle a {uvw};
-    return c->contains(r, a, on_surface);
+    Direction u {uvw};
+    return c->contains(r, u, on_surface);
   }
 
   void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface,
                      double *min_dist, int32_t *i_surf)
   {
     Position r {xyz};
-    Angle a {uvw};
-    std::pair<double, int32_t> out = c->distance(r, a, on_surface);
+    Direction u {uvw};
+    std::pair<double, int32_t> out = c->distance(r, u, on_surface);
     *min_dist = out.first;
     *i_surf = out.second;
   }
diff --git a/src/cell.h b/src/cell.h
index 08a6abf3f..10265beb1 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -8,7 +8,8 @@
 
 #include "hdf5.h"
 #include "pugixml.hpp"
-#include "geometry.h"
+
+#include "position.h"
 
 
 namespace openmc {
@@ -98,19 +99,19 @@ public:
   //!   known to be on.  This index takes precedence over surface sense
   //!   calculations.
   bool
-  contains(Position r, Angle a, int32_t on_surface) const;
+  contains(Position r, Direction u, int32_t on_surface) const;
 
   //! Find the oncoming boundary of this cell.
   std::pair<double, int32_t>
-  distance(Position r, Angle a, int32_t on_surface) const;
+  distance(Position r, Direction u, 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(Position r, Angle a, int32_t on_surface) const;
-  bool contains_complex(Position r, Angle a, int32_t on_surface) const;
+  bool contains_simple(Position r, Direction u, int32_t on_surface) const;
+  bool contains_complex(Position r, Direction u, int32_t on_surface) const;
 };
 
 } // namespace openmc
diff --git a/src/geometry.h b/src/geometry.h
index e069b9315..09e7eda8e 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -1,68 +1,10 @@
-#ifndef GEOMETRY_H
-#define GEOMETRY_H
+#ifndef OPENMC_GEOMETRY_H
+#define OPENMC_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(const 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(const double uvw[]) : Position{uvw} { };
-  Angle(Position r) : Position{r} { };
-};
-
 } // namespace openmc
 
-#endif // GEOMETRY_H
+#endif // OPENMC_GEOMETRY_H
diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h
index 8b1f7a434..b0466e613 100644
--- a/src/hdf5_interface.h
+++ b/src/hdf5_interface.h
@@ -11,7 +11,7 @@
 #include <vector>
 #include <complex.h>
 
-#include "geometry.h"
+#include "position.h"
 
 
 namespace openmc {
diff --git a/src/lattice.cpp b/src/lattice.cpp
index c67fc419f..73dcb9866 100644
--- a/src/lattice.cpp
+++ b/src/lattice.cpp
@@ -223,25 +223,23 @@ RectLattice::are_valid_indices(const int i_xyz[3]) const
 //==============================================================================
 
 std::pair<double, std::array<int, 3>>
-RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const
+RectLattice::distance(Position r, Direction u, const int i_xyz[3]) const
 {
   // Get short aliases to the coordinates.
   double x = r.x;
   double y = r.y;
   double z = r.z;
-  double u = a.u();
-  double v = a.v();
 
   // Determine the oncoming edge.
-  double x0 {copysign(0.5 * pitch[0], u)};
-  double y0 {copysign(0.5 * pitch[1], v)};
+  double x0 {copysign(0.5 * pitch[0], u.x)};
+  double y0 {copysign(0.5 * pitch[1], u.y)};
 
   // Left and right sides
   double d {INFTY};
   std::array<int, 3> lattice_trans;
-  if ((std::abs(x - x0) > FP_PRECISION) && u != 0) {
-    d = (x0 - x) / u;
-    if (u > 0) {
+  if ((std::abs(x - x0) > FP_PRECISION) && u.x != 0) {
+    d = (x0 - x) / u.x;
+    if (u.x > 0) {
       lattice_trans = {1, 0, 0};
     } else {
       lattice_trans = {-1, 0, 0};
@@ -249,11 +247,11 @@ RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const
   }
 
   // Front and back sides
-  if ((std::abs(y - y0) > FP_PRECISION) && v != 0) {
-    double this_d = (y0 - y) / v;
+  if ((std::abs(y - y0) > FP_PRECISION) && u.y != 0) {
+    double this_d = (y0 - y) / u.y;
     if (this_d < d) {
       d = this_d;
-      if (v > 0) {
+      if (u.y > 0) {
         lattice_trans = {0, 1, 0};
       } else {
         lattice_trans = {0, -1, 0};
@@ -263,13 +261,12 @@ RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const
 
   // Top and bottom sides
   if (is_3d) {
-    double w {a.w()};
-    double z0 {copysign(0.5 * pitch[2], w)};
-    if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
-      double this_d = (z0 - z) / w;
+    double z0 {copysign(0.5 * pitch[2], u.z)};
+    if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) {
+      double this_d = (z0 - z) / u.z;
       if (this_d < d) {
         d = this_d;
-        if (w > 0) {
+        if (u.z > 0) {
           lattice_trans = {0, 0, 1};
         } else {
           lattice_trans = {0, 0, -1};
@@ -579,11 +576,11 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const
 //==============================================================================
 
 std::pair<double, std::array<int, 3>>
-HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const
+HexLattice::distance(Position r, Direction u, const int i_xyz[3]) const
 {
   // Compute the direction on the hexagonal basis.
-  double beta_dir = a.u() * std::sqrt(3.0) / 2.0 + a.v() / 2.0;
-  double gamma_dir = a.u() * std::sqrt(3.0) / 2.0 - a.v() / 2.0;
+  double beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0;
+  double gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0;
 
   // Note that hexagonal lattice distance calculations are performed
   // using the particle's coordinates relative to the neighbor lattice
@@ -636,18 +633,18 @@ HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const
   }
 
   // Upper and lower sides.
-  edge = -copysign(0.5*pitch[0], a.v());
-  if (a.v() > 0) {
+  edge = -copysign(0.5*pitch[0], u.y);
+  if (u.y > 0) {
     const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]};
     r_t = get_local_position(r, i_xyz_t);
   } else {
     const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]};
     r_t = get_local_position(r, i_xyz_t);
   }
-  if ((std::abs(r_t.y - edge) > FP_PRECISION) && a.v() != 0) {
-    double this_d = (edge - r_t.y) / a.v();
+  if ((std::abs(r_t.y - edge) > FP_PRECISION) && u.y != 0) {
+    double this_d = (edge - r_t.y) / u.y;
     if (this_d < d) {
-      if (a.v() > 0) {
+      if (u.y > 0) {
         lattice_trans = {0, 1, 0};
       } else {
         lattice_trans = {0, -1, 0};
@@ -659,13 +656,12 @@ HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const
   // Top and bottom sides
   if (is_3d) {
     double z = r.z;
-    double w = a.w();
-    double z0 {copysign(0.5 * pitch[1], w)};
-    if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
-      double this_d = (z0 - z) / w;
+    double z0 {copysign(0.5 * pitch[1], u.z)};
+    if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) {
+      double this_d = (z0 - z) / u.z;
       if (this_d < d) {
         d = this_d;
-        if (w > 0) {
+        if (u.z > 0) {
           lattice_trans = {0, 0, 1};
         } else {
           lattice_trans = {0, 0, -1};
@@ -897,8 +893,8 @@ extern "C" {
                         const int i_xyz[3], double *d, int lattice_trans[3])
   {
     Position r {xyz};
-    Angle a {uvw};
-    std::pair<double, std::array<int, 3>> ld {lat->distance(r, a, i_xyz)};
+    Direction u {uvw};
+    std::pair<double, std::array<int, 3>> ld {lat->distance(r, u, i_xyz)};
     *d = ld.first;
     lattice_trans[0] = ld.second[0];
     lattice_trans[1] = ld.second[1];
diff --git a/src/lattice.h b/src/lattice.h
index a6c449662..5bd15a612 100644
--- a/src/lattice.h
+++ b/src/lattice.h
@@ -7,11 +7,12 @@
 #include <unordered_map>
 #include <vector>
 
-#include "constants.h"
-#include "geometry.h"
 #include "hdf5.h"
 #include "pugixml.hpp"
 
+#include "constants.h"
+#include "position.h"
+
 
 namespace openmc {
 
@@ -82,7 +83,7 @@ public:
   //! @return The distance to the next crossing and an array indicating how the
   //!   lattice indices would change after crossing that boundary.
   virtual std::pair<double, std::array<int, 3>>
-  distance(Position r, Angle a, const int i_xyz[3]) const
+  distance(Position r, Direction u, const int i_xyz[3]) const
   = 0;
 
   //! \brief Find the lattice tile indices for a given point.
@@ -194,7 +195,7 @@ public:
   bool are_valid_indices(const int i_xyz[3]) const;
 
   std::pair<double, std::array<int, 3>>
-  distance(Position r, Angle a, const int i_xyz[3]) const;
+  distance(Position r, Direction u, const int i_xyz[3]) const;
 
   std::array<int, 3> get_indices(Position r) const;
 
@@ -234,7 +235,7 @@ public:
   bool are_valid_indices(const int i_xyz[3]) const;
 
   std::pair<double, std::array<int, 3>>
-  distance(Position r, Angle a, const int i_xyz[3]) const;
+  distance(Position r, Direction u, const int i_xyz[3]) const;
 
   std::array<int, 3> get_indices(Position r) const;
 
diff --git a/src/geometry.cpp b/src/position.cpp
similarity index 75%
rename from src/geometry.cpp
rename to src/position.cpp
index 53bf32b5c..85e511041 100644
--- a/src/geometry.cpp
+++ b/src/position.cpp
@@ -1,7 +1,11 @@
-#include "geometry.h"
+#include "position.h"
 
 namespace openmc {
 
+//==============================================================================
+// Position implementation
+//==============================================================================
+
 Position&
 Position::operator+=(Position other)
 {
diff --git a/src/position.h b/src/position.h
new file mode 100644
index 000000000..808e0bcf6
--- /dev/null
+++ b/src/position.h
@@ -0,0 +1,74 @@
+#ifndef OPENMC_POSITION_H
+#define OPENMC_POSITION_H
+
+namespace openmc {
+
+//==============================================================================
+//! Type representing a position in Cartesian coordinates
+//==============================================================================
+
+struct Position {
+  // Constructors
+  Position() = default;
+  Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { };
+  Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
+
+  // Unary operators
+  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;
+    }
+  }
+
+  // Other member functions
+
+  //! Dot product of two vectors
+  //! \param[in] other Vector to take dot product with
+  //! \result Resulting dot product
+  inline double dot(Position other) {
+    return x*other.x + y*other.y + z*other.z;
+  }
+
+  // Data members
+  double x = 0.;
+  double y = 0.;
+  double z = 0.;
+};
+
+// Binary operators
+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; }
+
+//==============================================================================
+//! Type representing a vector direction in Cartesian coordinates
+//==============================================================================
+
+using Direction = Position;
+
+} // namespace openmc
+
+#endif // OPENMC_POSITION_H
\ No newline at end of file
diff --git a/src/surface.cpp b/src/surface.cpp
index aec0972b8..825886152 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -178,7 +178,7 @@ Surface::Surface(pugi::xml_node surf_node)
 }
 
 bool
-Surface::sense(Position r, Angle a) const
+Surface::sense(Position r, Direction u) const
 {
   // Evaluate the surface equation at the particle's coordinates to determine
   // which side the particle is on.
@@ -189,22 +189,22 @@ Surface::sense(Position r, Angle a) const
     // 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.
-    return a.dot(normal(r)) > 0.0;
+    return u.dot(normal(r)) > 0.0;
   }
   return f > 0.0;
 }
 
-Angle
-Surface::reflect(Position r, Angle a) const
+Direction
+Surface::reflect(Position r, Direction u) const
 {
   // Determine projection of direction onto normal and squared magnitude of
   // normal.
-  Angle n = normal(r);
-  const double projection = n.dot(a);
+  Direction n = normal(r);
+  const double projection = n.dot(u);
   const double magnitude = n.dot(n);
 
   // Reflect direction according to normal.
-  return a -= (2.0 * projection / magnitude) * n;
+  return u -= (2.0 * projection / magnitude) * n;
 }
 
 void
@@ -257,11 +257,11 @@ PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node)
 
 // The template parameter indicates the axis normal to the plane.
 template<int i> double
-axis_aligned_plane_distance(Position r, Angle a, bool coincident, double offset)
+axis_aligned_plane_distance(Position r, Direction u, bool coincident, double offset)
 {
   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 (coincident or std::abs(f) < FP_COINCIDENT or u[i] == 0.0) return INFTY;
+  const double d = f / u[i];
   if (d < 0.0) return INFTY;
   return d;
 }
@@ -281,12 +281,12 @@ double SurfaceXPlane::evaluate(Position r) const
   return r.x - x0;
 }
 
-double SurfaceXPlane::distance(Position r, Angle a, bool coincident) const
+double SurfaceXPlane::distance(Position r, Direction u, bool coincident) const
 {
-  return axis_aligned_plane_distance<0>(r, a, coincident, x0);
+  return axis_aligned_plane_distance<0>(r, u, coincident, x0);
 }
 
-Angle SurfaceXPlane::normal(Position r) const
+Direction SurfaceXPlane::normal(Position r) const
 {
   return {1., 0., 0.};
 }
@@ -299,23 +299,23 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
 }
 
 bool SurfaceXPlane::periodic_translate(const PeriodicSurface *other, Position& r,
-                                       Angle& a) const
+                                       Direction& u) const
 {
-  Angle other_n = other->normal(r);
-  if (other_n.u() == 1 and other_n.v() == 0 and other_n.w() == 0) {
+  Direction other_n = other->normal(r);
+  if (other_n.x == 1 and other_n.y == 0 and other_n.z == 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 position/angle for rotational
+    // evaluate function to find y0, then adjust position/Direction for rotational
     // symmetry.
     double y0 = -other->evaluate({0., 0., 0.});
     r.y = r.x - x0 + y0;
     r.x = x0;
 
-    double u = a.u();
-    a.u() = -a.v();
-    a.v() = u;
+    double ux = u.x;
+    u.x = -u.y;
+    u.y = ux;
 
     return true;
   }
@@ -342,12 +342,12 @@ double SurfaceYPlane::evaluate(Position r) const
   return r.y - y0;
 }
 
-double SurfaceYPlane::distance(Position r, Angle a, bool coincident) const
+double SurfaceYPlane::distance(Position r, Direction u, bool coincident) const
 {
-  return axis_aligned_plane_distance<1>(r, a, coincident, y0);
+  return axis_aligned_plane_distance<1>(r, u, coincident, y0);
 }
 
-Angle SurfaceYPlane::normal(Position r) const
+Direction SurfaceYPlane::normal(Position r) const
 {
   return {0., 1., 0.};
 }
@@ -360,24 +360,24 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
 }
 
 bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, Position& r,
-                                       Angle& a) const
+                                       Direction& u) const
 {
-  Angle other_n = other->normal(r);
-  if (other_n.u() == 0 and other_n.v() == 1 and other_n.w() == 0) {
+  Direction other_n = other->normal(r);
+  if (other_n.x == 0 and other_n.y == 1 and other_n.z == 0) {
     // The periodic partner is also aligned along y.  Just change the y coord.
     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 position/angle for rotational
+    // evaluate function to find x0, then adjust position/Direction for rotational
     // symmetry.
     double x0 = -other->evaluate({0., 0., 0.});
     r.x = r.y - y0 + x0;
     r.y = y0;
 
-    double u = a.u();
-    a.u() = a.v();
-    a.v() = -u;
+    double ux = u.x;
+    u.x = u.y;
+    u.y = -ux;
 
     return true;
   }
@@ -404,12 +404,12 @@ double SurfaceZPlane::evaluate(Position r) const
   return r.z - z0;
 }
 
-double SurfaceZPlane::distance(Position r, Angle a, bool coincident) const
+double SurfaceZPlane::distance(Position r, Direction u, bool coincident) const
 {
-  return axis_aligned_plane_distance<2>(r, a, coincident, z0);
+  return axis_aligned_plane_distance<2>(r, u, coincident, z0);
 }
 
-Angle SurfaceZPlane::normal(Position r) const
+Direction SurfaceZPlane::normal(Position r) const
 {
   return {0., 0., 1.};
 }
@@ -422,7 +422,7 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
 }
 
 bool SurfaceZPlane::periodic_translate(const PeriodicSurface *other, Position& r,
-                                       Angle& a) const
+                                       Direction& u) const
 {
   // Assume the other plane is aligned along z.  Just change the z coord.
   r.z = z0;
@@ -452,10 +452,10 @@ SurfacePlane::evaluate(Position r) const
 }
 
 double
-SurfacePlane::distance(Position r, Angle a, bool coincident) const
+SurfacePlane::distance(Position r, Direction u, bool coincident) const
 {
   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();
+  const double projection = A*u.x + B*u.y + C*u.z;
   if (coincident or std::abs(f) < FP_COINCIDENT or projection == 0.0) {
     return INFTY;
   } else {
@@ -465,7 +465,7 @@ SurfacePlane::distance(Position r, Angle a, bool coincident) const
   }
 }
 
-Angle
+Direction
 SurfacePlane::normal(Position r) const
 {
   return {A, B, C};
@@ -479,7 +479,7 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const
 }
 
 bool SurfacePlane::periodic_translate(const PeriodicSurface *other, Position& r,
-                                      Angle& a) const
+                                      Direction& u) const
 {
   // This function assumes the other plane shares this plane's normal direction.
 
@@ -520,7 +520,7 @@ axis_aligned_cylinder_evaluate(Position r, double offset1,
 // 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(Position r, Angle u,
+axis_aligned_cylinder_distance(Position r, Direction u,
      bool coincident, double offset1, double offset2, double radius)
 {
   const double a = 1.0 - u[i1]*u[i1];  // u^2 + v^2
@@ -565,14 +565,14 @@ axis_aligned_cylinder_distance(Position r, Angle u,
 // 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> Angle
+template<int i1, int i2, int i3> Direction
 axis_aligned_cylinder_normal(Position r, double offset1, double offset2)
 {
-  Angle a;
-  a[i2] = 2.0 * (r[i2] - offset1);
-  a[i3] = 2.0 * (r[i3] - offset2);
-  a[i1] = 0.0;
-  return a;
+  Direction u;
+  u[i2] = 2.0 * (r[i2] - offset1);
+  u[i3] = 2.0 * (r[i3] - offset2);
+  u[i1] = 0.0;
+  return u;
 }
 
 //==============================================================================
@@ -590,13 +590,13 @@ double SurfaceXCylinder::evaluate(Position r) const
   return axis_aligned_cylinder_evaluate<1, 2>(r, y0, z0, this->r);
 }
 
-double SurfaceXCylinder::distance(Position r, Angle a, bool coincident) const
+double SurfaceXCylinder::distance(Position r, Direction u, bool coincident) const
 {
-  return axis_aligned_cylinder_distance<0, 1, 2>(r, a, coincident, y0, z0,
+  return axis_aligned_cylinder_distance<0, 1, 2>(r, u, coincident, y0, z0,
                                                  this->r);
 }
 
-Angle SurfaceXCylinder::normal(Position r) const
+Direction SurfaceXCylinder::normal(Position r) const
 {
   return axis_aligned_cylinder_normal<0, 1, 2>(r, y0, z0);
 }
@@ -624,13 +624,13 @@ double SurfaceYCylinder::evaluate(Position r) const
   return axis_aligned_cylinder_evaluate<0, 2>(r, x0, z0, this->r);
 }
 
-double SurfaceYCylinder::distance(Position r, Angle a, bool coincident) const
+double SurfaceYCylinder::distance(Position r, Direction u, bool coincident) const
 {
-  return axis_aligned_cylinder_distance<1, 0, 2>(r, a, coincident, x0, z0,
+  return axis_aligned_cylinder_distance<1, 0, 2>(r, u, coincident, x0, z0,
                                                  this->r);
 }
 
-Angle SurfaceYCylinder::normal(Position r) const
+Direction SurfaceYCylinder::normal(Position r) const
 {
   return axis_aligned_cylinder_normal<1, 0, 2>(r, x0, z0);
 }
@@ -657,13 +657,13 @@ double SurfaceZCylinder::evaluate(Position r) const
   return axis_aligned_cylinder_evaluate<0, 1>(r, x0, y0, this->r);
 }
 
-double SurfaceZCylinder::distance(Position r, Angle a, bool coincident) const
+double SurfaceZCylinder::distance(Position r, Direction u, bool coincident) const
 {
-  return axis_aligned_cylinder_distance<2, 0, 1>(r, a, coincident, x0, y0,
+  return axis_aligned_cylinder_distance<2, 0, 1>(r, u, coincident, x0, y0,
                                                  this->r);
 }
 
-Angle SurfaceZCylinder::normal(Position r) const
+Direction SurfaceZCylinder::normal(Position r) const
 {
   return axis_aligned_cylinder_normal<2, 0, 1>(r, x0, y0);
 }
@@ -693,12 +693,12 @@ double SurfaceSphere::evaluate(Position r) const
   return x*x + y*y + z*z - this->r*this->r;
 }
 
-double SurfaceSphere::distance(Position r, Angle a, bool coincident) const
+double SurfaceSphere::distance(Position r, Direction u, bool coincident) const
 {
   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 k = x*u.x + y*u.y + z*u.z;
   const double c = x*x + y*y + z*z - this->r*this->r;
   const double quad = k*k - c;
 
@@ -731,7 +731,7 @@ double SurfaceSphere::distance(Position r, Angle a, bool coincident) const
   }
 }
 
-Angle SurfaceSphere::normal(Position r) const
+Direction SurfaceSphere::normal(Position r) const
 {
   return {2.0*(r.x - x0), 2.0*(r.y - y0), 2.0*(r.z - z0)};
 }
@@ -764,7 +764,7 @@ axis_aligned_cone_evaluate(Position r, double offset1,
 // 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(Position r, Angle u,
+axis_aligned_cone_distance(Position r, Direction u,
      bool coincident, double offset1, double offset2, double offset3,
      double radius_sq)
 {
@@ -817,15 +817,15 @@ axis_aligned_cone_distance(Position r, Angle u,
 // 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> Angle
+template<int i1, int i2, int i3> Direction
 axis_aligned_cone_normal(Position r, double offset1, double offset2,
                          double offset3, double radius_sq)
 {
-  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;
+  Direction u;
+  u[i1] = -2.0 * radius_sq * (r[i1] - offset1);
+  u[i2] = 2.0 * (r[i2] - offset2);
+  u[i3] = 2.0 * (r[i3] - offset3);
+  return u;
 }
 
 //==============================================================================
@@ -843,13 +843,13 @@ double SurfaceXCone::evaluate(Position r) const
   return axis_aligned_cone_evaluate<0, 1, 2>(r, x0, y0, z0, r_sq);
 }
 
-double SurfaceXCone::distance(Position r, Angle a, bool coincident) const
+double SurfaceXCone::distance(Position r, Direction u, bool coincident) const
 {
-  return axis_aligned_cone_distance<0, 1, 2>(r, a, coincident, x0, y0, z0,
+  return axis_aligned_cone_distance<0, 1, 2>(r, u, coincident, x0, y0, z0,
                                              r_sq);
 }
 
-Angle SurfaceXCone::normal(Position r) const
+Direction SurfaceXCone::normal(Position r) const
 {
   return axis_aligned_cone_normal<0, 1, 2>(r, x0, y0, z0, r_sq);
 }
@@ -876,13 +876,13 @@ double SurfaceYCone::evaluate(Position r) const
   return axis_aligned_cone_evaluate<1, 0, 2>(r, y0, x0, z0, r_sq);
 }
 
-double SurfaceYCone::distance(Position r, Angle a, bool coincident) const
+double SurfaceYCone::distance(Position r, Direction u, bool coincident) const
 {
-  return axis_aligned_cone_distance<1, 0, 2>(r, a, coincident, y0, x0, z0,
+  return axis_aligned_cone_distance<1, 0, 2>(r, u, coincident, y0, x0, z0,
                                              r_sq);
 }
 
-Angle SurfaceYCone::normal(Position r) const
+Direction SurfaceYCone::normal(Position r) const
 {
   return axis_aligned_cone_normal<1, 0, 2>(r, y0, x0, z0, r_sq);
 }
@@ -909,13 +909,13 @@ double SurfaceZCone::evaluate(Position r) const
   return axis_aligned_cone_evaluate<2, 0, 1>(r, z0, x0, y0, r_sq);
 }
 
-double SurfaceZCone::distance(Position r, Angle a, bool coincident) const
+double SurfaceZCone::distance(Position r, Direction u, bool coincident) const
 {
-  return axis_aligned_cone_distance<2, 0, 1>(r, a, coincident, z0, x0, y0,
+  return axis_aligned_cone_distance<2, 0, 1>(r, u, coincident, z0, x0, y0,
                                              r_sq);
 }
 
-Angle SurfaceZCone::normal(Position r) const
+Direction SurfaceZCone::normal(Position r) const
 {
   return axis_aligned_cone_normal<2, 0, 1>(r, z0, x0, y0, r_sq);
 }
@@ -949,14 +949,14 @@ SurfaceQuadric::evaluate(Position r) const
 }
 
 double
-SurfaceQuadric::distance(Position r, Angle ang, bool coincident) const
+SurfaceQuadric::distance(Position r, Direction ang, bool coincident) const
 {
   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 &u = ang.x;
+  const double &v = ang.y;
+  const double &w = ang.z;
 
   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) +
@@ -1002,7 +1002,7 @@ SurfaceQuadric::distance(Position r, Angle ang, bool coincident) const
   return d;
 }
 
-Angle
+Direction
 SurfaceQuadric::normal(Position r) const
 {
   const double &x = r.x;
@@ -1222,21 +1222,21 @@ extern "C" {
   void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
   {
     Position r {xyz};
-    Angle a {uvw};
-    a = surf->reflect(r, a);
+    Direction u {uvw};
+    u = surf->reflect(r, u);
 
-    uvw[0] = a.u();
-    uvw[1] = a.v();
-    uvw[2] = a.w();
+    uvw[0] = u.x;
+    uvw[1] = u.y;
+    uvw[2] = u.z;
   }
 
   void surface_normal(Surface *surf, double xyz[3], double uvw[3])
   {
     Position r {xyz};
-    Angle a = surf->normal(r);
-    uvw[0] = a.u();
-    uvw[1] = a.v();
-    uvw[2] = a.w();
+    Direction u = surf->normal(r);
+    uvw[0] = u.x;
+    uvw[1] = u.y;
+    uvw[2] = u.z;
   }
 
   void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);}
@@ -1248,16 +1248,16 @@ extern "C" {
                    double uvw[3])
   {
     Position r {xyz};
-    Angle a {uvw};
-    bool rotational = surf->periodic_translate(other, r, a);
+    Direction u {uvw};
+    bool rotational = surf->periodic_translate(other, r, u);
 
     // 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();
+    uvw[0] = u.x;
+    uvw[1] = u.y;
+    uvw[2] = u.z;
 
     return rotational;
   }
diff --git a/src/surface.h b/src/surface.h
index be3ad52ab..7d639c57a 100644
--- a/src/surface.h
+++ b/src/surface.h
@@ -9,7 +9,7 @@
 #include "pugixml.hpp"
 
 #include "constants.h"
-#include "geometry.h"
+#include "position.h"
 
 
 namespace openmc {
@@ -71,14 +71,14 @@ public:
   //!   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(Position r, Angle a) const;
+  bool sense(Position r, Direction u) const;
 
   //! Determine the direction of a ray reflected from the surface.
   //! @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.
-  Angle reflect(Position r, Angle a) const;
+  Direction reflect(Position r, Direction u) const;
 
   //! Evaluate the equation describing the surface.
   //!
@@ -92,12 +92,12 @@ public:
   //! @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(Position r, Angle a, bool coincident) const = 0;
+  virtual double distance(Position r, Direction u, bool coincident) const = 0;
 
   //! Compute the local outward normal direction of the surface.
   //! @param r A 3D Cartesian coordinate.
   //! @return Normal direction
-  virtual Angle normal(Position r) const = 0;
+  virtual Direction normal(Position r) const = 0;
 
   //! Write all information needed to reconstruct the surface to an HDF5 group.
   //! @param group_id An HDF5 group id.
@@ -132,7 +132,7 @@ public:
   //! @return true if this surface and its partner make a rotationally-periodic
   //!   boundary condition.
   virtual bool periodic_translate(const PeriodicSurface *other, Position& r,
-                                  Angle& a) const = 0;
+                                  Direction& u) const = 0;
 
   //! Get the bounding box for this surface.
   virtual BoundingBox bounding_box() const = 0;
@@ -150,10 +150,10 @@ class SurfaceXPlane : public PeriodicSurface
 public:
   explicit SurfaceXPlane(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
+  bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
        const;
   BoundingBox bounding_box() const;
 };
@@ -170,10 +170,10 @@ class SurfaceYPlane : public PeriodicSurface
 public:
   explicit SurfaceYPlane(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
+  bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
        const;
   BoundingBox bounding_box() const;
 };
@@ -190,10 +190,10 @@ class SurfaceZPlane : public PeriodicSurface
 public:
   explicit SurfaceZPlane(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
+  bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
        const;
   BoundingBox bounding_box() const;
 };
@@ -210,10 +210,10 @@ class SurfacePlane : public PeriodicSurface
 public:
   explicit SurfacePlane(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
-  bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a)
+  bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u)
        const;
   BoundingBox bounding_box() const;
 };
@@ -231,8 +231,8 @@ class SurfaceXCylinder : public Surface
 public:
   explicit SurfaceXCylinder(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -249,8 +249,8 @@ class SurfaceYCylinder : public Surface
 public:
   explicit SurfaceYCylinder(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -267,8 +267,8 @@ class SurfaceZCylinder : public Surface
 public:
   explicit SurfaceZCylinder(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -285,8 +285,8 @@ class SurfaceSphere : public Surface
 public:
   explicit SurfaceSphere(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -303,8 +303,8 @@ class SurfaceXCone : public Surface
 public:
   explicit SurfaceXCone(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -321,8 +321,8 @@ class SurfaceYCone : public Surface
 public:
   explicit SurfaceYCone(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -339,8 +339,8 @@ class SurfaceZCone : public Surface
 public:
   explicit SurfaceZCone(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -357,8 +357,8 @@ class SurfaceQuadric : public Surface
 public:
   explicit SurfaceQuadric(pugi::xml_node surf_node);
   double evaluate(Position r) const;
-  double distance(Position r, Angle a, bool coincident) const;
-  Angle normal(Position r) const;
+  double distance(Position r, Direction u, bool coincident) const;
+  Direction normal(Position r) const;
   void to_hdf5_inner(hid_t group_id) const;
 };
 
@@ -367,19 +367,19 @@ public:
 //==============================================================================
 
 extern "C" {
-    Surface* surface_pointer(int surf_ind);
-    int surface_id(Surface *surf);
-    int surface_bc(Surface *surf);
-    bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
-    void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
-    double surface_distance(Surface *surf, double xyz[3], double uvw[3],
-                            bool coincident);
-    void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
-    void surface_to_hdf5(Surface *surf, hid_t group);
-    int surface_i_periodic(PeriodicSurface *surf);
-    bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
-                          double xyz[3], double uvw[3]);
-    void free_memory_surfaces_c();
+  Surface* surface_pointer(int surf_ind);
+  int surface_id(Surface *surf);
+  int surface_bc(Surface *surf);
+  bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
+  void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
+  double surface_distance(Surface *surf, double xyz[3], double uvw[3],
+                          bool coincident);
+  void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
+  void surface_to_hdf5(Surface *surf, hid_t group);
+  int surface_i_periodic(PeriodicSurface *surf);
+  bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
+                        double xyz[3], double uvw[3]);
+  void free_memory_surfaces_c();
 }
 
 } // namespace openmc