From 8a07feebcb7e92d7ab346a9d4510354f7ae87ec4 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 24 Mar 2021 08:11:40 -0500
Subject: [PATCH] Sample isotropic direction consistently

---
 include/openmc/distribution_multi.h |  2 ++
 src/distribution_multi.cpp          |  7 ++++++-
 src/photon.cpp                      | 15 +++------------
 src/physics.cpp                     | 24 ++++++------------------
 4 files changed, 17 insertions(+), 31 deletions(-)

diff --git a/include/openmc/distribution_multi.h b/include/openmc/distribution_multi.h
index ae7768ac3..ddefa422c 100644
--- a/include/openmc/distribution_multi.h
+++ b/include/openmc/distribution_multi.h
@@ -57,6 +57,8 @@ private:
 //! Uniform distribution on the unit sphere
 //==============================================================================
 
+Direction isotropic_direction(uint64_t* seed);
+
 class Isotropic : public UnitSphereDistribution {
 public:
   Isotropic() { };
diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp
index 61d157024..511217e5f 100644
--- a/src/distribution_multi.cpp
+++ b/src/distribution_multi.cpp
@@ -69,7 +69,7 @@ Direction PolarAzimuthal::sample(uint64_t* seed) const
 // Isotropic implementation
 //==============================================================================
 
-Direction Isotropic::sample(uint64_t* seed) const
+Direction isotropic_direction(uint64_t* seed)
 {
   double phi = 2.0*PI*prn(seed);
   double mu = 2.0*prn(seed) - 1.0;
@@ -77,6 +77,11 @@ Direction Isotropic::sample(uint64_t* seed) const
       std::sqrt(1.0 - mu*mu) * std::sin(phi)};
 }
 
+Direction Isotropic::sample(uint64_t* seed) const
+{
+  return isotropic_direction(seed);
+}
+
 //==============================================================================
 // Monodirectional implementation
 //==============================================================================
diff --git a/src/photon.cpp b/src/photon.cpp
index 6333f9b29..ee89fa84d 100644
--- a/src/photon.cpp
+++ b/src/photon.cpp
@@ -2,6 +2,7 @@
 
 #include "openmc/bremsstrahlung.h"
 #include "openmc/constants.h"
+#include "openmc/distribution_multi.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/particle.h"
@@ -662,12 +663,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
 {
   // If no transitions, assume fluorescent photon from captured free electron
   if (shell.n_transitions == 0) {
-    double mu = 2.0*prn(p.current_seed()) - 1.0;
-    double phi = 2.0*PI*prn(p.current_seed());
-    Direction u;
-    u.x = mu;
-    u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
-    u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
+    Direction u = isotropic_direction(p.current_seed());
     double E = shell.binding_energy;
     p.create_secondary(p.wgt_, u, E, Particle::Type::photon);
     return;
@@ -687,12 +683,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
   int secondary = shell.transition_subshells(i_transition, 1);
 
   // Sample angle isotropically
-  double mu = 2.0*prn(p.current_seed()) - 1.0;
-  double phi = 2.0*PI*prn(p.current_seed());
-  Direction u;
-  u.x = mu;
-  u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
-  u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
+  Direction u = isotropic_direction(p.current_seed());
 
   // Get the transition energy
   double E = shell.transition_energy(i_transition);
diff --git a/src/physics.cpp b/src/physics.cpp
index a0be7f5c4..8e4b027ba 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -3,6 +3,7 @@
 #include "openmc/bank.h"
 #include "openmc/bremsstrahlung.h"
 #include "openmc/constants.h"
+#include "openmc/distribution_multi.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/endf.h"
 #include "openmc/error.h"
@@ -425,12 +426,7 @@ void sample_positron_reaction(Particle& p)
   }
 
   // Sample angle isotropically
-  double mu = 2.0*prn(p.current_seed()) - 1.0;
-  double phi = 2.0*PI*prn(p.current_seed());
-  Direction u;
-  u.x = mu;
-  u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
-  u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
+  Direction u = isotropic_direction(p.current_seed());
 
   // Create annihilation photon pair traveling in opposite directions
   p.create_secondary(p.wgt_, u, MASS_ELECTRON_EV, Particle::Type::photon);
@@ -712,13 +708,7 @@ void scatter(Particle& p, int i_nuclide)
     int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
     if (mat->p0_[i_nuc_mat]) {
       // Sample isotropic-in-lab outgoing direction
-      double mu = 2.0*prn(p.current_seed()) - 1.0;
-      double phi = 2.0*PI*prn(p.current_seed());
-
-      // Change direction of particle
-      p.u().x = mu;
-      p.u().y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
-      p.u().z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
+      p.u() = isotropic_direction(p.current_seed());
       p.mu_ = u_old.dot(p.u());
     }
   }
@@ -1068,11 +1058,9 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part
   }
 
   // Sample azimuthal angle uniformly in [0, 2*pi) and assign angle
-  // TODO: use rotate_angle instead
-  double phi = 2.0*PI*prn(seed);
-  site->u.x = mu;
-  site->u.y = std::sqrt(1.0 - mu*mu) * std::cos(phi);
-  site->u.z = std::sqrt(1.0 - mu*mu) * std::sin(phi);
+  // TODO: account for dependence on incident neutron?
+  Direction ref(1., 0., 0.);
+  site->u = rotate_angle(ref, mu, nullptr, seed);
 }
 
 void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)