wip

Author: GuySten

src/angle_energy.cpp

@@ -1,7 +1,9 @@
 #include "openmc/angle_energy.h"
 
-#include <cmath> // for sqrt
+#include <algorithm> // for clamp
+#include <cmath>     // for sqrt
 
+#include "openmc/constants.h"
 #include "openmc/random_lcg.h"
 
 namespace openmc {
@@ -23,9 +25,16 @@ double get_jac_and_transform_impl(
     return 0.0;
   D = std::sqrt(D);
 
-  if (mu_lab <= 0.0 && E_cm <= E_cm)
+  if ((mu_lab <= 0.0) && (E_cm <= E_com))
     return 0.0;
-  double E_out1 = E_com * (mu_lab + D) * (mu_lab + D);
+
+  double E_out1;
+  if (mu_lab > 0.0) {
+    E_out1 = E_com * (mu_lab + D) * (mu_lab + D);
+  } else {
+    E_out1 = E_com * ((E_cm / E_com - 1.0) / (D - mu_lab)) *
+             ((E_cm / E_com - 1.0) / (D - mu_lab));
+  }
   double mult;
   if (E_cm > E_com) {
     mult = 1.0;
@@ -39,6 +48,9 @@ double get_jac_and_transform_impl(
     }
   }
   mu = mu_lab * std::sqrt(E_out / E_cm) - std::sqrt(E_com / E_cm);
+
+  if (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));
 }
 

src/tallies/next_event_scoring.cpp

@@ -1,5 +1,7 @@
 #include "openmc/tallies/next_event_scoring.h"
 
+#include <algorithm> // for clamp
+
 #include "openmc/secondary_uncorrelated.h"
 #include "openmc/source.h"
 
@@ -14,11 +16,13 @@ void score_point_tally_elastic(
 
   // Neutron velocity in LAB
   Direction v_n = std::sqrt(p.E()) * p.u();
+  auto u_n = v_n / v_n.norm();
 
   // Velocity of center-of-mass
   Direction v_cm = (v_n + awr * v_t) / (awr + 1.0);
   auto u_cm = v_cm / v_cm.norm();
 
+  double E_in = p.E();
   double E_com = v_cm.dot(v_cm);
   double E_out = (v_n - v_cm).dot(v_n - v_cm);
 
@@ -27,11 +31,15 @@ void score_point_tally_elastic(
 
   auto pdf = [&](Direction u, double& E) {
     double mu = u.dot(u_cm);
+    double mu_l = u.dot(u_n);
     E = E_out;
     double jac =
       get_jac_and_transform_impl(E_com, mu, E, p.current_seed(), awr);
+    double mu_cm =
+      1.0 + mu_l * std::sqrt(E_in * E) / E_out - (E_in + E) / (2.0 * E_out);
+    mu_cm = std::clamp(mu_cm, -1.0, 1.0);
     if (!d_->angle().empty()) {
-      return jac * d_->angle().evaluate(p.E(), mu) / (2.0 * PI);
+      return jac * d_->angle().evaluate(p.E(), mu_cm) / (2.0 * PI);
     } else {
       return jac * 0.5 / (2.0 * PI);
     }