diff --git a/CMakeLists.txt b/CMakeLists.txt
index 10867384b14..20068d7e760 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -407,6 +407,7 @@ list(APPEND libopenmc_SOURCES
   src/random_ray/linear_source_domain.cpp
   src/random_ray/moment_matrix.cpp
   src/random_ray/source_region.cpp
+  src/ray.cpp
   src/reaction.cpp
   src/reaction_product.cpp
   src/scattdata.cpp
diff --git a/include/openmc/plot.h b/include/openmc/plot.h
index aa373945317..6c00fc2f6ba 100644
--- a/include/openmc/plot.h
+++ b/include/openmc/plot.h
@@ -17,6 +17,7 @@
 #include "openmc/particle.h"
 #include "openmc/position.h"
 #include "openmc/random_lcg.h"
+#include "openmc/ray.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -497,47 +498,6 @@ class SolidRayTracePlot : public RayTracePlot {
   Position light_location_;
 };
 
-// Base class that implements ray tracing logic, not necessarily through
-// defined regions of the geometry but also outside of it.
-class Ray : public GeometryState {
-
-public:
-  // Initialize from location and direction
-  Ray(Position r, Direction u) { init_from_r_u(r, u); }
-
-  // Initialize from known geometry state
-  Ray(const GeometryState& p) : GeometryState(p) {}
-
-  // Called at every surface intersection within the model
-  virtual void on_intersection() = 0;
-
-  /*
-   * Traces the ray through the geometry, calling on_intersection
-   * at every surface boundary.
-   */
-  void trace();
-
-  // Stops the ray and exits tracing when called from on_intersection
-  void stop() { stop_ = true; }
-
-  // Sets the dist_ variable
-  void compute_distance();
-
-protected:
-  // Records how far the ray has traveled
-  double traversal_distance_ {0.0};
-
-private:
-  // Max intersections before we assume ray tracing is caught in an infinite
-  // loop:
-  static const int MAX_INTERSECTIONS = 1000000;
-
-  bool hit_something_ {false};
-  bool stop_ {false};
-
-  unsigned event_counter_ {0};
-};
-
 class ProjectionRay : public Ray {
 public:
   ProjectionRay(Position r, Direction u, const WireframeRayTracePlot& plot,
diff --git a/include/openmc/ray.h b/include/openmc/ray.h
new file mode 100644
index 00000000000..62e86b0d90f
--- /dev/null
+++ b/include/openmc/ray.h
@@ -0,0 +1,50 @@
+#ifndef OPENMC_RAY_H
+#define OPENMC_RAY_H
+
+#include "openmc/particle_data.h"
+#include "openmc/position.h"
+
+namespace openmc {
+
+// Base class that implements ray tracing logic, not necessarily through
+// defined regions of the geometry but also outside of it.
+class Ray : public GeometryState {
+
+public:
+  // Initialize from location and direction
+  Ray(Position r, Direction u) { init_from_r_u(r, u); }
+
+  // Initialize from known geometry state
+  Ray(const GeometryState& p) : GeometryState(p) {}
+
+  // Called at every surface intersection within the model
+  virtual void on_intersection() = 0;
+
+  /*
+   * Traces the ray through the geometry, calling on_intersection
+   * at every surface boundary.
+   */
+  void trace();
+
+  // Stops the ray and exits tracing when called from on_intersection
+  void stop() { stop_ = true; }
+
+  // Sets the dist_ variable
+  void compute_distance();
+
+protected:
+  // Records how far the ray has traveled
+  double traversal_distance_ {0.0};
+
+private:
+  // Max intersections before we assume ray tracing is caught in an infinite
+  // loop:
+  static const int MAX_INTERSECTIONS = 1000000;
+
+  bool stop_ {false};
+
+  unsigned event_counter_ {0};
+};
+
+} // namespace openmc
+#endif // OPENMC_RAY_H
diff --git a/src/plot.cpp b/src/plot.cpp
index 5e3342dd685..4d48beaf97c 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -1649,165 +1649,6 @@ void SolidRayTracePlot::set_diffuse_fraction(pugi::xml_node node)
   }
 }
 
