From 79ef612cb5769ce2e3dee360dc9b38508741ba40 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 21 Dec 2018 14:30:28 -0600
Subject: [PATCH] Fix multipole and two other bugs

---
 include/openmc/nuclide.h |  2 ++
 src/nuclide.cpp          |  9 ++++-----
 src/nuclide_header.F90   | 15 +++++++++++++++
 3 files changed, 21 insertions(+), 5 deletions(-)

diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h
index 6deb6c01d..a8178b11c 100644
--- a/include/openmc/nuclide.h
+++ b/include/openmc/nuclide.h
@@ -193,6 +193,8 @@ extern "C" MaterialMacroXS material_xs;
 //==============================================================================
 
 extern "C" void set_micro_xs();
+extern "C" void nuclide_multipole_eval(int i_nuclide, double E, double sqrtkT,
+  double* sig_s, double* sig_a, double* sig_f);
 extern "C" bool nuclide_wmp_present(int i_nuclide);
 extern "C" double nuclide_wmp_emin(int i_nuclide);
 extern "C" double nuclide_wmp_emax(int i_nuclide);
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index 753c0af67..53579cc85 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -520,7 +520,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
   if (use_mp) {
     // Call multipole kernel
     double sig_s, sig_a, sig_f;
-    //multipole_eval(this % multipole, E, sqrtkT, sig_s, sig_a, sig_f)
+    nuclide_multipole_eval(i_nuclide_, E, sqrtkT, &sig_s, &sig_a, &sig_f);
 
     micro_xs.total = sig_s + sig_a;
     micro_xs.elastic = sig_s;
@@ -651,11 +651,10 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
         // If reaction is present and energy is greater than threshold, set the
         // reaction xs appropriately
         int i_rx = reaction_index_[DEPLETION_RX[j]];
-
-        const auto& rx = reactions_[i_rx];
-        const auto& rx_xs = rx->xs_[i_temp].value;
-
         if (i_rx >= 0) {
+          const auto& rx = reactions_[i_rx];
+          const auto& rx_xs = rx->xs_[i_temp].value;
+
           // Physics says that (n,gamma) is not a threshold reaction, so we don't
           // need to specifically check its threshold index
           if (j == 0) {
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 319ea79a7..cc6f1761c 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -1048,6 +1048,9 @@ contains
         ! Add entry to nuclide dictionary
         call nuclide_dict % set(to_lower(name_), n)
         n_nuclides = n
+
+        ! Initialize nuclide grid
+        call nuclides(n) % init_grid()
       else
         err = E_DATA
         call set_errmsg("Nuclide '" // trim(name_) // "' is not present &
@@ -1100,4 +1103,16 @@ contains
     E = nuclides(i_nuclide + 1) % multipole % E_max
   end function
 
+  subroutine nuclide_multipole_eval(i_nuclide, E, sqrtkT, sig_s, sig_a, sig_f) bind(C)
+    integer(C_INT), value :: i_nuclide
+    real(C_DOUBLE), value :: E
+    real(C_DOUBLE), value :: sqrtkT
+    real(C_DOUBLE), intent(out) :: sig_s
+    real(C_DOUBLE), intent(out) :: sig_a
+    real(C_DOUBLE), intent(out) :: sig_f
+
+    call multipole_eval(nuclides(i_nuclide + 1) % multipole, E, sqrtkT, &
+         sig_s, sig_a, sig_f)
+  end subroutine
+
 end module nuclide_header