From 8894696db13aa09159b3ebc6a793f6f3fe38b00d Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 8 Jan 2020 16:13:07 -0600
Subject: [PATCH] Scale fission photon production by (EGP + EGD)/EGP

---
 include/openmc/nuclide.h                      |  2 +-
 src/nuclide.cpp                               | 25 ++++++++++++++++---
 .../results_true.dat                          | 24 +++++++++---------
 3 files changed, 35 insertions(+), 16 deletions(-)

diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h
index a1571ea58e..dd9147c229 100644
--- a/include/openmc/nuclide.h
+++ b/include/openmc/nuclide.h
@@ -97,7 +97,7 @@ public:
   std::vector<int> index_inelastic_scatter_;
 
 private:
-  void create_derived();
+  void create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons);
 
   static int XS_TOTAL;
   static int XS_ABSORPTION;
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index 86f2ee1408..bffac7d2be 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -259,17 +259,24 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature, int i_nucl
   }
 
   // Read fission energy release data if present
+  std::unique_ptr<Function1D> prompt_photons;
+  std::unique_ptr<Function1D> delayed_photons;
   if (object_exists(group, "fission_energy_release")) {
     hid_t fer_group = open_group(group, "fission_energy_release");
     fission_q_prompt_ = read_function(fer_group, "q_prompt");
     fission_q_recov_ = read_function(fer_group, "q_recoverable");
+
+    // We need prompt/delayed photon energy release for scaling fission photon
+    // production
+    prompt_photons = read_function(fer_group, "prompt_photons");
+    delayed_photons = read_function(fer_group, "delayed_photons");
     close_group(fer_group);
   }
 
-  this->create_derived();
+  this->create_derived(prompt_photons.get(), delayed_photons.get());
 }
 
-void Nuclide::create_derived()
+void Nuclide::create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons)
 {
   for (const auto& grid : grid_) {
     // Allocate and initialize cross section
@@ -294,7 +301,19 @@ void Nuclide::create_derived()
           auto pprod = xt::view(xs_[t], xt::range(j, j+n), XS_PHOTON_PROD);
           for (int k = 0; k < n; ++k) {
             double E = grid_[t].energy[k+j];
-            pprod[k] += xs[k] * (*p.yield_)(E);
+
+            // For fission, artificially increase the photon yield to account
+            // for delayed photons
+            double f = 1.0;
+            if (is_fission(rx->mt_)) {
+              if (prompt_photons && delayed_photons) {
+                double energy_prompt = (*prompt_photons)(E);
+                double energy_delayed = (*delayed_photons)(E);
+                f = (energy_prompt + energy_delayed)/(energy_prompt);
+              }
+            }
+
+            pprod[k] += f * xs[k] * (*p.yield_)(E);
           }
         }
       }
diff --git a/tests/regression_tests/photon_production_fission/results_true.dat b/tests/regression_tests/photon_production_fission/results_true.dat
index 60e12bc3a0..ca4517b4ff 100644
--- a/tests/regression_tests/photon_production_fission/results_true.dat
+++ b/tests/regression_tests/photon_production_fission/results_true.dat
@@ -15,14 +15,14 @@ tally 1:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.491697E+07
-2.070245E+14
+4.764667E+07
+7.568912E+14
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.491697E+07
-2.070245E+14
+4.764667E+07
+7.568912E+14
 0.000000E+00
 0.000000E+00
 tally 2:
@@ -40,14 +40,14 @@ tally 2:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.508579E+07
-2.098299E+14
+4.749088E+07
+7.519602E+14
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.508579E+07
-2.098299E+14
+4.749088E+07
+7.519602E+14
 0.000000E+00
 0.000000E+00
 tally 3:
@@ -65,13 +65,13 @@ tally 3:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.495967E+07
-2.077339E+14
+4.729796E+07
+7.458273E+14
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-2.504536E+07
-2.091625E+14
+4.743827E+07
+7.502589E+14
 0.000000E+00
 0.000000E+00