-void Ray::compute_distance()
-{
-  boundary() = distance_to_boundary(*this);
-}
-
-void Ray::trace()
-{
-  // To trace the ray from its origin all the way through the model, we have
-  // to proceed in two phases. In the first, the ray may or may not be found
-  // inside the model. If the ray is already in the model, phase one can be
-  // skipped. Otherwise, the ray has to be advanced to the boundary of the
-  // model where all the cells are defined. Importantly, this is assuming that
-  // the model is convex, which is a very reasonable assumption for any
-  // radiation transport model.
-  //
-  // After phase one is done, we can starting tracing from cell to cell within
-  // the model. This step can use neighbor lists to accelerate the ray tracing.
-
-  bool inside_cell;
-  // Check for location if the particle is already known
-  if (lowest_coord().cell() == C_NONE) {
-    // The geometry position of the particle is either unknown or outside of the
-    // edge of the model.
-    if (lowest_coord().universe() == C_NONE) {
-      // Attempt to initialize the particle. We may have to
-      // enter a loop to move it up to the edge of the model.
-      inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
-    } else {
-      // It has been already calculated that the current position is outside of
-      // the edge of the model.
-      inside_cell = false;
-    }
-  } else {
-    // Availability of the cell means that the particle is located inside the
-    // edge.
-    inside_cell = true;
-  }
-
-  // Advance to the boundary of the model
-  while (!inside_cell) {
-    advance_to_boundary_from_void();
-    inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
-
-    // If true this means no surface was intersected. See cell.cpp and search
-    // for numeric_limits to see where we return it.
-    if (surface() == std::numeric_limits<int>::max()) {
-      warning(fmt::format("Lost a ray, r = {}, u = {}", r(), u()));
-      return;
-    }
-
-    // Exit this loop and enter into cell-to-cell ray tracing (which uses
-    // neighbor lists)
-    if (inside_cell)
-      break;
-
-    // if there is no intersection with the model, we're done
-    if (boundary().surface() == SURFACE_NONE)
-      return;
-
-    event_counter_++;
-    if (event_counter_ > MAX_INTERSECTIONS) {
-      warning("Likely infinite loop in ray traced plot");
-      return;
-    }
-  }
-
-  // Call the specialized logic for this type of ray. This is for the
-  // intersection for the first intersection if we had one.
-  if (boundary().surface() != SURFACE_NONE) {
-    // set the geometry state's surface attribute to be used for
-    // surface normal computation
-    surface() = boundary().surface();
-    on_intersection();
-    if (stop_)
-      return;
-  }
-
-  // reset surface attribute to zero after the first intersection so that it
-  // doesn't perturb surface crossing logic from here on out
-  surface() = 0;
-
-  // This is the ray tracing loop within the model. It exits after exiting
-  // the model, which is equivalent to assuming that the model is convex.
-  // It would be nice to factor out the on_intersection at the end of this
-  // loop and then do "while (inside_cell)", but we can't guarantee it's
-  // on a surface in that case. There might be some other way to set it
-  // up that is perhaps a little more elegant, but this is what works just
-  // fine.
-  while (true) {
-
-    compute_distance();
-
-    // There are no more intersections to process
-    // if we hit the edge of the model, so stop
-    // the particle in that case. Also, just exit
-    // if a negative distance was somehow computed.
-    if (boundary().distance() == INFTY || boundary().distance() == INFINITY ||
-        boundary().distance() < 0) {
-      return;
-    }
-
-    // See below comment where call_on_intersection is checked in an
-    // if statement for an explanation of this.
-    bool call_on_intersection {true};
-    if (boundary().distance() < 10 * TINY_BIT) {
-      call_on_intersection = false;
-    }
-
-    // DAGMC surfaces expect us to go a little bit further than the advance
-    // distance to properly check cell inclusion.
-    boundary().distance() += TINY_BIT;
-
-    // Advance particle, prepare for next intersection
-    for (int lev = 0; lev < n_coord(); ++lev) {
-      coord(lev).r() += boundary().distance() * coord(lev).u();
-    }
-    surface() = boundary().surface();
-    // Initialize last cells from the current cell, because the cell() variable
-    // does not contain the data for the case of a single-segment ray
-    for (int j = 0; j < n_coord(); ++j) {
-      cell_last(j) = coord(j).cell();
-    }
-    n_coord_last() = n_coord();
-    n_coord() = boundary().coord_level();
-    if (boundary().lattice_translation()[0] != 0 ||
-        boundary().lattice_translation()[1] != 0 ||
-        boundary().lattice_translation()[2] != 0) {
-      cross_lattice(*this, boundary(), settings::verbosity >= 10);
-    }
-
-    // Record how far the ray has traveled
-    traversal_distance_ += boundary().distance();
-    inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
-
-    // Call the specialized logic for this type of ray. Note that we do not
-    // call this if the advance distance is very small. Unfortunately, it seems
-    // darn near impossible to get the particle advanced to the model boundary
-    // and through it without sometimes accidentally calling on_intersection
-    // twice. This incorrectly shades the region as occluded when it might not
-    // actually be. By screening out intersection distances smaller than a
-    // threshold 10x larger than the scoot distance used to advance up to the
-    // model boundary, we can avoid that situation.
-    if (call_on_intersection) {
-      on_intersection();
-      if (stop_)
-        return;
-    }
-
-    if (!inside_cell)
-      return;
-
-    event_counter_++;
-    if (event_counter_ > MAX_INTERSECTIONS) {
-      warning("Likely infinite loop in ray traced plot");
-      return;
-    }
-  }
-}
-
 void ProjectionRay::on_intersection()
 {
   // This records a tuple with the following info
diff --git a/src/ray.cpp b/src/ray.cpp
new file mode 100644
index 00000000000..3d848e3a3a7
--- /dev/null
+++ b/src/ray.cpp
@@ -0,0 +1,168 @@
+#include "openmc/ray.h"
+
+#include "openmc/error.h"
+#include "openmc/geometry.h"
+#include "openmc/settings.h"
+
+namespace openmc {
+
+void Ray::compute_distance()
+{
+  boundary() = distance_to_boundary(*this);
+}
+
+void Ray::trace()
+{
+  // To trace the ray from its origin all the way through the model, we have
+  // to proceed in two phases. In the first, the ray may or may not be found
+  // inside the model. If the ray is already in the model, phase one can be
+  // skipped. Otherwise, the ray has to be advanced to the boundary of the
+  // model where all the cells are defined. Importantly, this is assuming that
+  // the model is convex, which is a very reasonable assumption for any
+  // radiation transport model.
+  //
+  // After phase one is done, we can starting tracing from cell to cell within
+  // the model. This step can use neighbor lists to accelerate the ray tracing.
+
+  bool inside_cell;
+  // Check for location if the particle is already known
+  if (lowest_coord().cell() == C_NONE) {
+    // The geometry position of the particle is either unknown or outside of the
+    // edge of the model.
+    if (lowest_coord().universe() == C_NONE) {
+      // Attempt to initialize the particle. We may have to
+      // enter a loop to move it up to the edge of the model.
+      inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
+    } else {
+      // It has been already calculated that the current position is outside of
+      // the edge of the model.
+      inside_cell = false;
+    }
+  } else {
+    // Availability of the cell means that the particle is located inside the
+    // edge.
+    inside_cell = true;
+  }
+
+  // Advance to the boundary of the model
+  while (!inside_cell) {
+    advance_to_boundary_from_void();
+    inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
+
+    // If true this means no surface was intersected. See cell.cpp and search
+    // for numeric_limits to see where we return it.
+    if (surface() == std::numeric_limits<int>::max()) {
+      warning(fmt::format("Lost a ray, r = {}, u = {}", r(), u()));
+      return;
+    }
+
+    // Exit this loop and enter into cell-to-cell ray tracing (which uses
+    // neighbor lists)
+    if (inside_cell)
+      break;
+
+    // if there is no intersection with the model, we're done
+    if (boundary().surface() == SURFACE_NONE)
+      return;
+
+    event_counter_++;
+    if (event_counter_ > MAX_INTERSECTIONS) {
+      warning("Likely infinite loop in ray traced plot");
+      return;
+    }
+  }
+
+  // Call the specialized logic for this type of ray. This is for the
+  // intersection for the first intersection if we had one.
+  if (boundary().surface() != SURFACE_NONE) {
+    // set the geometry state's surface attribute to be used for
+    // surface normal computation
+    surface() = boundary().surface();
+    on_intersection();
+    if (stop_)
+      return;
+  }
+
+  // reset surface attribute to zero after the first intersection so that it
+  // doesn't perturb surface crossing logic from here on out
+  surface() = 0;
+
+  // This is the ray tracing loop within the model. It exits after exiting
+  // the model, which is equivalent to assuming that the model is convex.
+  // It would be nice to factor out the on_intersection at the end of this
+  // loop and then do "while (inside_cell)", but we can't guarantee it's
+  // on a surface in that case. There might be some other way to set it
+  // up that is perhaps a little more elegant, but this is what works just
+  // fine.
+  while (true) {
+
+    compute_distance();
+
+    // There are no more intersections to process
+    // if we hit the edge of the model, so stop
+    // the particle in that case. Also, just exit
+    // if a negative distance was somehow computed.
+    if (boundary().distance() == INFTY || boundary().distance() == INFINITY ||
+        boundary().distance() < 0) {
+      return;
+    }
+
+    // See below comment where call_on_intersection is checked in an
+    // if statement for an explanation of this.
+    bool call_on_intersection {true};
+    if (boundary().distance() < 10 * TINY_BIT) {
+      call_on_intersection = false;
+    }
+
+    // DAGMC surfaces expect us to go a little bit further than the advance
+    // distance to properly check cell inclusion.
+    boundary().distance() += TINY_BIT;
+
+    // Advance particle, prepare for next intersection
+    for (int lev = 0; lev < n_coord(); ++lev) {
+      coord(lev).r() += boundary().distance() * coord(lev).u();
+    }
+    surface() = boundary().surface();
+    // Initialize last cells from the current cell, because the cell() variable
+    // does not contain the data for the case of a single-segment ray
+    for (int j = 0; j < n_coord(); ++j) {
+      cell_last(j) = coord(j).cell();
+    }
+    n_coord_last() = n_coord();
+    n_coord() = boundary().coord_level();
+    if (boundary().lattice_translation()[0] != 0 ||
+        boundary().lattice_translation()[1] != 0 ||
+        boundary().lattice_translation()[2] != 0) {
+      cross_lattice(*this, boundary(), settings::verbosity >= 10);
+    }
+
+    // Record how far the ray has traveled
+    traversal_distance_ += boundary().distance();
+    inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
+
+    // Call the specialized logic for this type of ray. Note that we do not
+    // call this if the advance distance is very small. Unfortunately, it seems
+    // darn near impossible to get the particle advanced to the model boundary
+    // and through it without sometimes accidentally calling on_intersection
+    // twice. This incorrectly shades the region as occluded when it might not
+    // actually be. By screening out intersection distances smaller than a
+    // threshold 10x larger than the scoot distance used to advance up to the
+    // model boundary, we can avoid that situation.
+    if (call_on_intersection) {
+      on_intersection();
+      if (stop_)
+        return;
+    }
+
+    if (!inside_cell)
+      return;
+
+    event_counter_++;
+    if (event_counter_ > MAX_INTERSECTIONS) {
+      warning("Likely infinite loop in ray traced plot");
+      return;
+    }
+  }
+}
+
+} // namespace openmc