diff --git a/CMakeLists.txt b/CMakeLists.txt
index 54f352f55b..34d71bc743 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -412,6 +412,7 @@ add_library(libopenmc SHARED
   src/position.cpp
   src/pugixml/pugixml_c.cpp
   src/random_lcg.cpp
+  src/secondary_nbody.cpp
   src/secondary_uncorrelated.cpp
   src/scattdata.cpp
   src/settings.cpp
diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp
new file mode 100644
index 0000000000..c843f8fc70
--- /dev/null
+++ b/src/secondary_nbody.cpp
@@ -0,0 +1,65 @@
+#include "secondary_nbody.h"
+
+#include <cmath> // for log
+
+#include "constants.h"
+#include "hdf5_interface.h"
+#include "math_functions.h"
+#include "random_lcg.h"
+
+namespace openmc {
+
+//==============================================================================
+// NBodyPhaseSpace implementation
+//==============================================================================
+
+NBodyPhaseSpace::NBodyPhaseSpace(hid_t group)
+{
+  read_attribute(group, "n_particles", n_bodies_);
+  read_attribute(group, "total_mass", mass_ratio_);
+  read_attribute(group, "atomic_weight_ratio", A_);
+  read_attribute(group, "q_value", Q_);
+}
+
+void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu) const
+{
+  // By definition, the distribution of the angle is isotropic for an N-body
+  // phase space distribution
+  mu = 2.0*prn() - 1.0;
+
+  // Determine E_max parameter
+  double Ap = mass_ratio_;
+  double E_max = (Ap - 1.0)/Ap * (A_/(A_ + 1.0)*E_in + Q_);
+
+  // x is essentially a Maxwellian distribution
+  double x = maxwell_spectrum_c(1.0);
+
+  double y;
+  double r1, r2, r3, r4, r5, r6;
+  switch (n_bodies_) {
+  case 3:
+    y = maxwell_spectrum_c(1.0);
+    break;
+  case 4:
+    r1 = prn();
+    r2 = prn();
+    r3 = prn();
+    y = -std::log(r1*r2*r3);
+    break;
+  case 5:
+    r1 = prn();
+    r2 = prn();
+    r3 = prn();
+    r4 = prn();
+    r5 = prn();
+    r6 = prn();
+    y = -std::log(r1*r2*r3*r4) - std::log(r5) * std::pow(std::cos(PI/2.0*r6), 2);
+    break;
+  }
+
+  // Now determine v and E_out
+  double v = x/(x + y);
+  E_out = E_max * v;
+}
+
+}
\ No newline at end of file
diff --git a/src/secondary_nbody.h b/src/secondary_nbody.h
new file mode 100644
index 0000000000..a739dde96d
--- /dev/null
+++ b/src/secondary_nbody.h
@@ -0,0 +1,23 @@
+#ifndef OPENMC_SECONDARY_NBODY_H
+#define OPENMC_SECONDARY_NBODY_H
+
+#include "hdf5.h"
+
+#include "angle_energy.h"
+
+namespace openmc {
+
+class NBodyPhaseSpace : public AngleEnergy {
+public:
+  explicit NBodyPhaseSpace(hid_t group);
+  void sample(double E_in, double& E_out, double& mu) const;
+private:
+  int n_bodies_;
+  double mass_ratio_;
+  double A_;
+  double Q_;
+};
+
+}
+
+#endif // OPENMC_SECONDARY_NBODY_H