From 2cecbd5d0f6b449fd4c6df79079a7f7c1f15a504 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <drewej@protonmail.com>
Date: Mon, 9 Sep 2019 17:11:45 -0500
Subject: [PATCH] Add Nuclide::modified_fission_q_ from Tabulated1D, Polynomial

The new Function1D attribute holds recoverable fission
energy minus energy from prompt neutrons. The motivation for this
is to treat the total energy deposition as energy from heating
and energy released from fission. There is overlap here, as
for fissile nuclides, the energy released from fission is included
in parts of the fission heating data include energy from the fission
event.

Nuclide::modified_fission_q_ currently includes the recoverable
energy from fission minus energy from prompt neutrons. This leaves
energy from prompt and delayed photons, fission fragments, betas,
and delayed neutrons.
---
 include/openmc/nuclide.h |  2 ++
 src/nuclide.cpp          | 48 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 50 insertions(+)

diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h
index a936ac0102..51d1277ad2 100644
--- a/include/openmc/nuclide.h
+++ b/include/openmc/nuclide.h
@@ -80,6 +80,8 @@ public:
   std::unique_ptr<Function1D> total_nu_; //!< Total neutron yield
   std::unique_ptr<Function1D> fission_q_prompt_; //!< Prompt fission energy release
   std::unique_ptr<Function1D> fission_q_recov_; //!< Recoverable fission energy release
+  //! Fission energy release without prompt neutrons
+  std::unique_ptr<Function1D> modified_fission_q_;
 
   // Resonance scattering information
   bool resonant_ {false};
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index 3b092ddf4b..1333ec0a7c 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -263,6 +263,54 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature, int i_nucl
     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");
+
+    // Check that fission q data is Tabulated1D or Polynomial
+    Tabulated1D* tab_recov = dynamic_cast<Tabulated1D*>(fission_q_recov_.get());
+    Polynomial* poly_recov = dynamic_cast<Polynomial*>(fission_q_recov_.get());
+    if (nullptr == tab_recov && nullptr == poly_recov) {
+      throw std::runtime_error{"Recoverable fission energy in " + object_name(fer_group)
+        + " is not Tabulated1D nor Polynomial."};
+    }
+
+    // Subtract away energy from prompt neutrons and photons
+    // needed for energy deposition tally, as neutron and photon
+    // heating are included in MT=301
+
+    auto prompt_neutrons = read_function(fer_group, "prompt_neutrons");
+    std::vector<double> mod_data;
+
+    // Check if fission_q and prompt_neutrons are tabulated 1D
+
+    if (tab_recov) {
+      Tabulated1D* neutrons = dynamic_cast<Tabulated1D*>(prompt_neutrons.get());
+      if (nullptr == neutrons) {
+        throw std::runtime_error{"prompt_neutrons not Tabulated1D in " + 
+          object_name(fer_group)};
+      }
+      auto recov_data = tab_recov->y();
+      auto recov_x = tab_recov->x();
+      int i=0;
+      for (auto yit = recov_data.cbegin(); yit != recov_data.cend(); ++yit) {
+        mod_data.push_back(*yit - (*neutrons)(recov_x[i]));
+        ++i;
+      }
+      modified_fission_q_ = std::make_unique<Tabulated1D>(tab_recov->x(), mod_data);
+    } else {
+        Polynomial* neutrons = dynamic_cast<Polynomial*>(prompt_neutrons.get());
+      if (nullptr == neutrons) {
+        throw std::runtime_error{"prompt_neutrons not Polynomial in " + 
+          object_name(fer_group)};
+      }
+      auto recov_data = poly_recov->coeffs();
+      auto neutron_data = neutrons->coeffs();
+      int i=0;
+      for (auto it = recov_data.cbegin(); it != recov_data.cend(); ++it) {
+        mod_data.push_back(*it - neutron_data.at(i));
+        ++i;
+      }
+      modified_fission_q_ = std::make_unique<Polynomial>(mod_data);
+    }
+    
     close_group(fer_group);
   }