update

CMakeLists.txt

@@ -355,6 +355,7 @@ endif()
 #===============================================================================
 
 list(APPEND libopenmc_SOURCES
+  src/angle_energy.cpp
   src/atomic_mass.cpp
   src/bank.cpp
   src/boundary_condition.cpp

include/openmc/angle_energy.h

@@ -1,11 +1,8 @@
 #ifndef OPENMC_ANGLE_ENERGY_H
 #define OPENMC_ANGLE_ENERGY_H
 
-#include <cmath> // for sqrt
 #include <cstdint>
 
-#include "openmc/random_lcg.h"
-
 namespace openmc {
 
 //==============================================================================
@@ -38,42 +35,11 @@ public:
   virtual ~AngleEnergy() = default;
 };
 
-inline double get_jac_and_transform(
-  double E_in, double& mu, double& E_out, uint64_t* seed, double awr)
-{
-  double E_cm = E_out;
-  double mu_lab = mu;
-  double D = mu_lab * mu_lab - 1.0 + (awr + 1.0) * (awr + 1.0) * E_cm / E_in;
-  if (D < 0.0)
-    return 0.0;
-  D = std::sqrt(D);
-  double E_out1 =
-    E_in * ((mu_lab + D) / (awr + 1.0)) * ((mu_lab + D) / (awr + 1.0));
-  double E_out2 =
-    E_in * ((mu_lab - D) / (awr + 1.0)) * ((mu_lab - D) / (awr + 1.0));
-  double mu_cm1 =
-    mu_lab * std::sqrt(E_out1 / E_cm) - std::sqrt(E_in / E_cm) / (awr + 1.0);
-  double mu_cm2 =
-    mu_lab * std::sqrt(E_out1 / E_cm) - std::sqrt(E_in / E_cm) / (awr + 1.0);
-  double mult = 1.0;
-  if (std::abs(mu_cm1) > 1.0) {
-    mu = mu_cm2;
-    E_out = E_out2;
-  } else if (std::abs(mu_cm2) > 1.0) {
-    mu = mu_cm1;
-    E_out = E_out1;
-  } else {
-    mult = 2.0;
-    if (prn(seed) < 0.5) {
-      mu = mu_cm2;
-      E_out = E_out2;
-    } else {
-      mu = mu_cm2;
-      E_out = E_out2;
-    }
-  }
-  return mult * E_out * (awr + 1.0) / (D * std::sqrt(E_cm * E_in));
-}
+double get_jac_and_transform(
+  double E_in, double& mu, double& E_out, uint64_t* seed, double awr);
+
+double get_jac_and_transform_impl(
+  double E_com, double& mu, double& E_out, uint64_t* seed, double awr);
 
 } // namespace openmc
 

src/angle_energy.cpp

@@ -0,0 +1,50 @@
+#include "openmc/angle_energy.h"
+
+#include <algorithm> // for clamp
+#include <cmath>     // for sqrt
+
+#include "openmc/constants.h"
+#include "openmc/random_lcg.h"
+
+namespace openmc {
+
+double get_jac_and_transform(
+  double E_in, double& mu, double& E_out, uint64_t* seed, double awr)
+{
+  double E_com = E_in / ((awr + 1.0) * (awr + 1.0));
+  return get_jac_and_transform_impl(E_com, mu, E_out, seed, awr);
+}
+
+double get_jac_and_transform_impl(
+  double E_com, double& mu, double& E_out, uint64_t* seed, double awr)
+{
+  double E_cm = E_out;
+  double mu_lab = mu;
+  double D = mu_lab * mu_lab - 1.0 + E_cm / E_com;
+  if (D <= 0.0)
+    return 0.0;
+  D = std::sqrt(D);
+
+  if ((mu_lab <= 0.0) && (E_cm <= E_com))
+    return 0.0;
+  double E_out1 = E_com * (mu_lab + D) * (mu_lab + D);
+  double mult;
+  if (E_cm > E_com) {
+    mult = 1.0;
+    E_out = E_out1;
+  } else {
+    mult = 2.0;
+    if (prn(seed) < 0.5) {
+      E_out = E_out1;
+    } else {
+      E_out = E_com * (mu_lab - D) * (mu_lab - D);
+    }
+  }
+  mu = mu_lab * std::sqrt(E_out / E_cm) - std::sqrt(E_com / E_cm);
+
+  if ((std::abs(mu) > 1.0) && (std::abs(mu) < 1.0 + FP_PRECISION))
+    mu = std::clamp(mu, -1.0, 1.0);
+  return mult * E_out / (D * std::sqrt(E_cm * E_com));
+}
+
+} // namespace openmc