From 714abb2c392c06106048dd94f82a320648a35129 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 9 Jul 2018 23:02:27 -0500
Subject: [PATCH] Convert UncorrelatedAngleEnergy

---
 CMakeLists.txt                 |  1 +
 src/angle_energy.h             | 14 +++++++
 src/distribution_angle.h       |  4 +-
 src/secondary_uncorrelated.cpp | 74 ++++++++++++++++++++++++++++++++++
 src/secondary_uncorrelated.h   | 26 ++++++++++++
 5 files changed, 118 insertions(+), 1 deletion(-)
 create mode 100644 src/angle_energy.h
 create mode 100644 src/secondary_uncorrelated.cpp
 create mode 100644 src/secondary_uncorrelated.h

diff --git a/CMakeLists.txt b/CMakeLists.txt
index db3f3f6f5..54f352f55 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_uncorrelated.cpp
   src/scattdata.cpp
   src/settings.cpp
   src/simulation.cpp
diff --git a/src/angle_energy.h b/src/angle_energy.h
new file mode 100644
index 000000000..005c62018
--- /dev/null
+++ b/src/angle_energy.h
@@ -0,0 +1,14 @@
+#ifndef OPENMC_ANGLE_ENERGY_H
+#define OPENMC_ANGLE_ENERGY_H
+
+namespace openmc {
+
+class AngleEnergy {
+public:
+  virtual void sample(double E_in, double& E_out, double& mu) const = 0;
+  virtual ~AngleEnergy() = default;
+};
+
+}
+
+#endif // OPENMC_ANGLE_ENERGY_H
diff --git a/src/distribution_angle.h b/src/distribution_angle.h
index 3e7dd51cc..3a2f4b62a 100644
--- a/src/distribution_angle.h
+++ b/src/distribution_angle.h
@@ -1,7 +1,6 @@
 #ifndef OPENMC_DISTRIBUTION_ANGLE_H
 #define OPENMC_DISTRIBUTION_ANGLE_H
 
-#include <memory> // for unique_ptr
 #include <vector> // for vector
 
 #include "distribution.h"
@@ -11,8 +10,11 @@ namespace openmc {
 
 class AngleDistribution {
 public:
+  AngleDistribution() = default;
   explicit AngleDistribution(hid_t group);
   double sample(double E) const;
+
+  bool empty() const { return energy_.empty(); }
 private:
   std::vector<double> energy_;
   std::vector<UPtrDist> distribution_;
diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp
new file mode 100644
index 000000000..99b68bad5
--- /dev/null
+++ b/src/secondary_uncorrelated.cpp
@@ -0,0 +1,74 @@
+#include "secondary_uncorrelated.h"
+
+#include <sstream> // for stringstream
+#include <string>  // for string
+
+#include "error.h"
+#include "hdf5_interface.h"
+#include "random_lcg.h"
+
+namespace openmc {
+
+//==============================================================================
+// UncorrelatedAngleEnergy implementation
+//==============================================================================
+
+UncorrelatedAngleEnergy::UncorrelatedAngleEnergy(hid_t group)
+{
+  // Check if angle group is present & read
+  if (object_exists(group, "angle")) {
+    hid_t angle_group = open_group(group, "angle");
+    angle_ = AngleDistribution{angle_group};
+    close_group(angle_group);
+  }
+
+  // Check if energy group is present & read
+  if (object_exists(group, "energy")) {
+    hid_t energy_group = open_group(group, "energy");
+
+    std::string type;
+    read_attribute(energy_group, "type", type);
+    using UPtrEDist = std::unique_ptr<EnergyDistribution>;
+    if (type == "discrete_photon") {
+      energy_ = UPtrEDist{new DiscretePhoton{energy_group}};
+    } else if (type == "level") {
+      energy_ = UPtrEDist{new LevelInelastic{energy_group}};
+    } else if (type == "continuous") {
+      energy_ = UPtrEDist{new ContinuousTabular{energy_group}};
+    } else if (type == "maxwell") {
+      energy_ = UPtrEDist{new MaxwellEnergy{energy_group}};
+    } else if (type == "evaporation") {
+      energy_ = UPtrEDist{new Evaporation{energy_group}};
+    } else if (type == "watt") {
+      energy_ = UPtrEDist{new WattEnergy{energy_group}};
+    } else {
+      std::stringstream msg;
+      msg << "Energy distribution type '" << type << "' not implemented.";
+      warning(msg);
+    }
+    close_group(energy_group);
+  }
+
+}
+
+void
+UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu) const
+{
+  // Sample cosine of scattering angle
+  if (fission_) {
+    // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+    // For fission, the angle is not used, so just assign a dummy value
+    mu = 1.0;
+    // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+  } else if (!angle_.empty()) {
+    mu = angle_.sample(E_in);
+  } else {
+    // no angle distribution given => assume isotropic for all energies
+    mu = 2.0*prn() - 1.0;
+  }
+
+  // Sample outgoing energy
+  E_out = energy_->sample(E_in);
+}
+
+}
diff --git a/src/secondary_uncorrelated.h b/src/secondary_uncorrelated.h
new file mode 100644
index 000000000..73c62c14a
--- /dev/null
+++ b/src/secondary_uncorrelated.h
@@ -0,0 +1,26 @@
+#ifndef OPENMC_SECONDARY_UNCORRELATED_H
+#define OPENMC_SECONDARY_UNCORRELATED_H
+
+#include <memory>
+#include <vector>
+
+#include "hdf5.h"
+#include "angle_energy.h"
+#include "distribution_angle.h"
+#include "distribution_energy.h"
+
+namespace openmc {
+
+class UncorrelatedAngleEnergy : public AngleEnergy {
+public:
+  explicit UncorrelatedAngleEnergy(hid_t group);
+  void sample(double E_in, double& E_out, double& mu) const;
+private:
+  bool fission_ {false};
+  AngleDistribution angle_;
+  std::unique_ptr<EnergyDistribution> energy_;
+};
+
+}
+
+#endif // OPENMC_SECONDARY_UNCORRELATED_H