From 7ee8cc267b3fc1135a1df2346f1136299e5113be Mon Sep 17 00:00:00 2001
From: GuySten <guyste@post.bgu.ac.il>
Date: Mon, 11 May 2026 11:02:37 +0300
Subject: [PATCH] update

---
 CMakeLists.txt                |  1 +
 include/openmc/angle_energy.h | 44 ++++--------------------------
 src/angle_energy.cpp          | 50 +++++++++++++++++++++++++++++++++++
 3 files changed, 56 insertions(+), 39 deletions(-)
 create mode 100644 src/angle_energy.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9fe133a22e..736d13dab6 100644
--- a/CMakeLists.txt
+++ b/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
diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h
index 1858d7a3b1..1ac8f6feb8 100644
--- a/include/openmc/angle_energy.h
+++ b/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
 
diff --git a/src/angle_energy.cpp b/src/angle_energy.cpp
new file mode 100644
index 0000000000..a3c9c83af7
--- /dev/null
+++ b/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