Merge pull request #1016 from nelsonag/mgxs_cpp

Conversion of Mgxs code to C++

CMakeLists.txt

@@ -363,7 +363,7 @@ set(LIBOPENMC_FORTRAN_SRC
   src/mesh_header.F90
   src/message_passing.F90
   src/mgxs_data.F90
-  src/mgxs_header.F90
+  src/mgxs_interface.F90
   src/multipole_header.F90
   src/nuclide_header.F90
   src/output.F90
@@ -380,7 +380,6 @@ set(LIBOPENMC_FORTRAN_SRC
   src/reaction_header.F90
   src/relaxng
   src/sab_header.F90
-  src/scattdata_header.F90
   src/secondary_correlated.F90
   src/secondary_kalbach.F90
   src/secondary_nbody.F90
@@ -440,12 +439,17 @@ set(LIBOPENMC_CXX_SRC
   src/lattice.cpp
   src/math_functions.cpp
   src/message_passing.cpp
+  src/mgxs.cpp
+  src/mgxs_interface.cpp
   src/plot.cpp
   src/random_lcg.cpp
+  src/scattdata.cpp
   src/simulation.cpp
   src/state_point.cpp
+  src/string_functions.cpp
   src/surface.cpp
   src/xml_interface.cpp
+  src/xsdata.cpp
   src/pugixml/pugixml.cpp)
 add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC} ${LIBOPENMC_CXX_SRC})
 set_target_properties(libopenmc PROPERTIES

openmc/mesh.py

@@ -208,12 +208,12 @@ class Mesh(IDManagerMixin):
 
         shape = np.array(lattice.shape)
         width = lattice.pitch*shape
-        
+
         mesh = cls(mesh_id, name)
         mesh.lower_left = lattice.lower_left
         mesh.upper_right = lattice.lower_left + width
         mesh.dimension = shape*division
-        
+
         return mesh
 
     def to_xml_element(self):
@@ -333,14 +333,16 @@ class Mesh(IDManagerMixin):
         if n_dim == 1:
             universe_array = np.array([universes])
         elif n_dim == 2:
-            universe_array = np.empty(self.dimension, dtype=openmc.Universe)
+            universe_array = np.empty(self.dimension[::-1],
+                                      dtype=openmc.Universe)
             i = 0
             for y in range(self.dimension[1] - 1, -1, -1):
                 for x in range(self.dimension[0]):
                     universe_array[y][x] = universes[i]
                     i += 1
         else:
-            universe_array = np.empty(self.dimension, dtype=openmc.Universe)
+            universe_array = np.empty(self.dimension[::-1],
+                                      dtype=openmc.Universe)
             i = 0
             for z in range(self.dimension[2]):
                 for y in range(self.dimension[1] - 1, -1, -1):

src/api.F90

@@ -129,7 +129,7 @@ contains
     index_ufs_mesh = -1
     keff = ONE
     legendre_to_tabular = .true.
-    legendre_to_tabular_points = 33
+    legendre_to_tabular_points = C_NONE
     n_batch_interval = 1
     n_lost_particles = 0
     n_particles = 0
@@ -310,7 +310,6 @@ contains
   subroutine free_memory()
 
     use cmfd_header
-    use mgxs_header
     use plot_header
     use sab_header
     use settings
@@ -328,7 +327,6 @@ contains
     call free_memory_simulation()
     call free_memory_nuclide()
     call free_memory_settings()
-    call free_memory_mgxs()
     call free_memory_sab()
     call free_memory_source()
     call free_memory_mesh()

src/cmfd_input.F90

@@ -3,8 +3,8 @@ module cmfd_input
   use, intrinsic :: ISO_C_BINDING
 
   use cmfd_header
-  use mesh_header, only: mesh_dict
-  use mgxs_header, only: energy_bins
+  use mesh_header,    only: mesh_dict
+  use mgxs_interface, only: energy_bins, num_energy_groups
   use tally
   use tally_header
   use timer_header

src/constants.F90

@@ -405,6 +405,25 @@ module constants
        DIFF_NUCLIDE_DENSITY = 2, &
        DIFF_TEMPERATURE = 3
 
+
+  ! Mgxs::get_xs enumerated types
+  integer(C_INT), parameter :: &
+       MG_GET_XS_TOTAL              = 0, &
+       MG_GET_XS_ABSORPTION         = 1, &
+       MG_GET_XS_INVERSE_VELOCITY   = 2, &
+       MG_GET_XS_DECAY_RATE         = 3, &
+       MG_GET_XS_SCATTER            = 4, &
+       MG_GET_XS_SCATTER_MULT       = 5, &
+       MG_GET_XS_SCATTER_FMU_MULT   = 6, &
+       MG_GET_XS_SCATTER_FMU        = 7, &
+       MG_GET_XS_FISSION            = 8, &
+       MG_GET_XS_KAPPA_FISSION      = 9, &
+       MG_GET_XS_PROMPT_NU_FISSION  = 10, &
+       MG_GET_XS_DELAYED_NU_FISSION = 11, &
+       MG_GET_XS_NU_FISSION         = 12, &
+       MG_GET_XS_CHI_PROMPT         = 13, &
+       MG_GET_XS_CHI_DELAYED        = 14
+
   ! ============================================================================
   ! RANDOM NUMBER STREAM CONSTANTS
 

src/constants.h

@@ -1,10 +1,85 @@
+//! \file constants.h
+//! A collection of constants
+
 #ifndef CONSTANTS_H
 #define CONSTANTS_H
 
+#include <cmath>
+#include <array>
 #include <limits>
+#include <vector>
 
+namespace openmc {
 
-namespace openmc{
+typedef std::vector<double> double_1dvec;
+typedef std::vector<std::vector<double> > double_2dvec;
+typedef std::vector<std::vector<std::vector<double> > > double_3dvec;
+typedef std::vector<std::vector<std::vector<std::vector<double> > > > double_4dvec;
+typedef std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > > double_5dvec;
+typedef std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > > > double_6dvec;
+typedef std::vector<int> int_1dvec;
+typedef std::vector<std::vector<int> > int_2dvec;
+typedef std::vector<std::vector<std::vector<int> > > int_3dvec;
+
+constexpr int MAX_SAMPLE {10000};
+
+constexpr std::array<int, 3> VERSION {0, 10, 0};
+constexpr std::array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
+
+// Maximum number of words in a single line, length of line, and length of
+// single word
+constexpr int MAX_WORDS {500};
+constexpr int MAX_LINE_LEN {250};
+constexpr int MAX_WORD_LEN {150};
+constexpr int MAX_FILE_LEN {255};
+
+// Physical Constants
+constexpr double K_BOLTZMANN {8.6173303e-5};  // Boltzmann constant in eV/K
+
+// Angular distribution type
+constexpr int ANGLE_ISOTROPIC {1};
+constexpr int ANGLE_32_EQUI {2};
+constexpr int ANGLE_TABULAR {3};
+constexpr int ANGLE_LEGENDRE {4};
+constexpr int ANGLE_HISTOGRAM {5};
+
+// MGXS Table Types
+constexpr int MGXS_ISOTROPIC {1}; // Isotroically weighted data
+constexpr int MGXS_ANGLE {2};     // Data by angular bins
+
+// Flag to denote this was a macroscopic data object
+constexpr double MACROSCOPIC_AWR {-2.};
+
+// Number of mu bins to use when converting Legendres to tabular type
+constexpr int DEFAULT_NMU {33};
+
+// Temperature treatment method
+constexpr int TEMPERATURE_NEAREST {1};
+constexpr int TEMPERATURE_INTERPOLATION {2};
+
+// TODO: cmath::M_PI has 3 more digits precision than the Fortran constant we
+// use so for now we will reuse the Fortran constant until we are OK with
+// modifying test results
+constexpr double PI {3.1415926535898};
+
+const double SQRT_PI {std::sqrt(PI)};
+
+// Mgxs::get_xs enumerated types
+constexpr int MG_GET_XS_TOTAL              {0};
+constexpr int MG_GET_XS_ABSORPTION         {1};
+constexpr int MG_GET_XS_INVERSE_VELOCITY   {2};
+constexpr int MG_GET_XS_DECAY_RATE         {3};
+constexpr int MG_GET_XS_SCATTER            {4};
+constexpr int MG_GET_XS_SCATTER_MULT       {5};
+constexpr int MG_GET_XS_SCATTER_FMU_MULT   {6};
+constexpr int MG_GET_XS_SCATTER_FMU        {7};
+constexpr int MG_GET_XS_FISSION            {8};
+constexpr int MG_GET_XS_KAPPA_FISSION      {9};
+constexpr int MG_GET_XS_PROMPT_NU_FISSION  {10};
+constexpr int MG_GET_XS_DELAYED_NU_FISSION {11};
+constexpr int MG_GET_XS_NU_FISSION         {12};
+constexpr int MG_GET_XS_CHI_PROMPT         {13};
+constexpr int MG_GET_XS_CHI_DELAYED        {14};
 
 extern "C" double FP_COINCIDENT;
 extern "C" double FP_PRECISION;

src/error.F90

@@ -265,4 +265,13 @@ contains
 
   end subroutine write_message
 
+  subroutine write_message_from_c(message, message_len, level) bind(C)
+    integer(C_INT),         intent(in), value :: message_len
+    character(kind=C_CHAR), intent(in)        :: message(message_len)
+    integer(C_INT),         intent(in), value :: level
+    character(message_len+1)                  :: message_out
+    write(message_out, *) message
+    call write_message(message_out, level)
+  end subroutine write_message_from_c
+
 end module error

src/error.h

@@ -11,7 +11,8 @@ namespace openmc {
 
 extern "C" void fatal_error_from_c(const char* message, int message_len);
 extern "C" void warning_from_c(const char* message, int message_len);
-
+extern "C" void write_message_from_c(const char* message, int message_len,
+                                     int level);
 
 inline
 void fatal_error(const char *message)
@@ -43,5 +44,23 @@ void warning(const std::stringstream& message)
   warning(message.str());
 }
 
+inline
+void write_message(const char* message, int level)
+{
+  write_message_from_c(message, strlen(message), level);
+}
+
+inline
+void write_message(const std::string& message, int level)
+{
+  write_message_from_c(message.c_str(), message.length(), level);
+}
+
+inline
+void write_message(const std::stringstream& message, int level)
+{
+  write_message(message.str(), level);
+}
+
 } // namespace openmc
 #endif // ERROR_H

src/hdf5_interface.cpp

@@ -447,6 +447,172 @@ read_complex(hid_t obj_id, const char* name, double _Complex* buffer, bool indep
 }
 
 
+void
+read_nd_vector(hid_t obj_id, const char* name, std::vector<double>& result,
+               bool must_have)
+{
+  if (object_exists(obj_id, name)) {
+    read_double(obj_id, name, result.data(), true);
+  } else if (must_have) {
+    fatal_error(std::string("Must provide " + std::string(name) + "!"));
+  }
+}
+
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<double> >& result, bool must_have)
+{
+  if (object_exists(obj_id, name)) {
+    int dim1 = result.size();
+    int dim2 = result[0].size();
+    double temp_arr[dim1 * dim2];
+    read_double(obj_id, name, temp_arr, true);
+
+    int temp_idx = 0;
+    for (int i = 0; i < dim1; i++) {
+      for (int j = 0; j < dim2; j++) {
+        result[i][j] = temp_arr[temp_idx++];
+      }
+    }
+  } else if (must_have) {
+    fatal_error(std::string("Must provide " + std::string(name) + "!"));
+  }
+}
+
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<int> >& result, bool must_have)
+{
+  if (object_exists(obj_id, name)) {
+    int dim1 = result.size();
+    int dim2 = result[0].size();
+    int temp_arr[dim1 * dim2];
+    read_int(obj_id, name, temp_arr, true);
+
+    int temp_idx = 0;
+    for (int i = 0; i < dim1; i++) {
+      for (int j = 0; j < dim2; j++) {
+        result[i][j] = temp_arr[temp_idx++];
+      }
+    }
+  } else if (must_have) {
+    fatal_error(std::string("Must provide " + std::string(name) + "!"));
+  }
+}
+
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<std::vector<double> > >& result,
+               bool must_have)
+{
+  if (object_exists(obj_id, name)) {
+    int dim1 = result.size();
+    int dim2 = result[0].size();
+    int dim3 = result[0][0].size();
+    double temp_arr[dim1 * dim2 * dim3];
+    read_double(obj_id, name, temp_arr, true);
+
+    int temp_idx = 0;
+    for (int i = 0; i < dim1; i++) {
+      for (int j = 0; j < dim2; j++) {
+        for (int k = 0; k < dim3; k++) {
+          result[i][j][k] = temp_arr[temp_idx++];
+        }
+      }
+    }
+  } else if (must_have) {
+    fatal_error(std::string("Must provide " + std::string(name) + "!"));
+  }
+}
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<std::vector<int> > >& result,
+               bool must_have)
+{
+  if (object_exists(obj_id, name)) {
+    int dim1 = result.size();
+    int dim2 = result[0].size();
+    int dim3 = result[0][0].size();
+    int temp_arr[dim1 * dim2 * dim3];
+    read_int(obj_id, name, temp_arr, true);
+
+    int temp_idx = 0;
+    for (int i = 0; i < dim1; i++) {
+      for (int j = 0; j < dim2; j++) {
+        for (int k = 0; k < dim3; k++) {
+          result[i][j][k] = temp_arr[temp_idx++];
+        }
+      }
+    }
+  } else if (must_have) {
+    fatal_error(std::string("Must provide " + std::string(name) + "!"));
+  }
+}
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<std::vector<std::vector<double> > > >& result,
+               bool must_have)
+{
+  if (object_exists(obj_id, name)) {
+    int dim1 = result.size();
+    int dim2 = result[0].size();
+    int dim3 = result[0][0].size();
+    int dim4 = result[0][0][0].size();
+    double temp_arr[dim1 * dim2 * dim3 * dim4];
+    read_double(obj_id, name, temp_arr, true);
+
+    int temp_idx = 0;
+    for (int i = 0; i < dim1; i++) {
+      for (int j = 0; j < dim2; j++) {
+        for (int k = 0; k < dim3; k++) {
+          for (int l = 0; l < dim4; l++) {
+            result[i][j][k][l] = temp_arr[temp_idx++];
+          }
+        }
+      }
+    }
+  } else if (must_have) {
+    fatal_error(std::string("Must provide " + std::string(name) + "!"));
+  }
+}
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > >& result,
+               bool must_have)
+{
+  if (object_exists(obj_id, name)) {
+    int dim1 = result.size();
+    int dim2 = result[0].size();
+    int dim3 = result[0][0].size();
+    int dim4 = result[0][0][0].size();
+    int dim5 = result[0][0][0][0].size();
+    double temp_arr[dim1 * dim2 * dim3 * dim4 * dim5];
+    read_double(obj_id, name, temp_arr, true);
+
+    int temp_idx = 0;
+    for (int i = 0; i < dim1; i++) {
+      for (int j = 0; j < dim2; j++) {
+        for (int k = 0; k < dim3; k++) {
+          for (int l = 0; l < dim4; l++) {
+            for (int m = 0; m < dim5; m++) {
+              result[i][j][k][l][m] = temp_arr[temp_idx++];
+            }
+          }
+        }
+      }
+    }
+  } else if (must_have) {
+    fatal_error(std::string("Must provide " + std::string(name) + "!"));
+  }
+}
+
+
 void
 read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results)
 {

src/hdf5_interface.h

@@ -7,6 +7,7 @@
 #include <array>
 #include <string>
 #include <sstream>
+#include <vector>
 #include <complex.h>
 
 
@@ -55,6 +56,39 @@ extern "C" void read_string(hid_t obj_id, const char* name, size_t slen,
 extern "C" void read_complex(hid_t obj_id, const char* name,
                              double _Complex* buffer, bool indep);
 
+void
+read_nd_vector(hid_t obj_id, const char* name, std::vector<double>& result,
+               bool must_have = false);
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<double> >& result,
+               bool must_have = false);
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<int> >& result, bool must_have = false);
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<std::vector<double> > >& result,
+               bool must_have = false);
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<std::vector<int> > >& result,
+               bool must_have = false);
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<std::vector<std::vector<double> > > >& result,
+               bool must_have = false);
+
+void
+read_nd_vector(hid_t obj_id, const char* name,
+               std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > >& result,
+               bool must_have = false);
+
 extern "C" void read_tally_results(hid_t group_id, hsize_t n_filter,
                                    hsize_t n_score, double* results);
 

src/input_xml.F90

@@ -19,7 +19,7 @@ module input_xml
   use mesh_header
   use message_passing
   use mgxs_data,        only: create_macro_xs, read_mgxs
-  use mgxs_header
+  use mgxs_interface
   use nuclide_header
   use output,           only: title, header, print_plot
   use plot_header
@@ -3458,7 +3458,7 @@ contains
           if (run_CE) then
             awr = nuclides(mat % nuclide(j)) % awr
           else
-            awr = nuclides_MG(mat % nuclide(j)) % obj % awr
+            awr = get_awr_c(mat % nuclide(j))
           end if
 
           ! if given weight percent, convert all values so that they are divided
@@ -3483,7 +3483,7 @@ contains
             if (run_CE) then
               awr = nuclides(mat % nuclide(j)) % awr
             else
-              awr = nuclides_MG(mat % nuclide(j)) % obj % awr
+              awr = get_awr_c(mat % nuclide(j))
             end if
             x = mat % atom_density(j)
             sum_percent = sum_percent + x*awr

src/material_header.F90

@@ -34,7 +34,7 @@ module material_header
     integer              :: n_nuclides = 0  ! number of nuclides
     integer, allocatable :: nuclide(:)      ! index in nuclides array
     real(8)              :: density         ! total atom density in atom/b-cm
-    real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
+    real(C_DOUBLE), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
     real(8)              :: density_gpcc    ! total density in g/cm^3
 
     ! To improve performance of tallying, we store an array (direct address

src/math_functions.cpp

@@ -97,7 +97,7 @@ void calc_pn_c(int n, double x, double pnx[]) {
   }
 
   // Use recursion relation to build the higher orders
-  for (int l = 1; l < n; l ++) {
+  for (int l = 1; l < n; l++) {
     pnx[l + 1] = ((2 * l + 1) * x * pnx[l] - l * pnx[l - 1]) / (l + 1);
   }
 }

src/math_functions.h

@@ -7,22 +7,12 @@
 #include <cmath>
 #include <cstdlib>
 
+#include "constants.h"
 #include "random_lcg.h"
 
 
 namespace openmc {
 
-//==============================================================================
-// Module constants.
-//==============================================================================
-
-// TODO: cmath::M_PI has 3 more digits precision than the Fortran constant we
-// use so for now we will reuse the Fortran constant until we are OK with
-// modifying test results
-extern "C" constexpr double PI {3.1415926535898};
-
-extern "C" const double SQRT_PI {std::sqrt(PI)};
-
 //==============================================================================
 //! Calculate the percentile of the standard normal distribution with a
 //! specified probability level.

src/mgxs.cpp

@@ -0,0 +1,712 @@
+#include <cmath>
+#include <cstdlib>
+#include <algorithm>
+#include <valarray>
+
+ #ifdef _OPENMP
+ # include <omp.h>
+ #endif
+
+#include "error.h"
+#include "math_functions.h"
+#include "random_lcg.h"
+#include "string_functions.h"
+#include "mgxs.h"
+
+
+namespace openmc {
+
+// Storage for the MGXS data
+std::vector<Mgxs> nuclides_MG;
+std::vector<Mgxs> macro_xs;
+
+
+//==============================================================================
+// Mgxs base-class methods
+//==============================================================================
+
+void
+Mgxs::init(const std::string& in_name, double in_awr,
+     const double_1dvec& in_kTs, bool in_fissionable, int in_scatter_format,
+     int in_num_groups, int in_num_delayed_groups, bool in_is_isotropic,
+     const double_1dvec& in_polar, const double_1dvec& in_azimuthal)
+{
+  // Set the metadata
+  name = in_name;
+  awr = in_awr;
+  kTs = in_kTs;
+  fissionable = in_fissionable;
+  scatter_format = in_scatter_format;
+  num_groups = in_num_groups;
+  num_delayed_groups = in_num_delayed_groups;
+  xs.resize(in_kTs.size());
+  is_isotropic = in_is_isotropic;
+  n_pol = in_polar.size();
+  n_azi = in_azimuthal.size();
+  polar = in_polar;
+  azimuthal = in_azimuthal;
+
+  // Set the cross section index cache
+#ifdef _OPENMP
+  int n_threads = omp_get_max_threads();
+#else
+  int n_threads = 1;
+#endif
+  cache.resize(n_threads);
+  // std::vector.resize() will value-initialize the members of cache[:]
+}
+
+//==============================================================================
+
+void
+Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
+     int in_num_delayed_groups, const double_1dvec& temperature,
+     double tolerance, int_1dvec& temps_to_read, int& order_dim, int& method)
+{
+  // get name
+  char char_name[MAX_WORD_LEN];
+  get_name(xs_id, char_name);
+  std::string in_name {char_name};
+  // remove the leading '/'
+  in_name = in_name.substr(1);
+
+  // Get the AWR
+  double in_awr;
+  if (attribute_exists(xs_id, "atomic_weight_ratio")) {
+    read_attr_double(xs_id, "atomic_weight_ratio", &in_awr);
+  } else {
+    in_awr = MACROSCOPIC_AWR;
+  }
+
+  // Determine the available temperatures
+  hid_t kT_group = open_group(xs_id, "kTs");
+  int num_temps = get_num_datasets(kT_group);
+  char* dset_names[num_temps];
+  for (int i = 0; i < num_temps; i++) {
+    dset_names[i] = new char[151];
+  }
+  get_datasets(kT_group, dset_names);
+  double_1dvec available_temps(num_temps);
+  for (int i = 0; i < num_temps; i++) {
+    read_double(kT_group, dset_names[i], &available_temps[i], true);
+
+    // convert eV to Kelvin
+    available_temps[i] /= K_BOLTZMANN;
+
+    // Done with dset_names, so delete it
+    delete[] dset_names[i];
+  }
+  std::sort(available_temps.begin(), available_temps.end());
+
+  // If only one temperature is available, lets just use nearest temperature
+  // interpolation
+  if ((num_temps == 1) && (method == TEMPERATURE_INTERPOLATION)) {
+    warning("Cross sections for " + strtrim(name) + " are only available " +
+            "at one temperature.  Reverting to the nearest temperature " +
+            "method.");
+    method = TEMPERATURE_NEAREST;
+  }
+
+  switch(method) {
+    case TEMPERATURE_NEAREST:
+      // Find the minimum difference
+      for (int i = 0; i < temperature.size(); i++) {
+        std::valarray<double> temp_diff(available_temps.data(),
+                                        available_temps.size());
+        temp_diff = std::abs(temp_diff - temperature[i]);
+        int i_closest = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
+             std::begin(temp_diff);
+        double temp_actual = available_temps[i_closest];
+
+        if (std::abs(temp_actual - temperature[i]) < tolerance) {
+          if (std::find(temps_to_read.begin(), temps_to_read.end(),
+                        std::round(temp_actual)) == temps_to_read.end()) {
+            temps_to_read.push_back(std::round(temp_actual));
+          } else {
+            fatal_error("MGXS Library does not contain cross section for " +
+                        in_name + " at or near " +
+                        std::to_string(std::round(temperature[i])) + " K.");
+          }
+        }
+      }
+      break;
+
+    case TEMPERATURE_INTERPOLATION:
+      for (int i = 0; i < temperature.size(); i++) {
+        for (int j = 0; j < num_temps - 1; j++) {
+          if ((available_temps[j] <= temperature[i]) &&
+              (temperature[i] < available_temps[j + 1])) {
+            if (std::find(temps_to_read.begin(),
+                          temps_to_read.end(),
+                          std::round(available_temps[j])) == temps_to_read.end()) {
+              temps_to_read.push_back(std::round((int)available_temps[j]));
+            }
+
+            if (std::find(temps_to_read.begin(), temps_to_read.end(),
+                          std::round(available_temps[j + 1])) == temps_to_read.end()) {
+              temps_to_read.push_back(std::round((int) available_temps[j + 1]));
+            }
+            continue;
+          }
+        }
+
+      fatal_error("MGXS Library does not contain cross sections for " +
+                  in_name + " at temperatures that bound " +
+                  std::to_string(std::round(temperature[i])));
+      }
+  }
+  std::sort(temps_to_read.begin(), temps_to_read.end());
+
+  // Get the library's temperatures
+  int n_temperature = temps_to_read.size();
+  double_1dvec in_kTs(n_temperature);
+  for (int i = 0; i < n_temperature; i++) {
+    std::string temp_str(std::to_string(temps_to_read[i]) + "K");
+
+    //read exact temperature value
+    read_double(kT_group, temp_str.c_str(), &in_kTs[i], true);
+  }
+  close_group(kT_group);
+
+  // Load the remaining metadata
+  int in_scatter_format;
+  if (attribute_exists(xs_id, "scatter_format")) {
+    std::string temp_str(MAX_WORD_LEN, ' ');
+    read_attr_string(xs_id, "scatter_format", MAX_WORD_LEN, &temp_str[0]);
+    to_lower(strtrim(temp_str));
+    if (temp_str.compare(0, 8, "legendre") == 0) {
+      in_scatter_format = ANGLE_LEGENDRE;
+    } else if (temp_str.compare(0, 9, "histogram") == 0) {
+      in_scatter_format = ANGLE_HISTOGRAM;
+    } else if (temp_str.compare(0, 7, "tabular") == 0) {
+      in_scatter_format = ANGLE_TABULAR;
+    } else {
+      fatal_error("Invalid scatter_format option!");
+    }
+  } else {
+    in_scatter_format = ANGLE_LEGENDRE;
+  }
+
+  if (attribute_exists(xs_id, "scatter_shape")) {
+    std::string temp_str(MAX_WORD_LEN, ' ');
+    read_attr_string(xs_id, "scatter_shape", MAX_WORD_LEN, &temp_str[0]);
+    to_lower(strtrim(temp_str));
+    if (temp_str.compare(0, 14, "[g][g\'][order]") != 0) {
+      fatal_error("Invalid scatter_shape option!");
+    }
+  }
+
+  bool in_fissionable = false;
+  if (attribute_exists(xs_id, "fissionable")) {
+    int int_fiss;
+    read_attr_int(xs_id, "fissionable", &int_fiss);
+    in_fissionable = int_fiss;
+  } else {
+    fatal_error("Fissionable element must be set!");
+  }
+
+  // Get the library's value for the order
+  if (attribute_exists(xs_id, "order")) {
+    read_attr_int(xs_id, "order", &order_dim);
+  } else {
+    fatal_error("Order must be provided!");
+  }
+
+  // Store the dimensionality of the data in order_dim.
+  // For Legendre data, we usually refer to it as Pn where n is the order.
+  // However Pn has n+1 sets of points (since you need to count the P0
+  // moment). Adjust for that. Histogram and Tabular formats dont need this
+  // adjustment.
+  if (in_scatter_format == ANGLE_LEGENDRE) {
+    ++order_dim;
+  }
+
+  // Get the angular information
+  int in_n_pol;
+  int in_n_azi;
+  bool in_is_isotropic = true;
+  if (attribute_exists(xs_id, "representation")) {
+    std::string temp_str(MAX_WORD_LEN, ' ');
+    read_attr_string(xs_id, "representation", MAX_WORD_LEN, &temp_str[0]);
+    to_lower(strtrim(temp_str));
+    if (temp_str.compare(0, 5, "angle") == 0) {
+      in_is_isotropic = false;
+    } else if (temp_str.compare(0, 9, "isotropic") != 0) {
+      fatal_error("Invalid Data Representation!");
+    }
+  }
+
+  if (!in_is_isotropic) {
+    if (attribute_exists(xs_id, "num_polar")) {
+      read_attr_int(xs_id, "num_polar", &in_n_pol);
+    } else {
+      fatal_error("num_polar must be provided!");
+    }
+    if (attribute_exists(xs_id, "num_azimuthal")) {
+      read_attr_int(xs_id, "num_azimuthal", &in_n_azi);
+    } else {
+      fatal_error("num_azimuthal must be provided!");
+    }
+  } else {
+    in_n_pol = 1;
+    in_n_azi = 1;
+  }
+
+  // Set the angular bins to use equally-spaced bins
+  double_1dvec in_polar(in_n_pol);
+  double dangle = PI / in_n_pol;
+  for (int p = 0; p < in_n_pol; p++) {
+    in_polar[p]  = (p + 0.5) * dangle;
+  }
+  double_1dvec in_azimuthal(in_n_azi);
+  dangle = 2. * PI / in_n_azi;
+  for (int a = 0; a < in_n_azi; a++) {
+    in_azimuthal[a] = (a + 0.5) * dangle - PI;
+  }
+
+  // Finally use this data to initialize the MGXS Object
+  init(in_name, in_awr, in_kTs, in_fissionable, in_scatter_format,
+       in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar,
+       in_azimuthal);
+}
+
+//==============================================================================
+
+Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups,
+     const double_1dvec& temperature, double tolerance, int max_order,
+     bool legendre_to_tabular, int legendre_to_tabular_points, int& method)
+{
+  // Call generic data gathering routine (will populate the metadata)
+  int order_data;
+  int_1dvec temps_to_read;
+  metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature,
+       tolerance, temps_to_read, order_data, method);
+
+  // Set number of energy and delayed groups
+  int final_scatter_format = scatter_format;
+  if (legendre_to_tabular) {
+    if (scatter_format == ANGLE_LEGENDRE) final_scatter_format = ANGLE_TABULAR;
+  }
+
+  // Load the more specific XsData information
+  for (int t = 0; t < temps_to_read.size(); t++) {
+    xs[t] = XsData(energy_groups, delayed_groups, fissionable,
+                   final_scatter_format, n_pol, n_azi);
+    // Get the temperature as a string and then open the HDF5 group
+    std::string temp_str = std::to_string(temps_to_read[t]) + "K";
+    hid_t xsdata_grp = open_group(xs_id, temp_str.c_str());
+
+    xs[t].from_hdf5(xsdata_grp, fissionable, scatter_format,
+                    final_scatter_format, order_data, max_order,
+                    legendre_to_tabular_points, is_isotropic, n_pol, n_azi);
+    close_group(xsdata_grp);
+
+  } // end temperature loop
+
+  // Make sure the scattering format is updated to the final case
+  scatter_format = final_scatter_format;
+}
+
+//==============================================================================
+
+Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
+     const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
+     double tolerance, int& method)
+{
+  // Get the minimum data needed to initialize:
+  // Dont need awr, but lets just initialize it anyways
+  double in_awr = -1.;
+  // start with the assumption it is not fissionable
+  bool in_fissionable = false;
+  for (int m = 0; m < micros.size(); m++) {
+    if (micros[m]->fissionable) in_fissionable = true;
+  }
+  // Force all of the following data to be the same; these will be verified
+  // to be true later
+  int in_scatter_format = micros[0]->scatter_format;
+  int in_num_groups = micros[0]->num_groups;
+  int in_num_delayed_groups = micros[0]->num_delayed_groups;
+  bool in_is_isotropic = micros[0]->is_isotropic;
+  double_1dvec in_polar = micros[0]->polar;
+  double_1dvec in_azimuthal = micros[0]->azimuthal;
+
+  init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format,
+       in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar,
+       in_azimuthal);
+
+  // Create the xs data for each temperature
+  for (int t = 0; t < mat_kTs.size(); t++) {
+    xs[t] = XsData(in_num_groups, in_num_delayed_groups, in_fissionable,
+       in_scatter_format, in_polar.size(), in_azimuthal.size());
+
+    // Find the right temperature index to use
+    double temp_desired = mat_kTs[t];
+
+    // Create the list of temperature indices and interpolation factors for
+    // each microscopic data at the material temperature
+    int_1dvec micro_t(micros.size(), 0);
+    double_1dvec micro_t_interp(micros.size(), 0.);
+    for (int m = 0; m < micros.size(); m++) {
+      switch(method) {
+      case TEMPERATURE_NEAREST:
+        {
+          // Find the nearest temperature
+          std::valarray<double> temp_diff(micros[m]->kTs.data(),
+                                          micros[m]->kTs.size());
+          temp_diff = std::abs(temp_diff - temp_desired);
+          micro_t[m] = std::min_element(std::begin(temp_diff),
+                                        std::end(temp_diff)) -
+               std::begin(temp_diff);
+          double temp_actual = micros[m]->kTs[micro_t[m]];
+
+          if (std::abs(temp_actual - temp_desired) >= K_BOLTZMANN * tolerance) {
+            fatal_error("MGXS Library does not contain cross section for " +
+                        name + " at or near " +
+                        std::to_string(std::round(temp_desired / K_BOLTZMANN))
+                        + " K.");
+          }
+        }
+        break;
+      case TEMPERATURE_INTERPOLATION:
+        // Get a list of bounding temperatures for each actual temperature
+        // present in the model
+        for (int k = 0; k < micros[m]->kTs.size() - 1; k++) {
+          if ((micros[m]->kTs[k] <= temp_desired) &&
+              (temp_desired < micros[m]->kTs[k + 1])) {
+            micro_t[m] = k;
+            if (k == 0) {
+              micro_t_interp[m] = (temp_desired - micros[m]->kTs[k]) /
+                   (micros[m]->kTs[k + 1] - micros[m]->kTs[k]);
+            } else {
+              micro_t_interp[m] = 1.;
+            }
+          }
+        }
+      } // end switch
+    } // end microscopic temperature loop
+
+    // We are about to loop through each of the microscopic objects
+    // and incorporate the contribution of each microscopic data at
+    // one of the two temperature interpolants to this macroscopic quantity.
+    // If we are doing nearest temperature interpolation, then we don't need
+    // to do the 2nd temperature
+    int num_interp_points = 2;
+    if (method == TEMPERATURE_NEAREST) num_interp_points = 1;
+    for (int interp_point = 0; interp_point < num_interp_points; interp_point++) {
+      double_1dvec interp(micros.size());
+      double_1dvec temp_indices(micros.size());
+      for (int m = 0; m < micros.size(); m++) {
+        interp[m] = (1. - micro_t_interp[m]) * atom_densities[m];
+        temp_indices[m] = micro_t[m] + interp_point;
+      }
+
+      combine(micros, interp, micro_t, t);
+    } // end loop to sum all micros across the temperatures
+  } // end temperature (t) loop
+}
+
+//==============================================================================
+
+void
+Mgxs::combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
+              const int_1dvec& micro_ts, int this_t)
+{
+  // Build the vector of pointers to the xs objects within micros
+  std::vector<XsData*> those_xs(micros.size());
+  for (int i = 0; i < micros.size(); i++) {
+    if (!xs[this_t].equiv(micros[i]->xs[micro_ts[i]])) {
+      fatal_error("Cannot combine the Mgxs objects!");
+    }
+    those_xs[i] = &(micros[i]->xs[micro_ts[i]]);
+  }
+
+  xs[this_t].combine(those_xs, scalars);
+}
+
+//==============================================================================
+
+double
+Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
+{
+  // This method assumes that the temperature and angle indices are set
+#ifdef _OPENMP
+  int tid = omp_get_thread_num();
+  XsData* xs_t = &xs[cache[tid].t];
+  int a = cache[tid].a;
+#else
+  XsData* xs_t = &xs[cache[0].t];
+  int a = cache[0].a;
+#endif
+  double val;
+  switch(xstype) {
+  case MG_GET_XS_TOTAL:
+    val = xs_t->total[a][gin];
+    break;
+  case MG_GET_XS_NU_FISSION:
+    val = fissionable ? xs_t->nu_fission[a][gin] : 0.;
+    break;
+  case MG_GET_XS_ABSORPTION:
+    val = xs_t->absorption[a][gin];
+    break;
+  case MG_GET_XS_FISSION:
+    val = fissionable ? xs_t->fission[a][gin] : 0.;
+    break;
+  case MG_GET_XS_KAPPA_FISSION:
+    val = fissionable ? xs_t->kappa_fission[a][gin] : 0.;
+    break;
+  case MG_GET_XS_SCATTER:
+  case MG_GET_XS_SCATTER_MULT:
+  case MG_GET_XS_SCATTER_FMU_MULT:
+  case MG_GET_XS_SCATTER_FMU:
+    val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
+    break;
+  case MG_GET_XS_PROMPT_NU_FISSION:
+    val = fissionable ? xs_t->prompt_nu_fission[a][gin] : 0.;
+    break;
+  case MG_GET_XS_DELAYED_NU_FISSION:
+    if (fissionable) {
+      if (dg != nullptr) {
+        val = xs_t->delayed_nu_fission[a][gin][*dg];
+      } else {
+        val = 0.;
+        for (auto& num : xs_t->delayed_nu_fission[a][gin]) {
+          val += num;
+        }
+      }
+    } else {
+      val = 0.;
+    }
+    break;
+  case MG_GET_XS_CHI_PROMPT:
+    if (fissionable) {
+      if (gout != nullptr) {
+        val = xs_t->chi_prompt[a][gin][*gout];
+      } else {
+        // provide an outgoing group-wise sum
+        val = 0.;
+        for (auto& num : xs_t->chi_prompt[a][gin]) {
+          val += num;
+        }
+      }
+    } else {
+      val = 0.;
+    }
+    break;
+  case MG_GET_XS_CHI_DELAYED:
+    if (fissionable) {
+      if (gout != nullptr) {
+        if (dg != nullptr) {
+          val = xs_t->chi_delayed[a][gin][*gout][*dg];
+        } else {
+          val = xs_t->chi_delayed[a][gin][*gout][0];
+        }
+      } else {
+        if (dg != nullptr) {
+          val = 0.;
+          for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
+            val += xs_t->chi_delayed[a][gin][i][*dg];
+          }
+        } else {
+          val = 0.;
+          for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
+            for (auto& num : xs_t->chi_delayed[a][gin][i]) {
+              val += num;
+            }
+          }
+        }
+      }
+    } else {
+      val = 0.;
+    }
+    break;
+  case MG_GET_XS_INVERSE_VELOCITY:
+    val = xs_t->inverse_velocity[a][gin];
+    break;
+  case MG_GET_XS_DECAY_RATE:
+    if (dg != nullptr) {
+      val = xs_t->decay_rate[a][*dg + 1];
+    } else {
+      val = xs_t->decay_rate[a][0];
+    }
+    break;
+  default:
+    val = 0.;
+  }
+  return val;
+}
+
+//==============================================================================
+
+void
+Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
+{
+  // This method assumes that the temperature and angle indices are set
+#ifdef _OPENMP
+  int tid = omp_get_thread_num();
+#else
+  int tid = 0;
+#endif
+  XsData* xs_t = &xs[cache[tid].t];
+  double nu_fission = xs_t->nu_fission[cache[tid].a][gin];
+
+  // Find the probability of having a prompt neutron
+  double prob_prompt =
+       xs_t->prompt_nu_fission[cache[tid].a][gin];
+
+  // sample random numbers
+  double xi_pd = prn() * nu_fission;
+  double xi_gout = prn();
+
+  // Select whether the neutron is prompt or delayed
+  if (xi_pd <= prob_prompt) {
+    // the neutron is prompt
+
+    // set the delayed group for the particle to be -1, indicating prompt
+    dg = -1;
+
+    // sample the outgoing energy group
+    gout = 0;
+    double prob_gout =
+         xs_t->chi_prompt[cache[tid].a][gin][gout];
+    while (prob_gout < xi_gout) {
+      gout++;
+      prob_gout += xs_t->chi_prompt[cache[tid].a][gin][gout];
+    }
+
+  } else {
+    // the neutron is delayed
+
+    // get the delayed group
+    dg = 0;
+    while (xi_pd >= prob_prompt) {
+      dg++;
+      prob_prompt +=
+           xs_t->delayed_nu_fission[cache[tid].a][gin][dg];
+    }
+
+    // adjust dg in case of round-off error
+    dg = std::min(dg, num_delayed_groups - 1);
+
+    // sample the outgoing energy group
+    gout = 0;
+    double prob_gout =
+         xs_t->chi_delayed[cache[tid].a][gin][gout][dg];
+    while (prob_gout < xi_gout) {
+      gout++;
+      prob_gout +=
+           xs_t->chi_delayed[cache[tid].a][gin][gout][dg];
+    }
+  }
+}
+
+//==============================================================================
+
+void
+Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt)
+{
+  // This method assumes that the temperature and angle indices are set
+  // Sample the data
+#ifdef _OPENMP
+  int tid = omp_get_thread_num();
+#else
+  int tid = 0;
+#endif
+  xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt);
+}
+
+//==============================================================================
+
+void
+Mgxs::calculate_xs(int gin, double sqrtkT, const double uvw[3],
+     double& total_xs, double& abs_xs, double& nu_fiss_xs)
+{
+  // Set our indices
+#ifdef _OPENMP
+  int tid = omp_get_thread_num();
+#else
+  int tid = 0;
+#endif
+  set_temperature_index(sqrtkT);
+  set_angle_index(uvw);
+  XsData* xs_t = &xs[cache[tid].t];
+  total_xs = xs_t->total[cache[tid].a][gin];
+  abs_xs = xs_t->absorption[cache[tid].a][gin];
+
+  nu_fiss_xs = fissionable ? xs_t->nu_fission[cache[tid].a][gin] : 0.;
+}
+
+//==============================================================================
+
+bool
+Mgxs::equiv(const Mgxs& that)
+{
+  return ((num_delayed_groups == that.num_delayed_groups) &&
+       (num_groups == that.num_groups) &&
+       (n_pol == that.n_pol) &&
+       (n_azi == that.n_azi) &&
+       (std::equal(polar.begin(), polar.end(), that.polar.begin())) &&
+       (std::equal(azimuthal.begin(), azimuthal.end(), that.azimuthal.begin())) &&
+       (scatter_format == that.scatter_format));
+}
+
+//==============================================================================
+
+void
+Mgxs::set_temperature_index(double sqrtkT)
+{
+  // See if we need to find the new index
+#ifdef _OPENMP
+  int tid = omp_get_thread_num();
+#else
+  int tid = 0;
+#endif
+  if (sqrtkT != cache[tid].sqrtkT) {
+    double kT = sqrtkT * sqrtkT;
+
+    // initialize vector for storage of the differences
+    std::valarray<double> temp_diff(kTs.data(), kTs.size());
+
+    // Find the minimum difference of kT and kTs
+    temp_diff = std::abs(temp_diff - kT);
+    cache[tid].t = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
+         std::begin(temp_diff);
+
+    // store this temperature as the last one used
+    cache[tid].sqrtkT = sqrtkT;
+  }
+}
+
+//==============================================================================
+
+void
+Mgxs::set_angle_index(const double uvw[3])
+{
+  // See if we need to find the new index
+#ifdef _OPENMP
+  int tid = omp_get_thread_num();
+#else
+  int tid = 0;
+#endif
+  if (!is_isotropic &&
+      ((uvw[0] != cache[tid].u) || (uvw[1] != cache[tid].v) ||
+       (uvw[2] != cache[tid].w))) {
+    // convert uvw to polar and azimuthal angles
+    double my_pol = std::acos(uvw[2]);
+    double my_azi = std::atan2(uvw[1], uvw[0]);
+
+    // Find the location, assuming equal-bin angles
+    double delta_angle = PI / n_pol;
+    int p = std::floor(my_pol / delta_angle);
+    delta_angle = 2. * PI / n_azi;
+    int a = std::floor((my_azi + PI) / delta_angle);
+
+    cache[tid].a = n_azi * p + a;
+
+    // store this direction as the last one used
+    cache[tid].u = uvw[0];
+    cache[tid].v = uvw[1];
+    cache[tid].w = uvw[2];
+  }
+}
+
+} // namespace openmc

src/mgxs.h

@@ -0,0 +1,205 @@
+//! \file mgxs.h
+//! A collection of classes for Multi-Group Cross Section data
+
+#ifndef MGXS_H
+#define MGXS_H
+
+#include <string>
+#include <vector>
+
+#include "constants.h"
+#include "hdf5_interface.h"
+#include "xsdata.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// Cache contains the cached data for an MGXS object
+//==============================================================================
+
+struct CacheData {
+  double sqrtkT; // last temperature corresponding to t
+  int t; // temperature index
+  int a; // angle index
+  // last angle that corresponds to a
+  double u;
+  double v;
+  double w;
+};
+
+//==============================================================================
+// MGXS contains the mgxs data for a nuclide/material
+//==============================================================================
+
+class Mgxs {
+  private:
+
+    double_1dvec kTs;   // temperature in eV (k * T)
+    int scatter_format; // flag for if this is legendre, histogram, or tabular
+    int num_delayed_groups; // number of delayed neutron groups
+    int num_groups;     // number of energy groups
+    std::vector<XsData> xs; // Cross section data
+    // MGXS Incoming Flux Angular grid information
+    bool is_isotropic; // used to skip search for angle indices if isotropic
+    int n_pol;
+    int n_azi;
+    double_1dvec polar;
+    double_1dvec azimuthal;
+
+    //! \brief Initializes the Mgxs object metadata
+    //!
+    //! @param in_name Name of the object.
+    //! @param in_awr atomic-weight ratio.
+    //! @param in_kTs temperatures (in units of eV) that data is available.
+    //! @param in_fissionable Is this item fissionable or not.
+    //! @param in_scatter_format Denotes whether Legendre, Tabular, or
+    //!   Histogram scattering is used.
+    //! @param in_num_groups Number of energy groups.
+    //! @param in_num_delayed_groups Number of delayed groups.
+    //! @param in_is_isotropic Is this an isotropic or angular with respect to
+    //!   the incoming particle.
+    //! @param in_polar Polar angle grid.
+    //! @param in_azimuthal Azimuthal angle grid.
+    void
+    init(const std::string& in_name, double in_awr, const double_1dvec& in_kTs,
+         bool in_fissionable, int in_scatter_format, int in_num_groups,
+         int in_num_delayed_groups, bool in_is_isotropic,
+         const double_1dvec& in_polar, const double_1dvec& in_azimuthal);
+
+    //! \brief Initializes the Mgxs object metadata from the HDF5 file
+    //!
+    //! @param xs_id HDF5 group id for the cross section data.
+    //! @param in_num_groups Number of energy groups.
+    //! @param in_num_delayed_groups Number of delayed groups.
+    //! @param temperature Temperatures to read.
+    //! @param tolerance Tolerance of temperature selection method.
+    //! @param temps_to_read Resultant list of temperatures in the library
+    //!   to read which correspond to the requested temperatures.
+    //! @param order_dim Resultant dimensionality of the scattering order.
+    //! @param method Method of choosing nearest temperatures.
+    void
+    metadata_from_hdf5(hid_t xs_id, int in_num_groups,
+         int in_num_delayed_groups, const double_1dvec& temperature,
+         double tolerance, int_1dvec& temps_to_read, int& order_dim,
+         int& method);
+
+    //! \brief Performs the actual act of combining the microscopic data for a
+    //!   single temperature.
+    //!
+    //! @param micros Microscopic objects to combine.
+    //! @param scalars Scalars to multiply the microscopic data by.
+    //! @param micro_ts The temperature index of the microscopic objects that
+    //!   corresponds to the temperature of interest.
+    //! @param this_t The temperature index of the macroscopic object.
+    void
+    combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
+         const int_1dvec& micro_ts, int this_t);
+
+    //! \brief Checks to see if this and that are able to be combined
+    //!
+    //! This comparison is used when building macroscopic cross sections
+    //! from microscopic cross sections.
+    //! @param that The other Mgxs to compare to this one.
+    //! @return True if they can be combined, False otherwise.
+    bool equiv(const Mgxs& that);
+
+  public:
+
+    std::string name;   // name of dataset, e.g., UO2
+    double awr;         // atomic weight ratio
+    bool fissionable;   // Is this fissionable
+    std::vector<CacheData> cache; // index and data cache
+
+    Mgxs() = default;
+
+    //! \brief Constructor that loads the Mgxs object from the HDF5 file
+    //!
+    //! @param xs_id HDF5 group id for the cross section data.
+    //! @param energy_groups Number of energy groups.
+    //! @param delayed_groups Number of delayed groups.
+    //! @param temperature Temperatures to read.
+    //! @param tolerance Tolerance of temperature selection method.
+    //! @param max_order Maximum order requested by the user;
+    //!   this is only used for Legendre scattering.
+    //! @param legendre_to_tabular Flag to denote if any Legendre provided
+    //!   should be converted to a Tabular representation.
+    //! @param legendre_to_tabular_points If a conversion is requested, this
+    //!   provides the number of points to use in the tabular representation.
+    //! @param method Method of choosing nearest temperatures.
+    Mgxs(hid_t xs_id, int energy_groups,
+         int delayed_groups, const double_1dvec& temperature, double tolerance,
+         int max_order, bool legendre_to_tabular,
+         int legendre_to_tabular_points, int& method);
+
+    //! \brief Constructor that initializes and populates all data to build a
+    //!   macroscopic cross section from microscopic cross section.
+    //!
+    //! @param in_name Name of the object.
+    //! @param mat_kTs temperatures (in units of eV) that data is needed.
+    //! @param micros Microscopic objects to combine.
+    //! @param atom_densities Atom densities of those microscopic quantities.
+    //! @param tolerance Tolerance of temperature selection method.
+    //! @param method Method of choosing nearest temperatures.
+    Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
+         const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
+         double tolerance, int& method);
+
+    //! \brief Provides a cross section value given certain parameters
+    //!
+    //! @param xstype Type of cross section requested, according to the
+    //!   enumerated constants.
+    //! @param gin Incoming energy group.
+    //! @param gout Outgoing energy group; use nullptr if irrelevant, or if a
+    //!   sum is requested.
+    //! @param mu Cosine of the change-in-angle, for scattering quantities;
+    //!   use nullptr if irrelevant.
+    //! @param dg delayed group index; use nullptr if irrelevant.
+    //! @return Requested cross section value.
+    double
+    get_xs(int xstype, int gin, int* gout, double* mu, int* dg);
+
+    //! \brief Samples the fission neutron energy and if prompt or delayed.
+    //!
+    //! @param gin Incoming energy group.
+    //! @param dg Sampled delayed group index.
+    //! @param gout Sampled outgoing energy group.
+    void
+    sample_fission_energy(int gin, int& dg, int& gout);
+
+    //! \brief Samples the outgoing energy and angle from a scatter event.
+    //!
+    //! @param gin Incoming energy group.
+    //! @param gout Sampled outgoing energy group.
+    //! @param mu Sampled cosine of the change-in-angle.
+    //! @param wgt Weight of the particle to be adjusted.
+    void
+    sample_scatter(int gin, int& gout, double& mu, double& wgt);
+
+    //! \brief Calculates cross section quantities needed for tracking.
+    //!
+    //! @param gin Incoming energy group.
+    //! @param sqrtkT Temperature of the material.
+    //! @param uvw Incoming particle direction.
+    //! @param total_xs Resultant total cross section.
+    //! @param abs_xs Resultant absorption cross section.
+    //! @param nu_fiss_xs Resultant nu-fission cross section.
+    void
+    calculate_xs(int gin, double sqrtkT, const double uvw[3],
+         double& total_xs, double& abs_xs, double& nu_fiss_xs);
+
+    //! \brief Sets the temperature index in cache given a temperature
+    //!
+    //! @param sqrtkT Temperature of the material.
+    void
+    set_temperature_index(double sqrtkT);
+
+    //! \brief Sets the angle index in cache given a direction
+    //!
+    //! @param uvw Incoming particle direction.
+    void
+    set_angle_index(const double uvw[3]);
+};
+
+} // namespace openmc
+#endif // MGXS_H

src/mgxs_data.F90

@@ -1,5 +1,7 @@
 module mgxs_data
 
+  use, intrinsic :: ISO_C_BINDING
+
   use constants
   use algorithm,       only: find
   use dict_header,     only: DictCharInt
@@ -7,11 +9,11 @@ module mgxs_data
   use geometry_header, only: get_temperatures, cells
   use hdf5_interface
   use material_header, only: Material, materials, n_materials
-  use mgxs_header
+  use mgxs_interface
   use nuclide_header,  only: n_nuclides
   use set_header,      only: SetChar
   use settings
-  use stl_vector,      only: VectorReal
+  use stl_vector,      only: VectorReal, VectorChar
   use string,          only: to_lower
   implicit none
 
@@ -27,14 +29,11 @@ contains
     integer                 :: j              ! index over nuclides in material
     integer                 :: i_nuclide      ! index in nuclides array
     character(20)           :: name           ! name of library to load
-    integer                 :: representation ! Data representation
-    character(MAX_LINE_LEN) :: temp_str
     type(Material), pointer :: mat
     type(SetChar)           :: already_read
     integer(HID_T)          :: file_id
-    integer(HID_T)          :: xsdata_group
     logical                 :: file_exists
-    type(VectorReal), allocatable :: temps(:)
+    type(VectorReal), allocatable, target :: temps(:)
     character(MAX_WORD_LEN) :: word
     integer, allocatable    :: array(:)
 
@@ -69,9 +68,6 @@ contains
                        &version supported by OpenMC.")
     end if
 
-    ! allocate arrays for MGXS storage and cross section cache
-    allocate(nuclides_MG(n_nuclides))
-
     ! ==========================================================================
     ! READ ALL MGXS CROSS SECTION TABLES
 
@@ -80,89 +76,29 @@ contains
       mat => materials(i)
 
       NUCLIDE_LOOP: do j = 1, mat % n_nuclides
-        name = mat % names(j)
+        name = trim(mat % names(j)) // C_NULL_CHAR
+        i_nuclide = mat % nuclide(j)
 
         if (.not. already_read % contains(name)) then
-          i_nuclide = mat % nuclide(j)
+          call add_mgxs_c(file_id, name, num_energy_groups, num_delayed_groups, &
+               temps(i_nuclide) % size(), temps(i_nuclide) % data, &
+               temperature_tolerance, max_order, &
+               logical(legendre_to_tabular, C_BOOL), &
+               legendre_to_tabular_points, temperature_method)
 
-          call write_message("Loading " // trim(name) // " data...", 6)
-
-          ! Check to make sure cross section set exists in the library
-          if (object_exists(file_id, trim(name))) then
-            xsdata_group = open_group(file_id, trim(name))
-          else
-            call fatal_error("Data for '" // trim(name) // "' does not exist in "&
-                 &// trim(path_cross_sections))
-          end if
-
-          ! First find out the data representation
-          if (attribute_exists(xsdata_group, "representation")) then
-
-            call read_attribute(temp_str, xsdata_group, "representation")
-
-            if (trim(temp_str) == 'isotropic') then
-              representation = MGXS_ISOTROPIC
-            else if (trim(temp_str) == 'angle') then
-              representation = MGXS_ANGLE
-            else
-              call fatal_error("Invalid Data Representation!")
-            end if
-          else
-            ! Default to isotropic representation
-            representation = MGXS_ISOTROPIC
-          end if
-
-          ! Now allocate accordingly
-          select case(representation)
-
-          case(MGXS_ISOTROPIC)
-            allocate(MgxsIso :: nuclides_MG(i_nuclide) % obj)
-
-          case(MGXS_ANGLE)
-            allocate(MgxsAngle :: nuclides_MG(i_nuclide) % obj)
-
-          end select
-
-          ! Now read in the data specific to the type we just declared
-          call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, &
-               num_energy_groups, num_delayed_groups, temps(i_nuclide), &
-               temperature_method, temperature_tolerance, max_order, &
-               legendre_to_tabular, legendre_to_tabular_points)
-
-          ! Add name to dictionary
           call already_read % add(name)
-
-          call close_group(xsdata_group)
-
         end if
       end do NUCLIDE_LOOP
+
+      mat % fissionable = query_fissionable_c(mat % n_nuclides, mat % nuclide)
+
     end do MATERIAL_LOOP
 
+    call file_close(file_id)
+
     ! Avoid some valgrind leak errors
     call already_read % clear()
 
-    ! Loop around material
-    MATERIAL_LOOP3: do i = 1, n_materials
-
-      ! Get material
-      mat => materials(i)
-
-      ! Loop around nuclides in material
-      NUCLIDE_LOOP2: do j = 1, mat % n_nuclides
-
-        ! Is this fissionable?
-        if (nuclides_MG(mat % nuclide(j)) % obj % fissionable) then
-          mat % fissionable = .true.
-        end if
-        if (mat % fissionable) then
-          exit NUCLIDE_LOOP2
-        end if
-
-      end do NUCLIDE_LOOP2
-    end do MATERIAL_LOOP3
-
-    call file_close(file_id)
-
   end subroutine read_mgxs
 
 !===============================================================================
@@ -173,8 +109,7 @@ contains
     integer                       :: i_mat ! index in materials array
     type(Material), pointer       :: mat   ! current material
     type(VectorReal), allocatable :: kTs(:)
-
-    allocate(macro_xs(n_materials))
+    character(MAX_WORD_LEN)       :: name  ! name of material
 
     ! Get temperatures to read for each material
     call get_mat_kTs(kTs)
@@ -188,19 +123,13 @@ contains
       ! Get the material
       mat => materials(i_mat)
 
-      ! Get the scattering type for the first nuclide
-      select type(nuc => nuclides_MG(mat % nuclide(1)) % obj)
-      type is (MgxsIso)
-        allocate(MgxsIso :: macro_xs(i_mat) % obj)
-      type is (MgxsAngle)
-        allocate(MgxsAngle :: macro_xs(i_mat) % obj)
-      end select
+      name = trim(mat % name) // C_NULL_CHAR
 
       ! Do not read materials which we do not actually use in the problem to
       ! reduce storage
       if (allocated(kTs(i_mat) % data)) then
-        call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, &
-             num_energy_groups, num_delayed_groups, max_order, &
+        call create_macro_xs_c(name, mat % n_nuclides, mat % nuclide, &
+             kTs(i_mat) % size(), kTs(i_mat) % data, mat % atom_density, &
              temperature_tolerance, temperature_method)
       end if
     end do

src/mgxs_interface.F90

@@ -0,0 +1,164 @@
+module mgxs_interface
+
+  use, intrinsic :: ISO_C_BINDING
+
+  use hdf5_interface
+
+  implicit none
+
+  interface
+
+    subroutine add_mgxs_c(file_id, name, energy_groups, delayed_groups, &
+         n_temps, temps, tolerance, max_order, legendre_to_tabular, &
+         legendre_to_tabular_points, method) bind(C)
+      use ISO_C_BINDING
+      import HID_T
+      implicit none
+      integer(HID_T), value, intent(in) :: file_id
+      character(kind=C_CHAR),intent(in) :: name(*)
+      integer(C_INT), value, intent(in) :: energy_groups
+      integer(C_INT), value, intent(in) :: delayed_groups
+      integer(C_INT), value, intent(in) :: n_temps
+      real(C_DOUBLE),        intent(in) :: temps(1:n_temps)
+      real(C_DOUBLE), value, intent(in) :: tolerance
+      integer(C_INT), value, intent(in) :: max_order
+      logical(C_BOOL),value, intent(in) :: legendre_to_tabular
+      integer(C_INT), value, intent(in) :: legendre_to_tabular_points
+      integer(C_INT),     intent(inout) :: method
+    end subroutine add_mgxs_c
+
+    function query_fissionable_c(n_nuclides, i_nuclides) result(result) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n_nuclides
+      integer(C_INT),        intent(in) :: i_nuclides(1:n_nuclides)
+      logical(C_BOOL)                   :: result
+    end function query_fissionable_c
+
+    subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, &
+         atom_densities, tolerance, method) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      character(kind=C_CHAR),intent(in) :: name(*)
+      integer(C_INT), value, intent(in) :: n_nuclides
+      integer(C_INT),        intent(in) :: i_nuclides(1:n_nuclides)
+      integer(C_INT), value, intent(in) :: n_temps
+      real(C_DOUBLE),        intent(in) :: temps(1:n_temps)
+      real(C_DOUBLE),        intent(in) :: atom_densities(1:n_nuclides)
+      real(C_DOUBLE), value, intent(in) :: tolerance
+      integer(C_INT),     intent(inout) :: method
+    end subroutine create_macro_xs_c
+
+    subroutine calculate_xs_c(i_mat, gin, sqrtkT, uvw, total_xs, abs_xs, &
+         nu_fiss_xs) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: i_mat
+      integer(C_INT), value, intent(in) :: gin
+      real(C_DOUBLE), value, intent(in) :: sqrtkT
+      real(C_DOUBLE),        intent(in) :: uvw(1:3)
+      real(C_DOUBLE),     intent(inout) :: total_xs
+      real(C_DOUBLE),     intent(inout) :: abs_xs
+      real(C_DOUBLE),     intent(inout) :: nu_fiss_xs
+    end subroutine calculate_xs_c
+
+    subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: i_mat
+      integer(C_INT), value, intent(in) :: gin
+      integer(C_INT),     intent(inout) :: gout
+      real(C_DOUBLE),     intent(inout) :: mu
+      real(C_DOUBLE),     intent(inout) :: wgt
+      real(C_DOUBLE),     intent(inout) :: uvw(1:3)
+    end subroutine sample_scatter_c
+
+    subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: i_mat
+      integer(C_INT), value, intent(in) :: gin
+      integer(C_INT),     intent(inout) :: dg
+      integer(C_INT),     intent(inout) :: gout
+    end subroutine sample_fission_energy_c
+
+    subroutine get_name_c(index, name_len, name) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value,  intent(in)    :: index
+      integer(C_INT), value,  intent(in)    :: name_len
+      character(kind=C_CHAR), intent(inout) :: name(name_len)
+    end subroutine get_name_c
+
+    function get_awr_c(index) result(awr) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: index
+      real(C_DOUBLE)                    :: awr
+    end function get_awr_c
+
+    function get_nuclide_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
+         bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT),  value,   intent(in) :: index
+      integer(C_INT),  value,   intent(in) :: xstype
+      integer(C_INT),  value,   intent(in) :: gin
+      integer(C_INT), optional, intent(in) :: gout
+      real(C_DOUBLE), optional, intent(in) :: mu
+      integer(C_INT), optional, intent(in) :: dg
+      real(C_DOUBLE)                       :: val
+    end function get_nuclide_xs_c
+
+    function get_macro_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
+         bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT),  value,   intent(in) :: index
+      integer(C_INT),  value,   intent(in) :: xstype
+      integer(C_INT),  value,   intent(in) :: gin
+      integer(C_INT), optional, intent(in) :: gout
+      real(C_DOUBLE), optional, intent(in) :: mu
+      integer(C_INT), optional, intent(in) :: dg
+      real(C_DOUBLE)                       :: val
+    end function get_macro_xs_c
+
+    subroutine set_nuclide_angle_index_c(index, uvw) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: index
+      real(C_DOUBLE),        intent(in) :: uvw(1:3)
+    end subroutine set_nuclide_angle_index_c
+
+    subroutine set_macro_angle_index_c(index, uvw) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: index
+      real(C_DOUBLE),        intent(in) :: uvw(1:3)
+    end subroutine set_macro_angle_index_c
+
+    subroutine set_nuclide_temperature_index_c(index, sqrtkT) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: index
+      real(C_DOUBLE), value, intent(in) :: sqrtkT
+    end subroutine set_nuclide_temperature_index_c
+
+  end interface
+
+  ! Number of energy groups
+  integer(C_INT) :: num_energy_groups
+
+  ! Number of delayed groups
+  integer(C_INT) :: num_delayed_groups
+
+  ! Energy group structure with decreasing energy
+  real(8), allocatable :: energy_bins(:)
+
+  ! Midpoint of the energy group structure
+  real(8), allocatable :: energy_bin_avg(:)
+
+  ! Energy group structure with increasing energy
+  real(8), allocatable :: rev_energy_bins(:)
+
+end module mgxs_interface

src/mgxs_interface.cpp

@@ -0,0 +1,229 @@
+#include <string>
+
+#include "error.h"
+#include "math_functions.h"
+#include "mgxs_interface.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// Mgxs data loading interface methods
+//==============================================================================
+
+void
+add_mgxs_c(hid_t file_id, const char* name, int energy_groups,
+     int delayed_groups, int n_temps, const double temps[],  double tolerance,
+     int max_order, bool legendre_to_tabular, int legendre_to_tabular_points,
+     int& method)
+{
+  // Convert temps to a vector for the from_hdf5 function
+  double_1dvec temperature(temps, temps + n_temps);
+
+  write_message("Loading " + std::string(name) + " data...", 6);
+
+  // Check to make sure cross section set exists in the library
+  hid_t xs_grp;
+  if (object_exists(file_id, name)) {
+    xs_grp = open_group(file_id, name);
+  } else {
+    fatal_error("Data for " + std::string(name) + " does not exist in "
+                + "provided MGXS Library");
+  }
+
+  Mgxs mg(xs_grp, energy_groups, delayed_groups, temperature, tolerance,
+          max_order, legendre_to_tabular, legendre_to_tabular_points, method);
+
+  nuclides_MG.push_back(mg);
+  close_group(xs_grp);
+}
+
+//==============================================================================
+
+bool
+query_fissionable_c(int n_nuclides, const int i_nuclides[])
+{
+  bool result = false;
+  for (int n = 0; n < n_nuclides; n++) {
+    if (nuclides_MG[i_nuclides[n] - 1].fissionable) result = true;
+  }
+  return result;
+}
+
+//==============================================================================
+
+void
+create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
+     int n_temps, const double temps[], const double atom_densities[],
+     double tolerance, int& method)
+{
+  if (n_temps > 0) {
+    // // Convert temps to a vector
+    double_1dvec temperature(temps, temps + n_temps);
+
+    // Convert atom_densities to a vector
+    double_1dvec atom_densities_vec(atom_densities,
+         atom_densities + n_nuclides);
+
+    // Build array of pointers to nuclides_MG's Mgxs objects needed for this
+    // material
+    std::vector<Mgxs*> mgxs_ptr(n_nuclides);
+    for (int n = 0; n < n_nuclides; n++) {
+      mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1];
+    }
+
+    Mgxs macro(mat_name, temperature, mgxs_ptr, atom_densities_vec,
+         tolerance, method);
+    macro_xs.emplace_back(macro);
+  } else {
+    // Preserve the ordering of materials by including a blank entry
+    Mgxs macro;
+    macro_xs.emplace_back(macro);
+  }
+}
+
+//==============================================================================
+// Mgxs tracking/transport/tallying interface methods
+//==============================================================================
+
+void
+calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
+     double& total_xs, double& abs_xs, double& nu_fiss_xs)
+{
+  macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
+       nu_fiss_xs);
+}
+
+//==============================================================================
+
+void
+sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
+                 double uvw[3])
+{
+  int gout_c = gout - 1;
+  macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
+
+  // adjust return value for fortran indexing
+  gout = gout_c + 1;
+
+  // Rotate the angle
+  rotate_angle_c(uvw, mu, nullptr);
+}
+
+//==============================================================================
+
+void
+sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout)
+{
+  int dg_c = 0;
+  int gout_c = 0;
+  macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
+
+  // adjust return values for fortran indexing
+  dg = dg_c + 1;
+  gout = gout_c + 1;
+}
+
+//==============================================================================
+
+double
+get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
+{
+  int gout_c;
+  int* gout_c_p;
+  int dg_c;
+  int* dg_c_p;
+  if (gout != nullptr) {
+    gout_c = *gout - 1;
+    gout_c_p = &gout_c;
+  } else {
+    gout_c_p = gout;
+  }
+  if (dg != nullptr) {
+    dg_c = *dg - 1;
+    dg_c_p = &dg_c;
+  } else {
+    dg_c_p = dg;
+  }
+  return nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
+}
+
+//==============================================================================
+
+double
+get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
+{
+  int gout_c;
+  int* gout_c_p;
+  int dg_c;
+  int* dg_c_p;
+  if (gout != nullptr) {
+    gout_c = *gout - 1;
+    gout_c_p = &gout_c;
+  } else {
+    gout_c_p = gout;
+  }
+  if (dg != nullptr) {
+    dg_c = *dg - 1;
+    dg_c_p = &dg_c;
+  } else {
+    dg_c_p = dg;
+  }
+  return macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
+}
+
+//==============================================================================
+
+void
+set_nuclide_angle_index_c(int index, const double uvw[3])
+{
+  // Update the values
+  nuclides_MG[index - 1].set_angle_index(uvw);
+}
+
+//==============================================================================
+
+void
+set_macro_angle_index_c(int index, const double uvw[3])
+{
+  // Update the values
+  macro_xs[index - 1].set_angle_index(uvw);
+}
+
+//==============================================================================
+
+void
+set_nuclide_temperature_index_c(int index, double sqrtkT)
+{
+  // Update the values
+  nuclides_MG[index - 1].set_temperature_index(sqrtkT);
+}
+
+//==============================================================================
+// General Mgxs methods
+//==============================================================================
+
+void
+get_name_c(int index, int name_len, char* name)
+{
+  // First blank out our input string
+  std::string str(name_len, ' ');
+  std::strcpy(name, str.c_str());
+
+  // Now get the data and copy to the C-string
+  str = nuclides_MG[index - 1].name;
+  std::strcpy(name, str.c_str());
+
+  // Finally, remove the null terminator
+  name[std::strlen(name)] = ' ';
+}
+
+//==============================================================================
+
+double
+get_awr_c(int index)
+{
+  return nuclides_MG[index - 1].awr;
+}
+
+} // namespace openmc

src/mgxs_interface.h

@@ -0,0 +1,79 @@
+//! \file mgxs_interface.h
+//! A collection of C interfaces to the C++ Mgxs class
+
+#ifndef MGXS_INTERFACE_H
+#define MGXS_INTERFACE_H
+
+#include "hdf5_interface.h"
+#include "mgxs.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+extern std::vector<Mgxs> nuclides_MG;
+extern std::vector<Mgxs> macro_xs;
+
+//==============================================================================
+// Mgxs data loading interface methods
+//==============================================================================
+
+extern "C" void
+add_mgxs_c(hid_t file_id, const char* name, int energy_groups,
+     int delayed_groups, int n_temps, const double temps[], double tolerance,
+     int max_order, bool legendre_to_tabular, int legendre_to_tabular_points,
+     int& method);
+
+extern "C" bool
+query_fissionable_c(int n_nuclides, const int i_nuclides[]);
+
+extern "C" void
+create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
+     int n_temps, const double temps[], const double atom_densities[],
+     double tolerance, int& method);
+
+//==============================================================================
+// Mgxs tracking/transport/tallying interface methods
+//==============================================================================
+
+extern "C" void
+calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
+     double& total_xs, double& abs_xs, double& nu_fiss_xs);
+
+extern "C" void
+sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
+     double uvw[3]);
+
+extern "C" void
+sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout);
+
+extern "C" double
+get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
+
+extern "C" double
+get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
+
+extern "C" void
+set_nuclide_angle_index_c(int index, const double uvw[3]);
+
+extern "C" void
+set_macro_angle_index_c(int index, const double uvw[3]);
+
+extern "C" void
+set_nuclide_temperature_index_c(int index, double sqrtkT);
+
+//==============================================================================
+// General Mgxs methods
+//==============================================================================
+
+extern "C" void
+get_name_c(int index, int name_len, char* name);
+
+extern "C" double
+get_awr_c(int index);
+
+} // namespace openmc
+#endif // MGXS_INTERFACE_H

src/nuclide_header.F90

@@ -165,10 +165,10 @@ module nuclide_header
 !===============================================================================
 
   type MaterialMacroXS
-    real(8) :: total         ! macroscopic total xs
-    real(8) :: absorption    ! macroscopic absorption xs
-    real(8) :: fission       ! macroscopic fission xs
-    real(8) :: nu_fission    ! macroscopic production xs
+    real(C_DOUBLE) :: total         ! macroscopic total xs
+    real(C_DOUBLE) :: absorption    ! macroscopic absorption xs
+    real(C_DOUBLE) :: fission       ! macroscopic fission xs
+    real(C_DOUBLE) :: nu_fission    ! macroscopic production xs
   end type MaterialMacroXS
 
 !===============================================================================

src/output.F90

@@ -12,7 +12,7 @@ module output
   use math,            only: t_percentile
   use mesh_header,     only: RegularMesh, meshes
   use message_passing, only: master, n_procs
-  use mgxs_header,     only: nuclides_MG
+  use mgxs_interface
   use nuclide_header
   use particle_header, only: LocalCoord, Particle
   use plot_header
@@ -680,6 +680,7 @@ contains
     character(36)           :: score_name                  ! names of scoring function
                                                            ! to be applied at write-time
     type(TallyFilterMatch), allocatable :: matches(:)
+    character(MAX_WORD_LEN) :: temp_name
 
     ! Skip if there are no tallies
     if (n_tallies == 0) return
@@ -843,8 +844,9 @@ contains
               write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
                    trim(nuclides(i_nuclide) % name)
             else
+              call get_name_c(i_nuclide, len(temp_name), temp_name)
               write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
-                   trim(nuclides_MG(i_nuclide) % obj % name)
+                   trim(temp_name)
             end if
           end if
 

src/particle_restart.F90

@@ -6,7 +6,7 @@ module particle_restart
   use constants
   use error,            only: write_message
   use hdf5_interface,   only: file_open, file_close, read_dataset, HID_T
-  use mgxs_header,      only: energy_bin_avg
+  use mgxs_interface,   only: energy_bin_avg
   use nuclide_header,   only: micro_xs, n_nuclides
   use output,           only: print_particle
   use particle_header,  only: Particle

src/physics_mg.F90

@@ -8,13 +8,12 @@ module physics_mg
   use material_header,        only: Material, materials
   use math,                   only: rotate_angle
   use mesh_header,            only: meshes
-  use mgxs_header
+  use mgxs_interface
   use message_passing
   use nuclide_header,         only: material_xs
   use particle_header,        only: Particle
   use physics_common
   use random_lcg,             only: prn
-  use scattdata_header
   use settings
   use simulation_header
   use string,                 only: to_str
@@ -143,16 +142,12 @@ contains
 
     type(Particle), intent(inout)  :: p
 
-    call macro_xs(p % material) % obj % sample_scatter(p % coord(1) % uvw, &
-                                                       p % last_g, p % g, &
-                                                       p % mu, p % wgt)
+    call sample_scatter_c(p % material, p % last_g, p % g, p % mu, &
+         p % wgt, p % coord(1) % uvw)
 
     ! Update energy value for downstream compatability (in tallying)
     p % E = energy_bin_avg(p % g)
 
-    ! Convert change in angle (mu) to new direction
-    p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu)
-
     ! Set event component
     p % event = EVENT_SCATTER
 
@@ -178,10 +173,6 @@ contains
     real(8) :: mu                       ! fission neutron angular cosine
     real(8) :: phi                      ! fission neutron azimuthal angle
     real(8) :: weight                   ! weight adjustment for ufs method
-    class(Mgxs), pointer :: xs
-
-    ! Get Pointers
-    xs => macro_xs(p % material) % obj
 
     ! TODO: Heat generation from fission
 
@@ -261,7 +252,7 @@ contains
 
       ! Sample secondary energy distribution for fission reaction and set energy
       ! in fission bank
-      call xs % sample_fission_energy(p % g, bank_array(i) % uvw, dg, gout)
+      call sample_fission_energy_c(p % material, p % g, dg, gout)
 
       bank_array(i) % E = real(gout, 8)
       bank_array(i) % delayed_group = dg

src/scattdata.cpp

@@ -0,0 +1,936 @@
+#include <algorithm>
+#include <numeric>
+#include <cmath>
+
+#include "constants.h"
+#include "math_functions.h"
+#include "random_lcg.h"
+#include "error.h"
+#include "scattdata.h"
+
+namespace openmc {
+
+//==============================================================================
+// ScattData base-class methods
+//==============================================================================
+
+void
+ScattData::base_init(int order, const int_1dvec& in_gmin,
+     const int_1dvec& in_gmax, const double_2dvec& in_energy,
+     const double_2dvec& in_mult)
+{
+  int groups = in_energy.size();
+
+  gmin = in_gmin;
+  gmax = in_gmax;
+  energy.resize(groups);
+  mult.resize(groups);
+  dist.resize(groups);
+
+  for (int gin = 0; gin < groups; gin++) {
+    // Store the inputted data
+    energy[gin] = in_energy[gin];
+    mult[gin] = in_mult[gin];
+
+    // Make sure the energy is normalized
+    double norm = std::accumulate(energy[gin].begin(), energy[gin].end(), 0.);
+
+    if (norm != 0.) {
+      for (auto& n : energy[gin]) n /= norm;
+    }
+
+    // Initialize the distribution data
+    dist[gin].resize(in_gmax[gin] - in_gmin[gin] + 1);
+    for (auto& v : dist[gin]) {
+      v.resize(order);
+    }
+  }
+}
+
+//==============================================================================
+
+void
+ScattData::base_combine(int max_order,
+     const std::vector<ScattData*>& those_scatts, const double_1dvec& scalars,
+     int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& sparse_mult,
+     double_3dvec& sparse_scatter)
+{
+  int groups = those_scatts[0] -> energy.size();
+
+  // Now allocate and zero our storage spaces
+  double_3dvec this_matrix = double_3dvec(groups, double_2dvec(groups,
+       double_1dvec(max_order, 0.)));
+  double_2dvec mult_numer(groups, double_1dvec(groups, 0.));
+  double_2dvec mult_denom(groups, double_1dvec(groups, 0.));
+
+  // Build the dense scattering and multiplicity matrices
+  // Get the multiplicity_matrix
+  // To combine from nuclidic data we need to use the final relationship
+  // mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
+  //              sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
+  // Developed as follows:
+  // mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'}
+  // mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'})
+  // mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
+  //              sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
+  // nuscatt_{i,g,g'} can be reconstructed from the energy and scattxs member
+  // variables
+  for (int i = 0; i < those_scatts.size(); i++) {
+    ScattData* that = those_scatts[i];
+
+    // Build the dense matrix for that object
+    double_3dvec that_matrix = that->get_matrix(max_order);
+
+    // Now add that to this for the scattering and multiplicity
+    for (int gin = 0; gin < groups; gin++) {
+      // Only spend time adding that's gmin to gmax data since the rest will
+      // be zeros
+      int i_gout = 0;
+      for (int gout = that->gmin[gin]; gout <= that->gmax[gin]; gout++) {
+        // Do the scattering matrix
+        for (int l = 0; l < max_order; l++) {
+          this_matrix[gin][gout][l] += scalars[i] * that_matrix[gin][gout][l];
+        }
+
+        // Incorporate that's contribution to the multiplicity matrix data
+        double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
+        mult_numer[gin][gout] += scalars[i] * nuscatt;
+        if (that->mult[gin][i_gout] > 0.) {
+          mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
+        } else {
+          mult_denom[gin][gout] += scalars[i];
+        }
+        i_gout++;
+      }
+    }
+  }
+
+  // Combine mult_numer and mult_denom into the combined multiplicity matrix
+  double_2dvec this_mult(groups, double_1dvec(groups, 1.));
+  for (int gin = 0; gin < groups; gin++) {
+    for (int gout = 0; gout < groups; gout++) {
+      if (mult_denom[gin][gout] > 0.) {
+        this_mult[gin][gout] = mult_numer[gin][gout] / mult_denom[gin][gout];
+      }
+    }
+  }
+  mult_numer.clear();
+  mult_denom.clear();
+
+  // We have the data, now we need to convert to a jagged array and then use
+  // the initialize function to store it on the object.
+  for (int gin = 0; gin < groups; gin++) {
+    // Find the minimum and maximum group boundaries
+    int gmin_;
+    for (gmin_ = 0; gmin_ < groups; gmin_++) {
+      bool non_zero = false;
+      for (int l = 0; l < this_matrix[gin][gmin_].size(); l++) {
+        if (this_matrix[gin][gmin_][l] != 0.) {
+          non_zero = true;
+          break;
+        }
+      }
+      if (non_zero) break;
+    }
+    int gmax_;
+    for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) {
+      bool non_zero = false;
+      for (int l = 0; l < this_matrix[gin][gmax_].size(); l++) {
+        if (this_matrix[gin][gmax_][l] != 0.) {
+          non_zero = true;
+          break;
+        }
+      }
+      if (non_zero) break;
+    }
+
+    // treat the case of all values being 0
+    if (gmin_ > gmax_) {
+      gmin_ = gin;
+      gmax_ = gin;
+    }
+
+    // Store the group bounds
+    in_gmin[gin] = gmin_;
+    in_gmax[gin] = gmax_;
+
+    // Store the data in the compressed format
+    sparse_scatter[gin].resize(gmax_ - gmin_ + 1);
+    sparse_mult[gin].resize(gmax_ - gmin_ + 1);
+    int i_gout = 0;
+    for (int gout = gmin_; gout <= gmax_; gout++) {
+      sparse_scatter[gin][i_gout] = this_matrix[gin][gout];
+      sparse_mult[gin][i_gout] = this_mult[gin][gout];
+      i_gout++;
+    }
+  }
+}
+
+//==============================================================================
+
+void
+ScattData::sample_energy(int gin, int& gout, int& i_gout)
+{
+  // Sample the outgoing group
+  double xi = prn();
+
+  i_gout = 0;
+  gout = gmin[gin];
+  double prob = energy[gin][i_gout];
+  while((prob < xi) && (gout < gmax[gin])) {
+    gout++;
+    i_gout++;
+    prob += energy[gin][i_gout];
+  }
+}
+
+//==============================================================================
+
+double
+ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu)
+{
+  // Set the outgoing group offset index as needed
+  int i_gout = 0;
+  if (gout != nullptr) {
+    // short circuit the function if gout is from a zero portion of the
+    // scattering matrix
+    if ((*gout < gmin[gin]) || (*gout > gmax[gin])) { // > gmax?
+      return 0.;
+    }
+    i_gout = *gout - gmin[gin];
+  }
+
+  double val = scattxs[gin];
+  switch(xstype) {
+  case MG_GET_XS_SCATTER:
+    if (gout != nullptr) val *= energy[gin][i_gout];
+    break;
+  case MG_GET_XS_SCATTER_MULT:
+    if (gout != nullptr) {
+      val *= energy[gin][i_gout] / mult[gin][i_gout];
+    } else {
+      val /= std::inner_product(mult[gin].begin(), mult[gin].end(),
+                                energy[gin].begin(), 0.0);
+    }
+    break;
+  case MG_GET_XS_SCATTER_FMU_MULT:
+    if ((gout != nullptr) && (mu != nullptr)) {
+      val *= energy[gin][i_gout] * calc_f(gin, *gout, *mu);
+    } else {
+      // This is not an expected path (asking for f_mu without asking for a
+      // group or mu is not useful
+      fatal_error("Invalid call to get_xs");
+    }
+    break;
+  case MG_GET_XS_SCATTER_FMU:
+    if ((gout != nullptr) && (mu != nullptr)) {
+      val *= energy[gin][i_gout] * calc_f(gin, *gout, *mu) / mult[gin][i_gout];
+    } else {
+      // This is not an expected path (asking for f_mu without asking for a
+      // group or mu is not useful
+      fatal_error("Invalid call to get_xs");
+    }
+    break;
+  }
+  return val;
+}
+
+//==============================================================================
+// ScattDataLegendre methods
+//==============================================================================
+
+void
+ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
+     const double_2dvec& in_mult, const double_3dvec& coeffs)
+{
+  int groups = coeffs.size();
+  int order = coeffs[0][0].size();
+
+  // make a copy of coeffs that we can use to both extract data and normalize
+  double_3dvec matrix = coeffs;
+
+  // Get the scattering cross section value by summing the un-normalized P0
+  // coefficient in the variable matrix over all outgoing groups.
+  scattxs.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
+    scattxs[gin] = 0.;
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      scattxs[gin] += matrix[gin][i_gout][0];
+    }
+  }
+
+  // Build the energy transfer matrix from data in the variable matrix while
+  // also normalizing the variable matrix itself
+  // (forcing the CDF of f(mu=1) == 1)
+  double_2dvec in_energy;
+  in_energy.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
+    in_energy[gin].resize(num_groups);
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      double norm = matrix[gin][i_gout][0];
+      in_energy[gin][i_gout] = norm;
+      if (norm != 0.) {
+        for (auto& n : matrix[gin][i_gout]) n /= norm;
+      }
+    }
+  }
+
+  // Initialize the base class attributes
+  ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult);
+
+  // Set the distribution (sdata.dist) values and initialize max_val
+  max_val.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = gmax[gin] - gmin[gin] + 1;
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      dist[gin][i_gout] = matrix[gin][i_gout];
+    }
+    max_val[gin].resize(num_groups);
+    for (auto& n : max_val[gin]) n = 0.;
+  }
+
+  // Now update the maximum value
+  update_max_val();
+}
+
+//==============================================================================
+
+void
+ScattDataLegendre::update_max_val()
+{
+  int groups = max_val.size();
+  // Step through the polynomial with fixed number of points to identify the
+  // maximal value
+  int Nmu = 1001;
+  double dmu = 2. / (Nmu - 1);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = gmax[gin] - gmin[gin] + 1;
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      for (int imu = 0; imu < Nmu; imu++) {
+        double mu;
+        if (imu == 0) {
+          mu = -1.;
+        } else if (imu == (Nmu - 1)) {
+          mu = 1.;
+        } else {
+          mu = -1. + (imu - 1) * dmu;
+        }
+
+        // Calculate probability
+        double f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
+                                       dist[gin][i_gout].data(), mu);
+
+        // if this is a new maximum, store it
+        if (f > max_val[gin][i_gout]) max_val[gin][i_gout] = f;
+      } // end imu loop
+
+      // Since we may not have caught the true max, add 10% margin
+      max_val[gin][i_gout] *= 1.1;
+    }
+  }
+}
+
+//==============================================================================
+
+double
+ScattDataLegendre::calc_f(int gin, int gout, double mu)
+{
+  double f;
+  if ((gout < gmin[gin]) || (gout > gmax[gin])) {
+    f = 0.;
+  } else {
+    int i_gout = gout - gmin[gin];
+    f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
+                            dist[gin][i_gout].data(), mu);
+  }
+  return f;
+}
+
+//==============================================================================
+
+void
+ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt)
+{
+  // Sample the outgoing energy using the base-class method
+  int i_gout;
+  sample_energy(gin, gout, i_gout);
+
+  // Now we can sample mu using the scattering kernel using rejection
+  // sampling from a rectangular bounding box
+  double M = max_val[gin][i_gout];
+  int samples = 0;
+
+  while(true) {
+    mu = 2. * prn() - 1.;
+    double f = calc_f(gin, gout, mu);
+    if (f > 0.) {
+      double u = prn() * M;
+      if (u <= f) break;
+    }
+    samples++;
+    if (samples > MAX_SAMPLE) {
+        fatal_error("Maximum number of Legendre expansion samples reached");
+    }
+  };
+
+  // Update the weight to reflect neutron multiplicity
+  wgt *= mult[gin][i_gout];
+}
+
+//==============================================================================
+
+void
+ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
+                           const double_1dvec& scalars)
+{
+  // Find the max order in the data set and make sure we can combine the sets
+  int max_order = 0;
+  for (int i = 0; i < those_scatts.size(); i++) {
+    // Lets also make sure these items are combineable
+    ScattDataLegendre* that = dynamic_cast<ScattDataLegendre*>(those_scatts[i]);
+    if (!that) {
+      fatal_error("Cannot combine the ScattData objects!");
+    }
+    int that_order = that->get_order();
+    if (that_order > max_order) max_order = that_order;
+  }
+  max_order++;  // Add one since this is a Legendre
+
+  int groups = those_scatts[0] -> energy.size();
+
+  int_1dvec in_gmin(groups);
+  int_1dvec in_gmax(groups);
+  double_3dvec sparse_scatter(groups);
+  double_2dvec sparse_mult(groups);
+
+  // The rest of the steps do not depend on the type of angular representation
+  // so we use a base class method to sum up xs and create new energy and mult
+  // matrices
+  ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
+                          sparse_mult, sparse_scatter);
+
+  // Got everything we need, store it.
+  init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
+}
+
+//==============================================================================
+
+double_3dvec
+ScattDataLegendre::get_matrix(int max_order)
+{
+  // Get the sizes and initialize the data to 0
+  int groups = energy.size();
+  int order_dim = max_order + 1;
+  double_3dvec matrix = double_3dvec(groups, double_2dvec(groups,
+       double_1dvec(order_dim, 0.)));
+
+  for (int gin = 0; gin < groups; gin++) {
+    for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) {
+      int gout = i_gout + gmin[gin];
+      for (int l = 0; l < order_dim; l++) {
+        matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] *
+             dist[gin][i_gout][l];
+      }
+    }
+  }
+  return matrix;
+}
+
+//==============================================================================
+// ScattDataHistogram methods
+//==============================================================================
+
+void
+ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
+     const double_2dvec& in_mult, const double_3dvec& coeffs)
+{
+  int groups = coeffs.size();
+  int order = coeffs[0][0].size();
+
+  // make a copy of coeffs that we can use to both extract data and normalize
+  double_3dvec matrix = coeffs;
+
+  // Get the scattering cross section value by summing the distribution
+  // over all the histogram bins in angle and outgoing energy groups
+  scattxs.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    scattxs[gin] = 0.;
+    for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) {
+      scattxs[gin] += std::accumulate(matrix[gin][i_gout].begin(),
+                                      matrix[gin][i_gout].end(), 0.);
+    }
+  }
+
+  // Build the energy transfer matrix from data in the variable matrix
+  double_2dvec in_energy;
+  in_energy.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
+    in_energy[gin].resize(num_groups);
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      double norm = std::accumulate(matrix[gin][i_gout].begin(),
+                                    matrix[gin][i_gout].end(), 0.);
+      in_energy[gin][i_gout] = norm;
+      if (norm != 0.) {
+        for (auto& n : matrix[gin][i_gout]) n /= norm;
+      }
+    }
+  }
+
+  // Initialize the base class attributes
+  ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult);
+
+  // Build the angular distribution mu values
+  mu = double_1dvec(order);
+  dmu = 2. / order;
+  mu[0] = -1.;
+  for (int imu = 1; imu < order; imu++) {
+    mu[imu] = -1. + imu * dmu;
+  }
+
+  // Calculate f(mu) and integrate it so we can avoid rejection sampling
+  fmu.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = gmax[gin] - gmin[gin] + 1;
+    fmu[gin].resize(num_groups);
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      fmu[gin][i_gout].resize(order);
+      // The variable matrix contains f(mu); so directly assign it
+      fmu[gin][i_gout] = matrix[gin][i_gout];
+
+      // Integrate the histogram
+      dist[gin][i_gout][0] = dmu * matrix[gin][i_gout][0];
+      for (int imu = 1; imu < order; imu++) {
+        dist[gin][i_gout][imu] = dmu * matrix[gin][i_gout][imu] +
+             dist[gin][i_gout][imu - 1];
+      }
+
+      // Now re-normalize for integral to unity
+      double norm = dist[gin][i_gout][order - 1];
+      if (norm > 0.) {
+        for (int imu = 0; imu < order; imu++) {
+          fmu[gin][i_gout][imu] /= norm;
+          dist[gin][i_gout][imu] /= norm;
+        }
+      }
+    }
+  }
+}
+
+//==============================================================================
+
+double
+ScattDataHistogram::calc_f(int gin, int gout, double mu)
+{
+  double f;
+  if ((gout < gmin[gin]) || (gout > gmax[gin])) {
+    f = 0.;
+  } else {
+    // Find mu bin
+    int i_gout = gout - gmin[gin];
+    int imu;
+    if (mu == 1.) {
+      // use size -2 to have the index one before the end
+      imu = this->mu.size() - 2;
+    } else {
+      imu = std::floor((mu + 1.) / dmu + 1.) - 1;
+    }
+
+    f = fmu[gin][i_gout][imu];
+  }
+  return f;
+}
+
+//==============================================================================
+
+void
+ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
+{
+  // Sample the outgoing energy using the base-class method
+  int i_gout;
+  sample_energy(gin, gout, i_gout);
+
+  // Determine the outgoing cosine bin
+  double xi = prn();
+
+  int imu;
+  if (xi < dist[gin][i_gout][0]) {
+    imu = 0;
+  } else {
+    // TODO lower_bound?  + 1?
+    imu = std::upper_bound(dist[gin][i_gout].begin(),
+                           dist[gin][i_gout].end(), xi) -
+         dist[gin][i_gout].begin();
+  }
+
+  // Randomly select mu within the imu bin
+  mu = prn() * dmu + this->mu[imu];
+
+  if (mu < -1.) {
+    mu = -1.;
+  } else if (mu > 1.) {
+    mu = 1.;
+  }
+
+  // Update the weight to reflect neutron multiplicity
+  wgt *= mult[gin][i_gout];
+}
+
+//==============================================================================
+
+double_3dvec
+ScattDataHistogram::get_matrix(int max_order)
+{
+  // Get the sizes and initialize the data to 0
+  int groups = energy.size();
+  // We ignore the requested order for Histogram and Tabular representations
+  int order_dim = get_order();
+  double_3dvec matrix = double_3dvec(groups, double_2dvec(groups,
+       double_1dvec(order_dim, 0.)));
+
+  for (int gin = 0; gin < groups; gin++) {
+    for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) {
+      int gout = i_gout + gmin[gin];
+      for (int l = 0; l < order_dim; l++) {
+        matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] *
+             fmu[gin][i_gout][l];
+      }
+    }
+  }
+  return matrix;
+}
+
+//==============================================================================
+
+void
+ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
+                            const double_1dvec& scalars)
+{
+  // Find the max order in the data set and make sure we can combine the sets
+  int max_order = those_scatts[0]->get_order();
+  for (int i = 0; i < those_scatts.size(); i++) {
+    // Lets also make sure these items are combineable
+    ScattDataHistogram* that = dynamic_cast<ScattDataHistogram*>(those_scatts[i]);
+    if (!that) {
+      fatal_error("Cannot combine the ScattData objects!");
+    }
+    if (max_order != that->get_order()) {
+      fatal_error("Cannot combine the ScattData objects!");
+    }
+  }
+
+  int groups = those_scatts[0] -> energy.size();
+
+  int_1dvec in_gmin(groups);
+  int_1dvec in_gmax(groups);
+  double_3dvec sparse_scatter(groups);
+  double_2dvec sparse_mult(groups);
+
+  // The rest of the steps do not depend on the type of angular representation
+  // so we use a base class method to sum up xs and create new energy and mult
+  // matrices
+  ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
+                             sparse_mult, sparse_scatter);
+
+  // Got everything we need, store it.
+  init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
+}
+
+//==============================================================================
+// ScattDataTabular methods
+//==============================================================================
+
+void
+ScattDataTabular::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
+     const double_2dvec& in_mult, const double_3dvec& coeffs)
+{
+  int groups = coeffs.size();
+  int order = coeffs[0][0].size();
+
+  // make a copy of coeffs that we can use to both extract data and normalize
+  double_3dvec matrix = coeffs;
+
+  // Build the angular distribution mu values
+  mu = double_1dvec(order);
+  dmu = 2. / (order - 1);
+  mu[0] = -1.;
+  for (int imu = 1; imu < order - 1; imu++) {
+    mu[imu] = -1. + imu * dmu;
+  }
+  mu[order - 1] = 1.;
+
+  // Get the scattering cross section value by integrating the distribution
+  // over all mu points and then combining over all outgoing groups
+  scattxs.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    scattxs[gin] = 0.;
+    for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) {
+      for (int imu = 1; imu < order; imu++) {
+        scattxs[gin] += 0.5 * dmu * (matrix[gin][i_gout][imu - 1] +
+                                     matrix[gin][i_gout][imu]);
+      }
+    }
+  }
+
+  // Build the energy transfer matrix from data in the variable matrix
+  double_2dvec in_energy(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
+    in_energy[gin].resize(num_groups);
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      double norm = 0.;
+      for (int imu = 1; imu < order; imu++) {
+        norm += 0.5 * dmu * (matrix[gin][i_gout][imu - 1] +
+                             matrix[gin][i_gout][imu]);
+      }
+      in_energy[gin][i_gout] = norm;
+    }
+  }
+
+  // Initialize the base class attributes
+  ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult);
+
+  // Calculate f(mu) and integrate it so we can avoid rejection sampling
+  fmu.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = gmax[gin] - gmin[gin] + 1;
+    fmu[gin].resize(num_groups);
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      fmu[gin][i_gout].resize(order);
+      // The variable matrix contains f(mu); so directly assign it
+      fmu[gin][i_gout] = matrix[gin][i_gout];
+
+      // Ensure positivity
+      for (auto& val : fmu[gin][i_gout]) {
+        if (val < 0.) val = 0.;
+      }
+
+      // Now re-normalize for numerical integration issues and to take care of
+      // the above negative fix-up.  Also accrue the CDF
+      double norm = 0.;
+      for (int imu = 1; imu < order; imu++) {
+        norm += 0.5 * dmu * (fmu[gin][i_gout][imu - 1] +
+                             fmu[gin][i_gout][imu]);
+        // incorporate to the CDF
+        dist[gin][i_gout][imu] = norm;
+      }
+
+      // now do the normalization
+      if (norm > 0.) {
+        for (int imu = 0; imu < order; imu++) {
+          fmu[gin][i_gout][imu] /= norm;
+          dist[gin][i_gout][imu] /= norm;
+        }
+      }
+    }
+  }
+}
+
+//==============================================================================
+
+double
+ScattDataTabular::calc_f(int gin, int gout, double mu)
+{
+  double f;
+  if ((gout < gmin[gin]) || (gout > gmax[gin])) {
+    f = 0.;
+  } else {
+    // Find mu bin
+    int i_gout = gout - gmin[gin];
+    int imu;
+    if (mu == 1.) {
+      // use size -2 to have the index one before the end
+      imu = this->mu.size() - 2;
+    } else {
+      imu = std::floor((mu + 1.) / dmu + 1.) - 1;
+    }
+
+    double r = (mu - this->mu[imu]) / (this->mu[imu + 1] - this->mu[imu]);
+    f = (1. - r) * fmu[gin][i_gout][imu] + r * fmu[gin][i_gout][imu + 1];
+  }
+  return f;
+}
+
+//==============================================================================
+
+void
+ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
+{
+  // Sample the outgoing energy using the base-class method
+  int i_gout;
+  sample_energy(gin, gout, i_gout);
+
+  // Determine the outgoing cosine bin
+  int NP = this->mu.size();
+  double xi = prn();
+
+  double c_k = dist[gin][i_gout][0];
+  int k;
+  for (k = 0; k < NP - 1; k++) {
+    double c_k1 = dist[gin][i_gout][k + 1];
+    if (xi < c_k1) break;
+    c_k = c_k1;
+  }
+
+  // Check to make sure k is <= NP - 1
+  k = std::min(k, NP - 2);
+
+  // Find the pdf values we want
+  double p0 = fmu[gin][i_gout][k];
+  double mu0 = this -> mu[k];
+  double p1 = fmu[gin][i_gout][k + 1];
+  double mu1 = this -> mu[k + 1];
+
+  if (p0 == p1) {
+    mu = mu0 + (xi - c_k) / p0;
+  } else {
+    double frac = (p1 - p0) / (mu1 - mu0);
+    mu = mu0 + (std::sqrt(std::max(0., p0 * p0 + 2. * frac * (xi - c_k)))
+                - p0) / frac;
+  }
+
+  if (mu < -1.) {
+    mu = -1.;
+  } else if (mu > 1.) {
+    mu = 1.;
+  }
+
+  // Update the weight to reflect neutron multiplicity
+  wgt *= mult[gin][i_gout];
+}
+
+//==============================================================================
+
+double_3dvec
+ScattDataTabular::get_matrix(int max_order)
+{
+  // Get the sizes and initialize the data to 0
+  int groups = energy.size();
+  // We ignore the requested order for Histogram and Tabular representations
+  int order_dim = get_order();
+  double_3dvec matrix = double_3dvec(groups, double_2dvec(groups,
+       double_1dvec(order_dim, 0.)));
+
+  for (int gin = 0; gin < groups; gin++) {
+    for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) {
+      int gout = i_gout + gmin[gin];
+      for (int l = 0; l < order_dim; l++) {
+        matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] *
+             fmu[gin][i_gout][l];
+      }
+    }
+  }
+  return matrix;
+}
+
+//==============================================================================
+
+void
+ScattDataTabular::combine(const std::vector<ScattData*>& those_scatts,
+                          const double_1dvec& scalars)
+{
+  // Find the max order in the data set and make sure we can combine the sets
+  int max_order = those_scatts[0]->get_order();
+  for (int i = 0; i < those_scatts.size(); i++) {
+    // Lets also make sure these items are combineable
+    ScattDataTabular* that = dynamic_cast<ScattDataTabular*>(those_scatts[i]);
+    if (!that) {
+      fatal_error("Cannot combine the ScattData objects!");
+    }
+    if (max_order != that->get_order()) {
+      fatal_error("Cannot combine the ScattData objects!");
+    }
+  }
+
+  int groups = those_scatts[0] -> energy.size();
+
+  int_1dvec in_gmin(groups);
+  int_1dvec in_gmax(groups);
+  double_3dvec sparse_scatter(groups);
+  double_2dvec sparse_mult(groups);
+
+  // The rest of the steps do not depend on the type of angular representation
+  // so we use a base class method to sum up xs and create new energy and mult
+  // matrices
+  ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
+                             sparse_mult, sparse_scatter);
+
+  // Got everything we need, store it.
+  init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
+}
+
+//==============================================================================
+// module-level methods
+//==============================================================================
+
+void
+convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
+                            int n_mu)
+{
+  // See if the user wants us to figure out how many points to use
+  if (n_mu == C_NONE) {
+    // then we will use 2 pts if its P0, or the default if a higher order
+    // TODO use an error minimization algorithm that also picks n_mu
+    if (leg.get_order() == 0) {
+      n_mu = 2;
+    } else {
+      n_mu = DEFAULT_NMU;
+    }
+  }
+
+  tab.base_init(n_mu, leg.gmin, leg.gmax, leg.energy, leg.mult);
+  tab.scattxs = leg.scattxs;
+
+  // Build mu and dmu
+  tab.mu = double_1dvec(n_mu);
+  tab.dmu = 2. / (n_mu - 1);
+  tab.mu[0] = -1.;
+  for (int imu = 1; imu < n_mu - 1; imu++) {
+    tab.mu[imu] = -1. + imu  * tab.dmu;
+  }
+  tab.mu[n_mu - 1] = 1.;
+
+  // Calculate f(mu) and integrate it so we can avoid rejection sampling
+  int groups = tab.energy.size();
+  tab.fmu.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = tab.gmax[gin] - tab.gmin[gin] + 1;
+    tab.fmu[gin].resize(num_groups);
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      tab.fmu[gin][i_gout].resize(n_mu);
+      for (int imu = 0; imu < n_mu; imu++) {
+        tab.fmu[gin][i_gout][imu] =
+             evaluate_legendre_c(leg.dist[gin][i_gout].size() - 1,
+                                 leg.dist[gin][i_gout].data(), tab.mu[imu]);
+      }
+
+      // Ensure positivity
+      for (auto& val : tab.fmu[gin][i_gout]) {
+        if (val < 0.) val = 0.;
+      }
+
+      // Now re-normalize for numerical integration issues and to take care of
+      // the above negative fix-up.  Also accrue the CDF
+      double norm = 0.;
+      tab.dist[gin][i_gout][0] = 0.;
+      for (int imu = 1; imu < n_mu; imu++) {
+        norm += 0.5 * tab.dmu * (tab.fmu[gin][i_gout][imu - 1] +
+                                 tab.fmu[gin][i_gout][imu]);
+        // incorporate to the CDF
+        tab.dist[gin][i_gout][imu] = norm;
+      }
+
+      // now do the normalization
+      if (norm > 0.) {
+        for (int imu = 0; imu < n_mu; imu++) {
+          tab.fmu[gin][i_gout][imu] /= norm;
+          tab.dist[gin][i_gout][imu] /= norm;
+        }
+      }
+    }
+  }
+}
+
+} // namespace openmc

src/scattdata.h

@@ -0,0 +1,258 @@
+//! \file scattdata.h
+//! A collection of multi-group scattering data classes
+
+#ifndef SCATTDATA_H
+#define SCATTDATA_H
+
+#include <vector>
+
+namespace openmc {
+
+// forward declarations so we can name our friend functions
+class ScattDataLegendre;
+class ScattDataTabular;
+
+//==============================================================================
+// SCATTDATA contains all the data needed to describe the scattering energy and
+// angular distribution data
+//==============================================================================
+
+class ScattData {
+  protected:
+    //! \brief Initializes the attributes of the base class.
+    void
+    base_init(int order, const int_1dvec& in_gmin, const int_1dvec& in_gmax,
+         const double_2dvec& in_energy, const double_2dvec& in_mult);
+
+    //! \brief Combines microscopic ScattDatas into a macroscopic one.
+    void
+    base_combine(int max_order, const std::vector<ScattData*>& those_scatts,
+         const double_1dvec& scalars, int_1dvec& in_gmin, int_1dvec& in_gmax,
+         double_2dvec& sparse_mult, double_3dvec& sparse_scatter);
+
+  public:
+
+    double_2dvec energy; // Normalized p0 matrix for sampling Eout
+    double_2dvec mult;   // nu-scatter multiplication (nu-scatt/scatt)
+    double_3dvec dist;   // Angular distribution
+    int_1dvec gmin;      // minimum outgoing group
+    int_1dvec gmax;      // maximum outgoing group
+    double_1dvec scattxs; // Isotropic Sigma_{s,g_{in}}
+
+    //! \brief Calculates the value of normalized f(mu).
+    //!
+    //! The value of f(mu) is normalized as in the integral of f(mu)dmu across
+    //! [-1,1] is 1.
+    //!
+    //! @param gin Incoming energy group of interest.
+    //! @param gout Outgoing energy group of interest.
+    //! @param mu Cosine of the change-in-angle of interest.
+    //! @return The value of f(mu).
+    virtual double
+    calc_f(int gin, int gout, double mu) = 0;
+
+    //! \brief Samples the outgoing energy and angle from the ScattData info.
+    //!
+    //! @param gin Incoming energy group.
+    //! @param gout Sampled outgoing energy group.
+    //! @param mu Sampled cosine of the change-in-angle.
+    //! @param wgt Weight of the particle to be adjusted.
+    virtual void
+    sample(int gin, int& gout, double& mu, double& wgt) = 0;
+
+    //! \brief Initializes the ScattData object from a given scatter and
+    //!   multiplicity matrix.
+    //!
+    //! @param in_gmin List of minimum outgoing groups for every incoming group
+    //! @param in_gmax List of maximum outgoing groups for every incoming group
+    //! @param in_mult Input sparse multiplicity matrix
+    //! @param coeffs Input sparse scattering matrix
+    virtual void
+    init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
+         const double_2dvec& in_mult, const double_3dvec& coeffs) = 0;
+
+    //! \brief Combines the microscopic data.
+    //!
+    //! @param those_scatts Microscopic objects to combine.
+    //! @param scalars Scalars to multiply the microscopic data by.
+    virtual void
+    combine(const std::vector<ScattData*>& those_scatts,
+         const double_1dvec& scalars) = 0;
+
+    //! \brief Getter for the dimensionality of the scattering order.
+    //!
+    //! If Legendre this is the "n" in "Pn"; for Tabular, this is the number
+    //! of points, and for Histogram this is the number of bins.
+    //!
+    //! @return The order.
+    virtual int
+    get_order() = 0;
+
+    //! \brief Builds a dense scattering matrix from the constituent parts
+    //!
+    //! @param max_order If Legendre this is the maximum value of "n" in "Pn"
+    //!   requested; ignored otherwise.
+    //! @return The dense scattering matrix.
+    virtual double_3dvec
+    get_matrix(int max_order) = 0;
+
+    //! \brief Samples the outgoing energy from the ScattData info.
+    //!
+    //! @param gin Incoming energy group.
+    //! @param gout Sampled outgoing energy group.
+    //! @param i_gout Sampled outgoing energy group index.
+    void
+    sample_energy(int gin, int& gout, int& i_gout);
+
+    //! \brief Provides a cross section value given certain parameters
+    //!
+    //! @param xstype Type of cross section requested, according to the
+    //!   enumerated constants.
+    //! @param gin Incoming energy group.
+    //! @param gout Outgoing energy group; use nullptr if irrelevant, or if a
+    //!   sum is requested.
+    //! @param mu Cosine of the change-in-angle, for scattering quantities;
+    //!   use nullptr if irrelevant.
+    //! @return Requested cross section value.
+    double
+    get_xs(int xstype, int gin, const int* gout, const double* mu);
+};
+
+//==============================================================================
+// ScattDataLegendre represents the angular distributions as Legendre kernels
+//==============================================================================
+
+class ScattDataLegendre: public ScattData {
+
+  protected:
+
+    // Maximal value for rejection sampling from a rectangle
+    double_2dvec max_val;
+
+    // Friend convert_legendre_to_tabular so it has access to protected
+    // parameters
+    friend void
+    convert_legendre_to_tabular(ScattDataLegendre& leg,
+         ScattDataTabular& tab, int n_mu);
+
+  public:
+
+    void
+    init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
+         const double_2dvec& in_mult, const double_3dvec& coeffs);
+
+    void
+    combine(const std::vector<ScattData*>& those_scatts,
+            const double_1dvec& scalars);
+
+    //! \brief Find the maximal value of the angular distribution to use as a
+    // bounding box with rejection sampling.
+    void
+    update_max_val();
+
+    double
+    calc_f(int gin, int gout, double mu);
+
+    void
+    sample(int gin, int& gout, double& mu, double& wgt);
+
+    int
+    get_order() {return dist[0][0].size() - 1;};
+
+    double_3dvec
+    get_matrix(int max_order);
+};
+
+//==============================================================================
+// ScattDataHistogram represents the angular distributions as a histogram, as it
+// would be if it came from a "mu" tally in OpenMC
+//==============================================================================
+
+class ScattDataHistogram: public ScattData {
+
+  protected:
+
+    double_1dvec mu;  // Angle distribution mu bin boundaries
+    double dmu;       // Quick storage of the spacing between the mu bin points
+    double_3dvec fmu; // The angular distribution histogram
+
+  public:
+
+    void
+    init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
+         const double_2dvec& in_mult, const double_3dvec& coeffs);
+
+    void
+    combine(const std::vector<ScattData*>& those_scatts,
+            const double_1dvec& scalars);
+
+    double
+    calc_f(int gin, int gout, double mu);
+
+    void
+    sample(int gin, int& gout, double& mu, double& wgt);
+
+    int
+    get_order() {return dist[0][0].size();};
+
+    double_3dvec
+    get_matrix(int max_order);
+};
+
+//==============================================================================
+// ScattDataTabular represents the angular distributions as a table of mu and
+// f(mu)
+//==============================================================================
+
+class ScattDataTabular: public ScattData {
+
+  protected:
+
+    double_1dvec mu;  // Angle distribution mu grid points
+    double dmu;       // Quick storage of the spacing between the mu points
+    double_3dvec fmu; // The angular distribution function
+
+    // Friend convert_legendre_to_tabular so it has access to protected
+    // parameters
+    friend void
+    convert_legendre_to_tabular(ScattDataLegendre& leg,
+         ScattDataTabular& tab, int n_mu);
+
+  public:
+
+    void
+    init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
+         const double_2dvec& in_mult, const double_3dvec& coeffs);
+
+    void
+    combine(const std::vector<ScattData*>& those_scatts,
+            const double_1dvec& scalars);
+
+    double
+    calc_f(int gin, int gout, double mu);
+
+    void
+    sample(int gin, int& gout, double& mu, double& wgt);
+
+    int
+    get_order() {return dist[0][0].size();};
+
+    double_3dvec get_matrix(int max_order);
+};
+
+//==============================================================================
+// Function to convert Legendre functions to tabular
+//==============================================================================
+
+//! \brief Converts a ScattDatalegendre to a ScattDataHistogram
+//!
+//! @param leg The initial ScattDataLegendre object.
+//! @param leg The resultant ScattDataTabular object.
+//! @param n_mu The number of mu points to use when building the
+//!   ScattDataTabular object.
+void
+convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
+                            int n_mu);
+
+} // namespace openmc
+#endif // SCATTDATA_H

src/scattdata_header.F90

@@ -1,852 +0,0 @@
-module scattdata_header
-
-  use algorithm,  only: binary_search
-  use constants
-  use error,      only: fatal_error
-  use math
-  use random_lcg, only: prn
-
-  implicit none
-
-
-!===============================================================================
-! JAGGED1D and JAGGED2D is a type which allows for jagged 1-D or 2-D array.
-!===============================================================================
-
-  type :: Jagged2D
-    real(8), allocatable :: data(:, :)
-  end type Jagged2D
-
-  type :: Jagged1D
-    real(8), allocatable :: data(:)
-  end type Jagged1D
-
-!===============================================================================
-! SCATTDATA contains all the data to describe the scattering energy and
-! angular distribution
-!===============================================================================
-
-  type, abstract :: ScattData
-    ! The data attribute of the energy, mult, and dist arrays
-    ! are not necessarily 1-indexed as they instead will be allocated
-    ! from a minimum outgoing group to an outgoing minimum group.
-    ! Normalized p0 matrix on its own for sampling energy
-    type(Jagged1D), allocatable :: energy(:) ! (Gin % data(Gout))
-    ! Nu-scatter multiplication (i.e. nu-scatt/scatt)
-    type(Jagged1D), allocatable :: mult(:)   ! (Gin % data(Gout))
-    ! Angular distribution
-    type(Jagged2D), allocatable :: dist(:)   ! (Gin % data(Order/Nmu, Gout)
-    integer, allocatable :: gmin(:)     ! Minimum outgoing group
-    integer, allocatable :: gmax(:)     ! Maximum outgoing group
-    real(8), allocatable :: scattxs(:)  ! Isotropic Sigma_{s,g_{in}}
-
-  contains
-    procedure(scattdata_init_),   deferred :: init   ! Initializes ScattData
-    procedure(scattdata_calc_f_), deferred :: calc_f ! Calculates f, given mu
-    procedure(scattdata_sample_), deferred :: sample ! sample the scatter event
-    procedure :: get_matrix => scattdata_get_matrix  ! Rebuild scattering matrix
-  end type ScattData
-
-  abstract interface
-    subroutine scattdata_init_(this, gmin, gmax, mult, coeffs)
-      import ScattData, Jagged1D, Jagged2D
-      class(ScattData), intent(inout) :: this      ! Object to work with
-      integer,          intent(in)    :: gmin(:)   ! Min Gout
-      integer,          intent(in)    :: gmax(:)   ! Max Gout
-      type(Jagged1D),   intent(in)    :: mult(:)   ! Scatter Prod'n Matrix
-      type(Jagged2D),   intent(in)    :: coeffs(:) ! Coefficients to use
-    end subroutine scattdata_init_
-
-    pure function scattdata_calc_f_(this, gin, gout, mu) result(f)
-      import ScattData
-      class(ScattData), intent(in) :: this ! Scattering Object to work with
-      integer, intent(in)          :: gin  ! Incoming Energy Group
-      integer, intent(in)          :: gout ! Outgoing Energy Group
-      real(8), intent(in)          :: mu   ! Angle of interest
-      real(8)                      :: f    ! Return value of f(mu)
-
-    end function scattdata_calc_f_
-
-    subroutine scattdata_sample_(this, gin, gout, mu, wgt)
-      import ScattData
-      class(ScattData), intent(in)    :: this ! Scattering Object to work with
-      integer,          intent(in)    :: gin  ! Incoming neutron group
-      integer,          intent(out)   :: gout ! Sampled outgoin group
-      real(8),          intent(out)   :: mu   ! Sampled change in angle
-      real(8),          intent(inout) :: wgt  ! Particle weight
-    end subroutine scattdata_sample_
-  end interface
-
-  type, extends(ScattData) :: ScattDataLegendre
-    ! Maximal value for rejection sampling from rectangle
-    type(Jagged1D), allocatable :: max_val(:) ! (Gin % data(Gout))
-  contains
-    procedure :: init => scattdatalegendre_init
-    procedure :: calc_f          => scattdatalegendre_calc_f
-    procedure :: sample          => scattdatalegendre_sample
-  end type ScattDataLegendre
-
-  type, extends(ScattData)      :: ScattDataHistogram
-    real(8), allocatable        :: mu(:) ! Mu bins
-    real(8)                     :: dmu   ! Mu spacing
-    ! Histogram of f(mu) (dist has CDF)
-    type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
-  contains
-    procedure :: init => scattdatahistogram_init
-    procedure :: calc_f          => scattdatahistogram_calc_f
-    procedure :: sample          => scattdatahistogram_sample
-    procedure :: get_matrix      => scattdatahistogram_get_matrix
-  end type ScattDataHistogram
-
-  type, extends(ScattData)      :: ScattDataTabular
-    real(8), allocatable        :: mu(:)       ! Mu bins
-    real(8)                     :: dmu         ! Mu spacing
-    ! PDF of f(mu) (dist has CDF)
-    type(Jagged2D), allocatable :: fmu(:) ! (Gin % data(Order/Nmu x Gout)
-  contains
-    procedure :: init => scattdatatabular_init
-    procedure :: calc_f          => scattdatatabular_calc_f
-    procedure :: sample          => scattdatatabular_sample
-    procedure :: get_matrix      => scattdatatabular_get_matrix
-  end type ScattDataTabular
-
-!===============================================================================
-! SCATTDATACONTAINER allocatable array for storing ScattData Objects (for angle)
-!===============================================================================
-
-  type ScattDataContainer
-    class(ScattData), allocatable :: obj
-  end type ScattDataContainer
-
-contains
-
-!===============================================================================
-! SCATTDATA*_INIT builds the scattdata object
-!===============================================================================
-
-    subroutine scattdata_init(this, order, gmin, gmax, energy, mult)
-      class(ScattData), intent(inout) :: this      ! Object to work on
-      integer,          intent(in)    :: order     ! Data Order
-      integer,          intent(in)    :: gmin(:)   ! Min Gout
-      integer,          intent(in)    :: gmax(:)   ! Max Gout
-      type(Jagged1D),   intent(inout) :: energy(:) ! Energy Transfer Matrix
-      type(Jagged1D),   intent(in)    :: mult(:)   ! Scatter Prod'n Matrix
-
-      integer :: groups, gin
-      real(8) :: norm
-
-      groups = size(energy, dim=1)
-
-      allocate(this % gmin(groups))
-      allocate(this % gmax(groups))
-      allocate(this % energy(groups))
-      allocate(this % mult(groups))
-      allocate(this % dist(groups))
-
-      this % gmin = gmin
-      this % gmax = gmax
-
-      ! Set the outgoing energy PDF values
-      do gin = 1, groups
-        ! Make sure energy is normalized (i.e., CDF is 1)
-        norm = sum(energy(gin) % data(:))
-        if (norm /= ZERO) energy(gin) % data(:) = energy(gin) % data(:) / norm
-        ! Set the values
-        allocate(this % energy(gin) % data(gmin(gin):gmax(gin)))
-        this % energy(gin) % data(:) = energy(gin) % data(:)
-        allocate(this % mult(gin) % data(gmin(gin):gmax(gin)))
-        this % mult(gin) % data(gmin(gin):gmax(gin)) = &
-             mult(gin) % data(gmin(gin):gmax(gin))
-        allocate(this % dist(gin) % data(order, gmin(gin):gmax(gin)))
-        this % dist(gin) % data = ZERO
-      end do
-    end subroutine scattdata_init
-
-    subroutine scattdatalegendre_init(this, gmin, gmax, mult, coeffs)
-      class(ScattDataLegendre), intent(inout) :: this      ! Object to work on
-      integer,                  intent(in)    :: gmin(:)   ! Min Gout
-      integer,                  intent(in)    :: gmax(:)   ! Max Gout
-      type(Jagged1D),           intent(in)    :: mult(:)   ! Scatter Prod'n Matrix
-      type(Jagged2D),           intent(in)    :: coeffs(:) ! Coefficients to use
-
-      real(8) :: dmu, mu, f, norm
-      integer :: imu, Nmu, gout, gin, groups, order
-      type(Jagged1D), allocatable :: energy(:)
-      type(Jagged2D), allocatable :: matrix(:)
-
-      groups = size(coeffs)
-      order = size(coeffs(1) % data, dim=1)
-
-      ! make a copy of coeffs that we can use to extract data and normalize
-      allocate(matrix(groups))
-      do gin = 1, groups
-        allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
-        matrix(gin) % data = coeffs(gin) % data
-      end do
-
-      ! Get scattxs value
-      allocate(this % scattxs(groups))
-      ! Get this by summing the un-normalized P0 coefficient in matrix
-      ! over all outgoing groups
-      do gin = 1, groups
-        this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
-      end do
-
-      allocate(energy(groups))
-      ! Build energy transfer probability matrix from data in matrix
-      ! while also normalizing matrix itself (making CDF of f(mu=1)=1)
-      do gin = 1, groups
-        allocate(energy(gin) % data(gmin(gin):gmax(gin)))
-        energy(gin) % data = ZERO
-        do gout = gmin(gin), gmax(gin)
-          norm = matrix(gin) % data(1, gout)
-          energy(gin) % data(gout) = norm
-          if (norm /= ZERO) then
-            matrix(gin) % data(:, gout) = matrix(gin) % data(:, gout) / norm
-          end if
-        end do
-      end do
-
-      call scattdata_init(this, order, gmin, gmax, energy, mult)
-
-      allocate(this % max_val(groups))
-      ! Set dist values from matrix and initialize max_val
-      do gin = 1, groups
-        do gout = gmin(gin), gmax(gin)
-          this % dist(gin) % data(:, gout) = matrix(gin) % data(:, gout)
-        end do
-        allocate(this % max_val(gin) % data(gmin(gin):gmax(gin)))
-        this % max_val(gin) % data(:) = ZERO
-      end do
-
-      ! Step through the polynomial with fixed number of points to identify
-      ! the maximal value.
-      Nmu = 1001
-      dmu = TWO / real(Nmu - 1, 8)
-      do gin = 1, groups
-        do gout = gmin(gin), gmax(gin)
-          do imu = 1, Nmu
-            ! Update mu. Do first and last seperate to avoid float errors
-            if (imu == 1) then
-              mu = -ONE
-            else if (imu == Nmu) then
-              mu = ONE
-            else
-              mu = -ONE + real(imu - 1, 8) * dmu
-            end if
-            ! Calculate probability
-            f = this % calc_f(gin,gout,mu)
-            ! If this is a new max, store it.
-            if (f > this % max_val(gin) % data(gout)) &
-                 this % max_val(gin) % data(gout) = f
-          end do
-          ! Finally, since we may not have caught the exact max, add 10% margin
-          this % max_val(gin) % data(gout) = &
-               this % max_val(gin) % data(gout) * 1.1_8
-        end do
-      end do
-    end subroutine scattdatalegendre_init
-
-    subroutine scattdatahistogram_init(this, gmin, gmax, mult, coeffs)
-      class(ScattDataHistogram), intent(inout) :: this      ! Object to work on
-      integer,                   intent(in)    :: gmin(:)   ! Min Gout
-      integer,                   intent(in)    :: gmax(:)   ! Max Gout
-      type(Jagged1D),            intent(in)    :: mult(:)   ! Scatter Prod'n Matrix
-      type(Jagged2D),            intent(in)    :: coeffs(:) ! Coefficients to use
-
-      integer :: imu, gin, gout, groups, order
-      real(8) :: norm
-      type(Jagged1D), allocatable :: energy(:)
-      type(Jagged2D), allocatable :: matrix(:)
-
-      groups = size(coeffs)
-      order = size(coeffs(1) % data, dim=1)
-
-      ! make a copy of coeffs that we can use to extract data and normalize
-      allocate(matrix(groups))
-      do gin = 1, groups
-        allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
-        matrix(gin) % data(:, :) = coeffs(gin) % data(:, :)
-      end do
-
-      ! Get scattxs value
-      allocate(this % scattxs(groups))
-      ! Get this by summing the un-normalized angular distribution in matrix
-      ! over all outgoing groups
-      do gin = 1, groups
-        this % scattxs(gin) = sum(matrix(gin) % data(:, :))
-      end do
-
-      allocate(energy(groups))
-      ! Build energy transfer probability matrix from data in matrix
-      ! while also normalizing matrix itself (making CDF of f(mu=1)=1)
-      do gin = 1, groups
-        allocate(energy(gin) % data(gmin(gin):gmax(gin)))
-        do gout = gmin(gin), gmax(gin)
-          norm = sum(matrix(gin) % data(:, gout))
-          energy(gin) % data(gout) = norm
-          if (norm /= ZERO) then
-            matrix(gin) % data(:, gout) = matrix(gin) % data(:, gout) / norm
-          end if
-        end do
-      end do
-
-      call scattdata_init(this, order, gmin, gmax, energy, mult)
-
-      allocate(this % mu(order))
-      this % dmu = TWO / real(order, 8)
-      this % mu(1) = -ONE
-      do imu = 2, order
-        this % mu(imu) = -ONE + real(imu - 1, 8) * this % dmu
-      end do
-
-      ! Integrate this histogram so we can avoid rejection sampling while
-      ! also saving the original histogram in fmu
-      allocate(this % fmu(groups))
-      do gin = 1, groups
-        allocate(this % fmu(gin) % data(order, gmin(gin):gmax(gin)))
-        do gout = gmin(gin), gmax(gin)
-          ! Store the histogram
-          this % fmu(gin) % data(:, gout) = matrix(gin) % data(:, gout)
-          ! Integrate the histogram
-          this % dist(gin) % data(1, gout) = &
-               this % dmu * matrix(gin) % data(1, gout)
-          do imu = 2, order
-            this % dist(gin) % data(imu, gout) = &
-                 this % dmu * matrix(gin) % data(imu, gout) + &
-                 this % dist(gin) % data(imu - 1, gout)
-          end do
-
-          ! Normalize the integral to unity
-          norm = this % dist(gin) % data(order, gout)
-          if (norm > ZERO) then
-            this % fmu(gin) % data(:, gout) = &
-                 this % fmu(gin) % data(:, gout) / norm
-            this % dist(gin) % data(:, gout) = &
-                 this % dist(gin) % data(:, gout) / norm
-          end if
-        end do
-      end do
-
-    end subroutine scattdatahistogram_init
-
-    subroutine scattdatatabular_init(this, gmin, gmax, mult, coeffs)
-      class(ScattDataTabular), intent(inout) :: this      ! Object to work on
-      integer,                 intent(in)    :: gmin(:)   ! Min Gout
-      integer,                 intent(in)    :: gmax(:)   ! Max Gout
-      type(Jagged1D),          intent(in)    :: mult(:)   ! Scatter Prod'n Matrix
-      type(Jagged2D),          intent(in)    :: coeffs(:) ! Coefficients to use
-
-      integer :: imu, gin, gout, groups, order
-      real(8) :: norm
-      type(Jagged1D), allocatable :: energy(:)
-      type(Jagged2D), allocatable :: matrix(:)
-
-      groups = size(coeffs)
-      order = size(coeffs(1) % data, dim=1)
-
-      ! make a copy of coeffs that we can use to extract data and normalize
-      allocate(matrix(groups))
-      do gin = 1, groups
-        allocate(matrix(gin) % data(order, gmin(gin):gmax(gin)))
-        matrix(gin) % data = coeffs(gin) % data
-      end do
-
-      ! Build the angular distribution mu values
-      allocate(this % mu(order))
-      this % dmu = TWO / real(order - 1, 8)
-      this % mu(1) = -ONE
-      do imu = 2, order - 1
-        this % mu(imu) = -ONE + real(imu - 1, 8) * this % dmu
-      end do
-      this % mu(order) = ONE
-
-      ! Get scattxs
-      allocate(this % scattxs(groups))
-      ! Get this by integrating the scattering distribution over all mu points
-      ! and then combining over all outgoing groups
-      ! over all outgoing groups
-      do gin = 1, groups
-        norm = ZERO
-        do gout = gmin(gin), gmax(gin)
-          do imu = 2, order
-            norm = norm + HALF * this % dmu * &
-                 (matrix(gin) % data(imu - 1, gout) + &
-                  matrix(gin) % data(imu, gout))
-          end do
-        end do
-        this % scattxs(gin) = norm
-      end do
-
-      allocate(energy(groups))
-      ! Build energy transfer probability matrix from data in matrix
-      do gin = 1, groups
-        allocate(energy(gin) % data(gmin(gin):gmax(gin)))
-        do gout = gmin(gin), gmax(gin)
-          norm = ZERO
-          do imu = 2, order
-            norm = norm + HALF * this % dmu * &
-                 (matrix(gin) % data(imu - 1, gout) + &
-                  matrix(gin) % data(imu, gout))
-          end do
-          energy(gin) % data(gout) = norm
-        end do
-      end do
-      call scattdata_init(this, order, gmin, gmax, energy, mult)
-
-      ! Calculate f(mu) and integrate it so we can avoid rejection sampling
-      allocate(this % fmu(groups))
-      do gin = 1, groups
-        allocate(this % fmu(gin) % data(order, gmin(gin):gmax(gin)))
-        do gout = gmin(gin), gmax(gin)
-          ! Coeffs contain f(mu), put in f(mu) as that is where the
-          ! PDF lives
-          this % fmu(gin) % data(:, gout) = matrix(gin) % data(:, gout)
-
-          ! Force positivity
-          do imu = 1, order
-            if (this % fmu(gin) % data(imu, gout) < ZERO) then
-              this % fmu(gin) % data(imu, gout) = ZERO
-            end if
-          end do
-
-          ! Re-normalize fmu for numerical integration issues and in case
-          ! the negative fix-up introduced un-normalized data while
-          ! accruing the CDF
-          norm = ZERO
-          do imu = 2, order
-            norm = norm + HALF * this % dmu * &
-                 (this % fmu(gin) % data(imu - 1, gout) + &
-                  this % fmu(gin) % data(imu, gout))
-            this % dist(gin) % data(imu, gout) = norm
-          end do
-          if (norm > ZERO) then
-            this % fmu(gin) % data(:, gout) = &
-                 this % fmu(gin) % data(:, gout) / norm
-            this % dist(gin) % data(:, gout) = &
-                 this % dist(gin) % data(:, gout) / norm
-          end if
-        end do
-      end do
-    end subroutine scattdatatabular_init
-
-!===============================================================================
-! SCATTDATA_*_CALC_F Calculates the value of f given mu (and gin,gout pair)
-!===============================================================================
-
-    pure function scattdatalegendre_calc_f(this, gin, gout, mu) result(f)
-      class(ScattDataLegendre), intent(in) :: this ! The ScattData to evaluate
-      integer, intent(in)                  :: gin  ! Incoming Energy Group
-      integer, intent(in)                  :: gout ! Outgoing Energy Group
-      real(8), intent(in)                  :: mu   ! Angle of interest
-      real(8)                              :: f    ! Return value of f(mu)
-
-      ! Plug mu in to the legendre expansion and go from there
-      if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then
-        f = ZERO
-      else
-        f = evaluate_legendre(this % dist(gin) % data(:, gout), mu)
-      end if
-
-    end function scattdatalegendre_calc_f
-
-    pure function scattdatahistogram_calc_f(this, gin, gout, mu) result(f)
-      class(ScattDataHistogram), intent(in) :: this ! The ScattData to evaluate
-      integer, intent(in)                   :: gin  ! Incoming Energy Group
-      integer, intent(in)                   :: gout ! Outgoing Energy Group
-      real(8), intent(in)                   :: mu   ! Angle of interest
-      real(8)                               :: f    ! Return value of f(mu)
-
-      integer :: imu
-
-      if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then
-        f = ZERO
-      else
-        ! Find mu bin
-        if (mu == ONE) then
-          imu = size(this % fmu(gin) % data, dim=1)
-        else
-          imu = floor((mu + ONE) / this % dmu + ONE)
-        end if
-
-        f = this % fmu(gin) % data(imu, gout)
-      end if
-
-    end function scattdatahistogram_calc_f
-
-    pure function scattdatatabular_calc_f(this, gin, gout, mu) result(f)
-      class(ScattDataTabular), intent(in) :: this ! The ScattData to evaluate
-      integer, intent(in)                 :: gin  ! Incoming Energy Group
-      integer, intent(in)                 :: gout ! Outgoing Energy Group
-      real(8), intent(in)                 :: mu   ! Angle of interest
-      real(8)                             :: f    ! Return value of f(mu)
-
-      integer :: imu
-      real(8) :: r
-
-      if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then
-        f = ZERO
-      else
-        ! Find mu bin
-        if (mu == ONE) then
-          imu = size(this % fmu(gin) % data, dim=1) - 1
-        else
-          imu = floor((mu + ONE) / this % dmu + ONE)
-        end if
-
-        ! Now interpolate to find f(mu)
-        r = (mu - this % mu(imu)) / (this % mu(imu + 1) - this % mu(imu))
-        f = (ONE - r) * this % fmu(gin) % data(imu, gout) + &
-             r * this % fmu(gin) % data(imu + 1, gout)
-      end if
-
-    end function scattdatatabular_calc_f
-
-!===============================================================================
-! SCATTDATA*_SCATTER Samples the outgoing energy and change in angle.
-!===============================================================================
-
-  subroutine scattdatalegendre_sample(this, gin, gout, mu, wgt)
-    class(ScattDataLegendre), intent(in)    :: this ! Scattering object to use
-    integer,                  intent(in)    :: gin  ! Incoming neutron group
-    integer,                  intent(out)   :: gout ! Sampled outgoin group
-    real(8),                  intent(out)   :: mu   ! Sampled change in angle
-    real(8),                  intent(inout) :: wgt  ! Particle weight
-
-    real(8) :: xi     ! Our random number
-    real(8) :: prob   ! Running probability
-    real(8) :: u, f, M
-    integer :: samples
-
-    xi = prn()
-    gout = this % gmin(gin)
-    prob = this % energy(gin) % data(gout)
-
-    do while ((prob < xi) .and. (gout < this % gmax(gin)))
-      gout = gout + 1
-      prob = prob + this % energy(gin) % data(gout)
-    end do
-
-    ! Now we can sample mu using the legendre representation of the scattering
-    ! kernel in data(1:this % order)
-
-    ! Do with rejection sampling from a rectangular bounding box
-    ! Set maximal value
-    M = this % max_val(gin) % data(gout)
-    samples = 0
-    do
-      mu = TWO * prn() - ONE
-      f = this % calc_f(gin, gout, mu)
-      if (f > ZERO) then
-        u = prn() * M
-        if (u <= f) then
-          exit
-        end if
-      end if
-      samples = samples + 1
-      if (samples > MAX_SAMPLE) then
-        call fatal_error("Maximum number of Legendre expansion samples reached!")
-      end if
-    end do
-
-    wgt = wgt * this % mult(gin) % data(gout)
-
-  end subroutine scattdatalegendre_sample
-
-  subroutine scattdatahistogram_sample(this, gin, gout, mu, wgt)
-    class(ScattDataHistogram), intent(in)    :: this ! Scattering object to use
-    integer,                   intent(in)    :: gin  ! Incoming neutron group
-    integer,                   intent(out)   :: gout ! Sampled outgoin group
-    real(8),                   intent(out)   :: mu   ! Sampled change in angle
-    real(8),                   intent(inout) :: wgt  ! Particle weight
-
-    real(8) :: xi     ! Our random number
-    real(8) :: prob   ! Running probability
-    integer :: imu
-
-    xi = prn()
-    gout = this % gmin(gin)
-    prob = this % energy(gin) % data(gout)
-
-    do while ((prob < xi) .and. (gout < this % gmax(gin)))
-      gout = gout + 1
-      prob = prob + this % energy(gin) % data(gout)
-    end do
-
-    xi = prn()
-    if (xi < this % dist(gin) % data(1, gout)) then
-      imu = 1
-    else
-      imu = binary_search(this % dist(gin) % data(:, gout), &
-                          size(this % dist(gin) % data(:, gout)), xi) + 1
-    end if
-
-    ! Randomly select a mu in this bin.
-    mu = prn() * this % dmu + this % mu(imu)
-
-    wgt = wgt * this % mult(gin) % data(gout)
-
-  end subroutine scattdatahistogram_sample
-
-  subroutine scattdatatabular_sample(this, gin, gout, mu, wgt)
-    class(ScattDataTabular), intent(in)    :: this ! Scattering object to use
-    integer,                  intent(in)    :: gin  ! Incoming neutron group
-    integer,                  intent(out)   :: gout ! Sampled outgoin group
-    real(8),                  intent(out)   :: mu   ! Sampled change in angle
-    real(8),                  intent(inout) :: wgt  ! Particle weight
-
-    real(8) :: xi     ! Our random number
-    real(8) :: prob   ! Running probability
-    real(8) :: mu0, frac, mu1
-    real(8) :: c_k, c_k1, p0, p1
-    integer :: k, NP
-
-    xi = prn()
-    gout = this % gmin(gin)
-    prob = this % energy(gin) % data(gout)
-
-    do while ((prob < xi) .and. (gout < this % gmax(gin)))
-      gout = gout + 1
-      prob = prob + this % energy(gin) % data(gout)
-    end do
-
-    ! determine outgoing cosine bin
-    NP = size(this % dist(gin) % data(:, gout))
-    xi = prn()
-
-    c_k = this % dist(gin) % data(1, gout)
-    do k = 1, NP - 1
-      c_k1 = this % dist(gin) % data(k + 1, gout)
-      if (xi < c_k1) exit
-      c_k = c_k1
-    end do
-
-    ! check to make sure k is <= NP - 1
-    k = min(k, NP - 1)
-
-    p0  = this % fmu(gin) % data(k, gout)
-    mu0 = this % mu(k)
-    ! Linear-linear interpolation to find mu value w/in bin.
-    p1  = this % fmu(gin) % data(k + 1, gout)
-    mu1 = this % mu(k + 1)
-
-    if (p0 == p1) then
-      mu = mu0 + (xi - c_k) / p0
-    else
-      frac = (p1 - p0) / (mu1 - mu0)
-      mu = mu0 + &
-           (sqrt(max(ZERO, p0 * p0 + TWO * frac * (xi - c_k))) - p0) / frac
-    end if
-
-    if (mu <= -ONE) then
-      mu = -ONE
-    else if (mu >= ONE) then
-      mu = ONE
-    end if
-
-    wgt = wgt * this % mult(gin) % data(gout)
-
-  end subroutine scattdatatabular_sample
-
-!===============================================================================
-! SCATTDATA*_GET_MATRIX Reproduces the original scattering matrix (densely)
-! using ScattData's information of fmu/dist, energy, and scattxs
-!===============================================================================
-
-    subroutine scattdata_get_matrix(this, req_order, matrix)
-      class(ScattData), intent(in) :: this      ! Scattering Object to work with
-      integer, intent(in)          :: req_order ! Requested order of matrix
-      type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
-
-      integer :: order, groups, gin, gout
-
-      groups = size(this % energy)
-      ! Set gin and gout for getting the order
-      order = min(req_order, size(this % dist(1) % data, dim=1))
-
-      if (allocated(matrix)) deallocate(matrix)
-      allocate(matrix(groups))
-      ! Initialize to 0; this way the zero entries in the dense matrix dont
-      ! need to be explicitly set, requiring a significant increase in the
-      ! lines of code.
-      do gin = 1, groups
-        allocate(matrix(gin) % data(order, groups))
-        do gout = this % gmin(gin), this % gmax(gin)
-          matrix(gin) % data(:, gout) = this % scattxs(gin) * &
-               this % energy(gin) % data(gout) * &
-               this % dist(gin) % data(1:order, gout)
-        end do
-      end do
-    end subroutine scattdata_get_matrix
-
-    subroutine scattdatahistogram_get_matrix(this, req_order, matrix)
-      class(ScattDataHistogram), intent(in) :: this      ! Scattering Object to work with
-      integer, intent(in)                   :: req_order ! Requested order of matrix
-      type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
-
-      integer :: order, groups, gin, gout
-
-      groups = size(this % energy)
-      order = min(req_order, size(this % dist(1) % data, dim=1))
-
-      if (allocated(matrix)) deallocate(matrix)
-      allocate(matrix(groups))
-      ! Initialize to 0; this way the zero entries in the dense matrix dont
-      ! need to be explicitly set, requiring a significant increase in the
-      ! lines of code.
-      do gin = 1, groups
-        allocate(matrix(gin) % data(order, groups))
-        do gout = this % gmin(gin), this % gmax(gin)
-          matrix(gin) % data(:, gout) = this % scattxs(gin) * &
-               this % energy(gin) % data(gout) * &
-               this % fmu(gin) % data(1:order, gout)
-        end do
-      end do
-    end subroutine scattdatahistogram_get_matrix
-
-    subroutine scattdatatabular_get_matrix(this, req_order, matrix)
-      class(ScattDataTabular), intent(in) :: this      ! Scattering Object to work with
-      integer, intent(in)                 :: req_order ! Requested order of matrix
-      type(Jagged2D), allocatable, intent(inout) :: matrix(:) ! Resultant matrix just built
-
-      integer :: order, groups, gin, gout
-
-      groups = size(this % energy)
-      order = min(req_order, size(this % dist(1) % data, dim=1))
-
-      if (allocated(matrix)) deallocate(matrix)
-      allocate(matrix(groups))
-      ! Initialize to 0; this way the zero entries in the dense matrix dont
-      ! need to be explicitly set, requiring a significant increase in the
-      ! lines of code.
-      do gin = 1, groups
-        allocate(matrix(gin) % data(order, groups))
-        do gout = this % gmin(gin), this % gmax(gin)
-          matrix(gin) % data(:, gout) = this % scattxs(gin) * &
-               this % energy(gin) % data(gout) * &
-               this % fmu(gin) % data(1:order, gout)
-        end do
-      end do
-    end subroutine scattdatatabular_get_matrix
-
-!===============================================================================
-! JAGGED_FROM_DENSE_*D Creates a jagged array from a sparse dense matrix.
-! The user can supply a key which indicates the values to remove, but the
-! default is ZERO
-!===============================================================================
-
-    subroutine jagged_from_dense_1D(dense, jagged, lo_bounds_, hi_bounds_, key_)
-      real(8), intent(in)                        :: dense(:, :)
-      type(Jagged1D), allocatable, intent(inout) :: jagged(:)
-      real(8), intent(in), optional              :: key_
-      integer, intent(inout), allocatable, optional :: lo_bounds_(:)
-      integer, intent(inout), allocatable, optional :: hi_bounds_(:)
-
-      real(8) :: key
-      integer :: i, jmin, jmax
-      integer, allocatable :: lo_bounds(:), hi_bounds(:)
-
-      if (present(key_)) then
-        key = key_
-      else
-        key = ZERO
-      end if
-
-      allocate(lo_bounds(size(dense, dim=2)))
-      allocate(hi_bounds(size(dense, dim=2)))
-
-      if (allocated(jagged)) deallocate(jagged)
-      allocate(jagged(size(dense, dim=2)))
-      do i = 1, size(dense, dim=2)
-        ! Find the min and max j values
-        do jmin = 1, size(dense, dim=1)
-          if (dense(jmin, i) /= key) exit
-        end do
-        do jmax = size(dense, dim=1), 1, -1
-          if (dense(jmax, i) /= key) exit
-        end do
-        ! Treat the case of all values matching the key
-        if (jmin > jmax) then
-          jmin = i
-          jmax = i
-        end if
-
-        ! Now store the jagged row
-        allocate(jagged(i) % data(jmin:jmax))
-        jagged(i) % data(jmin:jmax) = dense(jmin:jmax, i)
-
-        lo_bounds(i) = jmin
-        hi_bounds(i) = jmax
-      end do
-
-      if (present(lo_bounds_)) then
-        if (allocated(lo_bounds_)) deallocate(lo_bounds_)
-        allocate(lo_bounds_(size(dense, dim=2)))
-        lo_bounds_ = lo_bounds
-      end if
-      if (present(hi_bounds_)) then
-        if (allocated(hi_bounds_)) deallocate(hi_bounds_)
-        allocate(hi_bounds_(size(dense, dim=2)))
-        hi_bounds_ = hi_bounds
-      end if
-
-    end subroutine jagged_from_dense_1D
-
-    subroutine jagged_from_dense_2D(dense, jagged, lo_bounds_, hi_bounds_, key_)
-      real(8), intent(in)                        :: dense(:, :, :)
-      type(Jagged2D), allocatable, intent(inout) :: jagged(:)
-      real(8), intent(in), optional              :: key_
-      integer, intent(inout), allocatable, optional :: lo_bounds_(:)
-      integer, intent(inout), allocatable, optional :: hi_bounds_(:)
-
-      real(8) :: key
-      integer :: i, jmin, jmax
-      integer, allocatable :: lo_bounds(:), hi_bounds(:)
-
-      if (present(key_)) then
-        key = key_
-      else
-        key = ZERO
-      end if
-
-      allocate(lo_bounds(size(dense, dim=3)))
-      allocate(hi_bounds(size(dense, dim=3)))
-
-      if (allocated(jagged)) deallocate(jagged)
-      allocate(jagged(size(dense, dim=3)))
-      do i = 1, size(dense, dim=3)
-        ! Find the min and max j values
-        do jmin = 1, size(dense, dim=2)
-          if (any(dense(:, jmin, i) /= key)) exit
-        end do
-        do jmax = size(dense, dim=2), 1, -1
-          if (any(dense(:, jmax, i) /= key)) exit
-        end do
-        ! Treat the case of all values matching the key
-        if (jmin > jmax) then
-          jmin = i
-          jmax = i
-        end if
-
-        ! Now store the jagged row
-        allocate(jagged(i) % data(size(dense, dim=1), jmin:jmax))
-        jagged(i) % data(:, jmin:jmax) = dense(:, jmin:jmax, i)
-
-        lo_bounds(i) = jmin
-        hi_bounds(i) = jmax
-      end do
-
-      if (present(lo_bounds_)) then
-        if (allocated(lo_bounds_)) deallocate(lo_bounds_)
-        allocate(lo_bounds_(size(dense, dim=3)))
-        lo_bounds_ = lo_bounds
-      end if
-      if (present(hi_bounds_)) then
-        if (allocated(hi_bounds_)) deallocate(hi_bounds_)
-        allocate(hi_bounds_(size(dense, dim=3)))
-        hi_bounds_ = hi_bounds
-      end if
-
-    end subroutine jagged_from_dense_2D
-
-end module scattdata_header

src/settings.F90

@@ -18,9 +18,9 @@ module settings
   logical :: urr_ptables_on = .true.
 
   ! Default temperature and method for choosing temperatures
-  integer :: temperature_method = TEMPERATURE_NEAREST
+  integer(C_INT) :: temperature_method = TEMPERATURE_NEAREST
   logical :: temperature_multipole = .false.
-  real(8) :: temperature_tolerance = 10.0_8
+  real(C_DOUBLE) :: temperature_tolerance = 10.0_8
   real(8) :: temperature_default = 293.6_8
   real(8) :: temperature_range(2) = [ZERO, ZERO]
 
@@ -30,13 +30,16 @@ module settings
   ! MULTI-GROUP CROSS SECTION RELATED VARIABLES
 
   ! Maximum Data Order
-  integer :: max_order
+  integer(C_INT) :: max_order
 
   ! Whether or not to convert Legendres to tabulars
   logical :: legendre_to_tabular = .true.
 
   ! Number of points to use in the Legendre to tabular conversion
-  integer :: legendre_to_tabular_points = 33
+  integer(C_INT) :: legendre_to_tabular_points = C_NONE
+
+  ! ============================================================================
+  ! SIMULATION VARIABLES
 
   ! Assume all tallies are spatially distinct
   logical :: assume_separate = .false.
@@ -44,9 +47,6 @@ module settings
   ! Use confidence intervals for results instead of standard deviations
   logical :: confidence_intervals = .false.
 
-  ! ============================================================================
-  ! SIMULATION VARIABLES
-
   integer(C_INT64_T), bind(C) :: n_particles = 0   ! # of particles per generation
   integer(C_INT32_T), bind(C) :: n_batches         ! # of batches
   integer(C_INT32_T), bind(C) :: n_inactive        ! # of inactive batches

src/simulation.F90

@@ -20,7 +20,7 @@ module simulation
   use geometry_header, only: n_cells
   use material_header, only: n_materials, materials
   use message_passing
-  use mgxs_header,     only: energy_bins, energy_bin_avg
+  use mgxs_interface,  only: energy_bins, energy_bin_avg
   use nuclide_header,  only: micro_xs, n_nuclides
   use output,          only: header, print_columns, &
                              print_batch_keff, print_generation, print_runtime, &

src/source.F90

@@ -14,7 +14,7 @@ module source
   use hdf5_interface
   use math
   use message_passing,  only: rank
-  use mgxs_header,      only: rev_energy_bins, num_energy_groups
+  use mgxs_interface,   only: rev_energy_bins, num_energy_groups
   use output,           only: write_message
   use particle_header,  only: Particle
   use random_lcg,       only: prn, set_particle_seed, prn_set_stream

src/state_point.F90

@@ -22,7 +22,7 @@ module state_point
   use hdf5_interface
   use mesh_header,        only: RegularMesh, meshes, n_meshes
   use message_passing
-  use mgxs_header,        only: nuclides_MG
+  use mgxs_interface
   use nuclide_header,     only: nuclides
   use output,             only: time_stamp
   use random_lcg,         only: openmc_get_seed, openmc_set_seed
@@ -73,6 +73,7 @@ contains
     character(MAX_WORD_LEN), allocatable :: str_array(:)
     character(C_CHAR), pointer :: string(:)
     character(len=:, kind=C_CHAR), allocatable :: filename_
+    character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name
 
     err = 0
     if (present(filename)) then
@@ -307,11 +308,13 @@ contains
                   str_array(j) = nuclides(tally % nuclide_bins(j)) % name
                 end if
               else
-                i_xs = index(nuclides_MG(tally % nuclide_bins(j)) % obj % name, '.')
+                call get_name_c(tally % nuclide_bins(j), len(temp_name), &
+                                temp_name)
+                i_xs = index(temp_name, '.')
                 if (i_xs > 0) then
-                  str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name(1 : i_xs-1)
+                  str_array(j) = trim(temp_name(1 : i_xs-1))
                 else
-                  str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name
+                  str_array(j) = trim(temp_name)
                 end if
               end if
             else

src/string_functions.cpp

@@ -0,0 +1,30 @@
+#include "string_functions.h"
+
+namespace openmc {
+
+std::string& strtrim(std::string& s)
+{
+  const char* t = " \t\n\r\f\v";
+  s.erase(s.find_last_not_of(t) + 1);
+  s.erase(0, s.find_first_not_of(t));
+  return s;
+}
+
+
+char* strtrim(char* c_str)
+{
+  std::string std_str;
+  std_str.assign(c_str);
+  strtrim(std_str);
+  int length = std_str.copy(c_str, std_str.size());
+  c_str[length] = '\0';
+  return c_str;
+}
+
+
+void to_lower(std::string& str)
+{
+  for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
+}
+
+} // namespace openmc

src/string_functions.h

@@ -0,0 +1,18 @@
+//! \file string_functions.h
+//! A collection of helper routines for C-strings and STL strings
+
+#ifndef STRING_FUNCTIONS_H
+#define STRING_FUNCTIONS_H
+
+#include <string>
+
+namespace openmc {
+
+std::string& strtrim(std::string& s);
+
+char* strtrim(char* c_str);
+
+void to_lower(std::string& str);
+
+} // namespace openmc
+#endif // STRING_FUNCTIONS_H

src/summary.F90

@@ -8,7 +8,7 @@ module summary
   use material_header, only: Material, n_materials
   use mesh_header,     only: RegularMesh
   use message_passing
-  use mgxs_header,     only: nuclides_MG
+  use mgxs_interface
   use nuclide_header
   use output,          only: time_stamp
   use settings,        only: run_CE
@@ -93,7 +93,7 @@ contains
       num_nuclides = 0
       num_macros = 0
       do i = 1, n_nuclides
-        if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then
+        if (get_awr_c(i) /= MACROSCOPIC_AWR) then
           num_nuclides = num_nuclides + 1
         else
           num_macros = num_macros + 1
@@ -118,12 +118,14 @@ contains
         nuc_names(i) = nuclides(i) % name
         awrs(i)     = nuclides(i) % awr
       else
-        if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then
-          nuc_names(j) = nuclides_MG(i) % obj % name
-          awrs(j)     = nuclides_MG(i) % obj % awr
+        if (get_awr_c(i) /= MACROSCOPIC_AWR) then
+          call get_name_c(i, len(nuc_names(j)), nuc_names(j))
+          nuc_names(j) = trim(nuc_names(j))
+          awrs(j)      = get_awr_c(i)
           j = j + 1
         else
-          macro_names(k) = nuclides_MG(i) % obj % name
+          call get_name_c(i, len(macro_names(k)), macro_names(k))
+          macro_names(k) = trim(macro_names(k))
           k = k + 1
         end if
       end if
@@ -353,7 +355,7 @@ contains
         num_nuclides  = 0
         num_macros = 0
         do j = 1, m % n_nuclides
-          if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then
+          if (get_awr_c(m % nuclide(j)) /= MACROSCOPIC_AWR) then
             num_nuclides = num_nuclides + 1
           else
             num_macros = num_macros + 1
@@ -379,12 +381,14 @@ contains
         k = 1
         n = 1
         do j = 1, m % n_nuclides
-          if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then
-            nuc_names(k) = nuclides_MG(m % nuclide(j)) % obj % name
+          if (get_awr_c(m % nuclide(j)) /= MACROSCOPIC_AWR) then
+            call get_name_c(m % nuclide(j), len(nuc_names(k)), nuc_names(k))
+            nuc_names(k) = trim(nuc_names(k))
             nuc_densities(k) = m % atom_density(j)
             k = k + 1
           else
-            macro_names(n) = nuclides_MG(m % nuclide(j)) % obj % name
+            call get_name_c(m % nuclide(j), len(macro_names(n)), macro_names(n))
+            macro_names(n) = trim(macro_names(n))
             n = n + 1
           end if
         end do

src/tallies/tally.F90

@@ -10,7 +10,7 @@ module tally
   use math,             only: t_percentile
   use mesh_header,      only: RegularMesh, meshes
   use message_passing
-  use mgxs_header
+  use mgxs_interface
   use nuclide_header
   use output,           only: header
   use particle_header,  only: LocalCoord, Particle
@@ -1229,8 +1229,6 @@ contains
     real(8) :: p_uvw(3)             ! Particle's current uvw
     integer :: p_g                  ! Particle group to use for getting info
                                     ! to tally with.
-    class(Mgxs), pointer :: matxs
-    class(Mgxs), pointer :: nucxs
 
     ! Set the direction and group to use with get_xs
     if (t % estimator == ESTIMATOR_ANALOG .or. &
@@ -1268,13 +1266,13 @@ contains
 
     ! To significantly reduce de-referencing, point matxs to the
     ! macroscopic Mgxs for the material of interest
-    matxs => macro_xs(p % material) % obj
+    call set_macro_angle_index_c(p % material, p_uvw)
 
     ! Do same for nucxs, point it to the microscopic nuclide data of interest
     if (i_nuclide > 0) then
-      nucxs => nuclides_MG(i_nuclide) % obj
       ! And since we haven't calculated this temperature index yet, do so now
-      call nucxs % find_temperature(p % sqrtkT)
+      call set_nuclide_temperature_index_c(i_nuclide, p % sqrtkT)
+      call set_nuclide_angle_index_c(i_nuclide, p_uvw)
     end if
 
     i = 0
@@ -1328,14 +1326,14 @@ contains
           end if
 
           if (i_nuclide > 0) then
-            score = score * atom_density * &
-                 nucxs % get_xs('total', p_g, UVW=p_uvw) / &
-                 matxs % get_xs('total', p_g, UVW=p_uvw) * flux
+            score = score * flux * atom_density * &
+                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_TOTAL, p_g)
           end if
 
         else
           if (i_nuclide > 0) then
-            score = nucxs % get_xs('total', p_g, UVW=p_uvw) * &
+            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) * &
                  atom_density * flux
           else
             score = material_xs % total * flux
@@ -1358,19 +1356,23 @@ contains
           end if
 
           if (i_nuclide > 0) then
-            score = score * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) &
-                 / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
+            score = score * flux * get_nuclide_xs_c(i_nuclide, &
+                 MG_GET_XS_INVERSE_VELOCITY, p_g) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
           else
-            score = score * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) &
-                 / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
+            score = score * flux * get_macro_xs_c(p % material, &
+                 MG_GET_XS_INVERSE_VELOCITY, p_g) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
           end if
 
         else
 
           if (i_nuclide > 0) then
-            score = flux * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw)
+            score = flux * get_nuclide_xs_c(i_nuclide, &
+                 MG_GET_XS_INVERSE_VELOCITY, p_g)
           else
-            score = flux * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw)
+            score = flux * get_macro_xs_c(p % material, &
+                 MG_GET_XS_INVERSE_VELOCITY, p_g)
           end if
         end if
 
@@ -1392,21 +1394,23 @@ contains
           ! adjust the score by the actual probability for that nuclide.
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 nucxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, &
-                                UVW=p_uvw, MU=p % mu) / &
-                 matxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, &
-                                UVW=p_uvw, MU=p % mu)
+                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU_MULT, &
+                                  p % last_g, p % g, MU=p % mu) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU_MULT, &
+                                p % last_g, p % g, MU=p % mu)
           end if
 
         else
           if (i_nuclide > 0) then
             score = atom_density * flux * &
-                 nucxs % get_xs('scatter/mult', p_g, UVW=p_uvw)
+                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_MULT, &
+                                  p_g, MU=p % mu)
           else
             ! Get the scattering x/s and take away
             ! the multiplication baked in to sigS
             score = flux * &
-                 matxs % get_xs('scatter/mult', p_g, UVW=p_uvw)
+                 get_macro_xs_c(p % material, MG_GET_XS_SCATTER_MULT, &
+                                p_g, MU=p % mu)
           end if
         end if
 
@@ -1428,19 +1432,20 @@ contains
           ! adjust the score by the actual probability for that nuclide.
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 nucxs % get_xs('scatter*f_mu', p % last_g, p % g, &
-                                UVW=p_uvw, MU=p % mu) / &
-                 matxs % get_xs('scatter*f_mu', p % last_g, p % g, &
-                                UVW=p_uvw, MU=p % mu)
+                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU, &
+                                p % last_g, p % g, MU=p % mu) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU, &
+                                p % last_g, p % g, MU=p % mu)
           end if
 
         else
           if (i_nuclide > 0) then
-              score = nucxs % get_xs('scatter', p_g, UVW=p_uvw) * &
-                   atom_density * flux
+            score = atom_density * flux * &
+                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER, p_g)
           else
             ! Get the scattering x/s, which includes multiplication
-            score = matxs % get_xs('scatter', p_g, UVW=p_uvw) * flux
+            score = flux * &
+                   get_macro_xs_c(p % material, MG_GET_XS_SCATTER, p_g)
           end if
         end if
 
@@ -1460,13 +1465,13 @@ contains
           end if
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 nucxs % get_xs('absorption', p_g, UVW=p_uvw) / &
-                 matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
           end if
         else
           if (i_nuclide > 0) then
-            score = nucxs % get_xs('absorption', p_g, UVW=p_uvw) * &
-                 atom_density * flux
+            score = atom_density * flux * &
+                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g)
           else
             score = material_xs % absorption * flux
           end if
@@ -1491,19 +1496,19 @@ contains
           end if
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
-                 matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
           else
             score = score * &
-                 matxs % get_xs('fission', p_g, UVW=p_uvw) / &
-                 matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                 get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
           end if
         else
           if (i_nuclide > 0) then
-            score = nucxs % get_xs('fission', p_g, UVW=p_uvw) * &
+            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) * &
                  atom_density * flux
           else
-            score = matxs % get_xs('fission', p_g, UVW=p_uvw) * flux
+            score = get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
           end if
         end if
 
@@ -1529,12 +1534,12 @@ contains
             score = p % absorb_wgt * flux
             if (i_nuclide > 0) then
               score = score * atom_density * &
-                   nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) / &
-                   matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
             else
               score = score * &
-                   matxs % get_xs('nu-fission', p_g, UVW=p_uvw) / &
-                   matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                   get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
             end if
           else
             ! Skip any non-fission events
@@ -1547,17 +1552,17 @@ contains
             score = keff * p % wgt_bank * flux
             if (i_nuclide > 0) then
               score = score * atom_density * &
-                   nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
-                   matxs % get_xs('fission', p_g, UVW=p_uvw)
+                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
             end if
           end if
 
         else
           if (i_nuclide > 0) then
-            score = nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) * &
+            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * &
                  atom_density * flux
           else
-            score = matxs % get_xs('nu-fission', p_g, UVW=p_uvw) * flux
+            score = get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) * flux
           end if
         end if
 
@@ -1583,12 +1588,12 @@ contains
             score = p % absorb_wgt * flux
             if (i_nuclide > 0) then
               score = score * atom_density * &
-                   nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / &
-                   matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
             else
               score = score * &
-                   matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / &
-                   matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                   get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
             end if
           else
             ! Skip any non-fission events
@@ -1602,17 +1607,17 @@ contains
                  / real(p % n_bank, 8)) * flux
             if (i_nuclide > 0) then
               score = score * atom_density * &
-                   nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
-                   matxs % get_xs('fission', p_g, UVW=p_uvw)
+                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
             end if
           end if
 
         else
           if (i_nuclide > 0) then
-            score = nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * &
+            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * &
                  atom_density * flux
           else
-            score = matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * flux
+            score = get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux
           end if
         end if
 
@@ -1638,7 +1643,7 @@ contains
             ! No fission events occur if survival biasing is on -- need to
             ! calculate fraction of absorptions that would have resulted in
             ! nu-fission
-            if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then
+            if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
 
               if (dg_filter > 0) then
                 select type(filt => filters(t % filter(dg_filter)) % obj)
@@ -1653,13 +1658,13 @@ contains
 
                     score = p % absorb_wgt * flux
                     if (i_nuclide > 0) then
-                      score = score * nucxs % get_xs('delayed-nu-fission', &
-                           p_g, UVW=p_uvw, dg=d) / &
-                           matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                      score = score * &
+                           get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                           get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
                     else
-                      score = score * matxs % get_xs('delayed-nu-fission', &
-                           p_g, UVW=p_uvw, dg=d) / &
-                           matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                      score = score * &
+                           get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                           get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
                     end if
 
                     call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1669,11 +1674,13 @@ contains
               else
                 score = p % absorb_wgt * flux
                 if (i_nuclide > 0) then
-                  score = score * nucxs % get_xs('delayed-nu-fission', p_g, &
-                       UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                  score = score * &
+                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
+                       get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
                 else
-                  score = score * matxs % get_xs('delayed-nu-fission', p_g, &
-                       UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                  score = score * &
+                       get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
+                       get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
                 end if
               end if
             end if
@@ -1703,8 +1710,8 @@ contains
 
                   if (i_nuclide > 0) then
                     score = score * atom_density * &
-                         nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
-                         matxs % get_xs('fission', p_g, UVW=p_uvw)
+                         get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
+                         get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
                   end if
 
                   call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1715,8 +1722,8 @@ contains
               score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank) * flux
               if (i_nuclide > 0) then
                 score = score * atom_density * &
-                     nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
-                     matxs % get_xs('fission', p_g, UVW=p_uvw)
+                     get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
+                     get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
               end if
             end if
           end if
@@ -1735,11 +1742,11 @@ contains
                 d = filt % groups(d_bin)
 
                 if (i_nuclide > 0) then
-                  score = nucxs % get_xs('delayed-nu-fission', p_g, &
-                       UVW=p_uvw, dg=d) * atom_density * flux
+                  score = atom_density * flux * &
+                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
                 else
-                  score = matxs % get_xs('delayed-nu-fission', p_g, &
-                       UVW=p_uvw, dg=d) * flux
+                  score = flux * &
+                       get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
                 end if
 
                 call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1748,11 +1755,12 @@ contains
             end select
           else
             if (i_nuclide > 0) then
-              score = nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) &
-                   * atom_density * flux
+              score = atom_density * flux * &
+                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g)
+
             else
-              score = matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) &
-                   * flux
+              score = flux * &
+                   get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g)
             end if
           end if
         end if
@@ -1767,7 +1775,7 @@ contains
             ! No fission events occur if survival biasing is on -- need to
             ! calculate fraction of absorptions that would have resulted in
             ! nu-fission
-            if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then
+            if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
 
               if (dg_filter > 0) then
                 select type(filt => filters(t % filter(dg_filter)) % obj)
@@ -1782,17 +1790,15 @@ contains
 
                     score = p % absorb_wgt * flux
                     if (i_nuclide > 0) then
-                      score = score * nucxs % get_xs('decay rate', p_g, &
-                           UVW=p_uvw, dg=d) * &
-                           nucxs % get_xs('delayed-nu-fission', p_g, &
-                           UVW=p_uvw, dg=d) / matxs % get_xs('absorption', &
-                           p_g, UVW=p_uvw)
+                      score = score * &
+                           get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                           get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                           get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
                     else
-                      score = score * matxs % get_xs('decay rate', p_g, &
-                           UVW=p_uvw, dg=d) * &
-                           matxs % get_xs('delayed-nu-fission', p_g, &
-                           UVW=p_uvw, dg=d) / matxs % get_xs('absorption', &
-                           p_g, UVW=p_uvw)
+                      score = score * &
+                           get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                           get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                           get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
                     end if
 
                     call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1809,15 +1815,15 @@ contains
                 ! for all delayed groups.
                 do d = 1, num_delayed_groups
                   if (i_nuclide > 0) then
-                    score = score + p % absorb_wgt * &
-                         nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
-                         nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
-                         dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
+                    score = score + p % absorb_wgt * flux * &
+                         get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                         get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                         get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
                   else
-                    score = score + p % absorb_wgt * &
-                         matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
-                         matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
-                         dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
+                    score = score + p % absorb_wgt * flux * &
+                         get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                         get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                         get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
                   end if
                 end do
               end if
@@ -1846,13 +1852,13 @@ contains
                 if (i_nuclide > 0) then
                   score = score + keff * atom_density * &
                        fission_bank(n_bank - p % n_bank + k) % wgt * &
-                       nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) * &
-                       nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
-                       matxs % get_xs('fission', p_g, UVW=p_uvw) * flux
+                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
+                       get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
                 else
                   score = score + keff * &
                        fission_bank(n_bank - p % n_bank + k) % wgt * &
-                       matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) * flux
+                       get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
                 end if
 
                 ! if the delayed group filter is present, tally to corresponding
@@ -1904,13 +1910,13 @@ contains
                 d = filt % groups(d_bin)
 
                 if (i_nuclide > 0) then
-                  score = nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
-                       nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
-                       dg=d) * atom_density * flux
+                  score = atom_density * flux * &
+                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
                 else
-                  score = matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
-                       matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
-                       dg=d) * flux
+                  score = flux * &
+                       get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                       get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
                 end if
 
                 call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1927,12 +1933,12 @@ contains
             do d = 1, num_delayed_groups
               if (i_nuclide > 0) then
                 score = score + atom_density * flux * &
-                     nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
-                     nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, dg=d)
+                     get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                     get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
               else
                 score = score + flux * &
-                     matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
-                     matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, dg=d)
+                     get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                     get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
               end if
             end do
           end if
@@ -1957,19 +1963,20 @@ contains
           end if
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / &
-                 matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
           else
             score = score * &
-                 matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / &
-                 matxs % get_xs('absorption', p_g, UVW=p_uvw)
+                 get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g) / &
+                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
           end if
         else
           if (i_nuclide > 0) then
-            score = nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * &
+            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * &
                  atom_density * flux
           else
-            score = matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * flux
+            score = flux * &
+                 get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g)
 
           end if
         end if
@@ -1988,8 +1995,6 @@ contains
            t % results(RESULT_VALUE, score_index, filter_index) + score
 
     end do SCORE_LOOP
-
-    nullify(matxs, nucxs)
   end subroutine score_general_mg
 
 !===============================================================================

src/tallies/tally_filter_energy.F90

@@ -6,7 +6,7 @@ module tally_filter_energy
   use constants
   use error
   use hdf5_interface
-  use mgxs_header,         only: num_energy_groups, rev_energy_bins
+  use mgxs_interface,      only: num_energy_groups, rev_energy_bins
   use particle_header,     only: Particle
   use settings,            only: run_CE
   use string,              only: to_str

src/tracking.F90

@@ -1,5 +1,7 @@
 module tracking
 
+  use, intrinsic :: ISO_C_BINDING
+
   use constants
   use error,              only: warning, write_message
   use geometry_header,    only: cells
@@ -7,7 +9,7 @@ module tracking
                                 check_cell_overlap
   use material_header,    only: materials, Material
   use message_passing
-  use mgxs_header
+  use mgxs_interface
   use nuclide_header
   use particle_header,    only: LocalCoord, Particle
   use physics,            only: collision
@@ -112,8 +114,10 @@ contains
           end if
         else
           ! Get the MG data
-          call macro_xs(p % material) % obj % calculate_xs(p % g, p % sqrtkT, &
-               p % coord(p % n_coord) % uvw, material_xs)
+          call calculate_xs_c(p % material, p % g, p % sqrtkT, &
+               p % coord(p % n_coord) % uvw, material_xs % total, &
+               material_xs % absorption, material_xs % nu_fission)
+
 
           ! Finally, update the particle group while we have already checked
           ! for if multi-group

src/xsdata.cpp

@@ -0,0 +1,713 @@
+#include <cmath>
+#include <cstdlib>
+#include <algorithm>
+#include <numeric>
+
+#include "constants.h"
+#include "error.h"
+#include "math_functions.h"
+#include "random_lcg.h"
+#include "xsdata.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// XsData class methods
+//==============================================================================
+
+XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable,
+     int scatter_format, int n_pol, int n_azi)
+{
+  int n_ang = n_pol * n_azi;
+
+  // check to make sure scatter format is OK before we allocate
+  if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR &&
+      scatter_format != ANGLE_LEGENDRE) {
+    fatal_error("Invalid scatter_format!");
+  }
+  // allocate all [temperature][phi][theta][in group] quantities
+  total = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
+  absorption = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
+  inverse_velocity = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
+  if (fissionable) {
+    fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
+    nu_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
+    prompt_nu_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
+    kappa_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
+  }
+
+  // allocate decay_rate; [temperature][phi][theta][delayed group]
+  decay_rate = double_2dvec(n_ang, double_1dvec(num_delayed_groups, 0.));
+
+  if (fissionable) {
+    // allocate delayed_nu_fission; [temperature][phi][theta][in group][delay group]
+    delayed_nu_fission = double_3dvec(n_ang, double_2dvec(energy_groups,
+         double_1dvec(num_delayed_groups, 0.)));
+
+    // chi_prompt; [temperature][phi][theta][in group][delayed group]
+    chi_prompt = double_3dvec(n_ang, double_2dvec(energy_groups,
+         double_1dvec(energy_groups, 0.)));
+
+    // chi_delayed; [temperature][phi][theta][in group][out group][delay group]
+    chi_delayed = double_4dvec(n_ang, double_3dvec(energy_groups,
+         double_2dvec(energy_groups, double_1dvec(num_delayed_groups, 0.))));
+  }
+
+
+  for (int a = 0; a < n_ang; a++) {
+    if (scatter_format == ANGLE_HISTOGRAM) {
+      // scatter[a] = std::make_unique(ScattDataHistogram);
+      scatter.emplace_back(new ScattDataHistogram);
+    } else if (scatter_format == ANGLE_TABULAR) {
+      // scatter[a] = std::make_unique(ScattDataTabular);
+      scatter.emplace_back(new ScattDataTabular);
+    } else if (scatter_format == ANGLE_LEGENDRE) {
+      // scatter[a] = std::make_unique(ScattDataLegendre);
+      scatter.emplace_back(new ScattDataLegendre);
+    }
+  }
+}
+
+//==============================================================================
+
+void
+XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
+     int final_scatter_format, int order_data, int max_order,
+     int legendre_to_tabular_points, bool is_isotropic, int n_pol, int n_azi)
+{
+  // Reconstruct the dimension information so it doesn't need to be passed
+  int n_ang = n_pol * n_azi;
+  int energy_groups = total[0].size();
+  int delayed_groups = decay_rate[0].size();
+
+  // Set the fissionable-specific data
+  if (fissionable) {
+    fission_from_hdf5(xsdata_grp, n_pol, n_azi, energy_groups, delayed_groups,
+                      is_isotropic);
+  }
+  // Get the non-fission-specific data
+  read_nd_vector(xsdata_grp, "decay_rate", decay_rate);
+  read_nd_vector(xsdata_grp, "absorption", absorption, true);
+  read_nd_vector(xsdata_grp, "inverse-velocity", inverse_velocity);
+
+  // Get scattering data
+  scatter_from_hdf5(xsdata_grp, n_pol, n_azi, energy_groups, scatter_format,
+       final_scatter_format, order_data, max_order, legendre_to_tabular_points);
+
+  // Check absorption to ensure it is not 0 since it is often the
+  // denominator in tally methods
+  for (int a = 0; a < n_ang; a++) {
+    for (int gin = 0; gin < energy_groups; gin++) {
+      if (absorption[a][gin] == 0.) absorption[a][gin] = 1.e-10;
+    }
+  }
+
+  // Get or calculate the total x/s
+  if (object_exists(xsdata_grp, "total")) {
+    read_nd_vector(xsdata_grp, "total", total);
+  } else {
+    for (int a = 0; a < n_ang; a++) {
+      for (int gin = 0; gin < energy_groups; gin++) {
+        total[a][gin] = absorption[a][gin] + scatter[a]->scattxs[gin];
+      }
+    }
+  }
+
+  // Fix if total is 0, since it is in the denominator when tallying
+  for (int a = 0; a < n_ang; a++) {
+    for (int gin = 0; gin < energy_groups; gin++) {
+      if (total[a][gin] == 0.) total[a][gin] = 1.e-10;
+    }
+  }
+}
+
+//==============================================================================
+
+void
+XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
+     int energy_groups, int delayed_groups, bool is_isotropic)
+{
+  int n_ang = n_pol * n_azi;
+  // Get the fission and kappa_fission data xs; these are optional
+  read_nd_vector(xsdata_grp, "fission", fission);
+  read_nd_vector(xsdata_grp, "kappa-fission", kappa_fission);
+
+  // Set/get beta
+  double_3dvec temp_beta =double_3dvec(n_ang, double_2dvec(energy_groups,
+       double_1dvec(delayed_groups, 0.)));
+  if (object_exists(xsdata_grp, "beta")) {
+    hid_t xsdata = open_dataset(xsdata_grp, "beta");
+    int ndims = dataset_ndims(xsdata);
+
+    // raise ndims to make the isotropic ndims the same as angular
+    if (is_isotropic) ndims += 2;
+
+    if (ndims == 3) {
+      // Beta is input as [delayed group]
+      double_1dvec temp_arr(n_pol * n_azi * delayed_groups);
+      read_nd_vector(xsdata_grp, "beta", temp_arr);
+
+      // Broadcast to all incoming groups
+      int temp_idx = 0;
+      for (int a = 0; a < n_ang; a++) {
+        for (int dg = 0; dg < delayed_groups; dg++) {
+          // Set the first group index and copy the rest
+          temp_beta[a][0][dg] = temp_arr[temp_idx++];
+          for (int gin = 1; gin < energy_groups; gin++) {
+            temp_beta[a][gin] = temp_beta[a][0];
+          }
+        }
+      }
+    } else if (ndims == 4) {
+      // Beta is input as [in group][delayed group]
+      read_nd_vector(xsdata_grp, "beta", temp_beta);
+    } else {
+      fatal_error("beta must be provided as a 3D or 4D array!");
+    }
+  }
+
+  // If chi is provided, set chi-prompt and chi-delayed
+  if (object_exists(xsdata_grp, "chi")) {
+    double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
+    read_nd_vector(xsdata_grp, "chi", temp_arr);
+
+    for (int a = 0; a < n_ang; a++) {
+      // First set the first group
+      for (int gout = 0; gout < energy_groups; gout++) {
+        chi_prompt[a][0][gout] = temp_arr[a][gout];
+      }
+
+      // Now normalize this data
+      double chi_sum = std::accumulate(chi_prompt[a][0].begin(),
+                                       chi_prompt[a][0].end(),
+                                       0.);
+      if (chi_sum <= 0.) {
+        fatal_error("Encountered chi for a group that is <= 0!");
+      }
+      for (int gout = 0; gout < energy_groups; gout++) {
+        chi_prompt[a][0][gout] /= chi_sum;
+      }
+
+      // And extend to the remaining incoming groups
+      for (int gin = 1; gin < energy_groups; gin++) {
+        chi_prompt[a][gin] = chi_prompt[a][0];
+      }
+
+      // Finally set chi-delayed equal to chi-prompt
+      // Set chi-delayed to chi-prompt
+      for(int gin = 0; gin < energy_groups; gin++) {
+        for (int gout = 0; gout < energy_groups; gout++) {
+          for (int dg = 0; dg < delayed_groups; dg++) {
+            chi_delayed[a][gin][gout][dg] =
+                 chi_prompt[a][gin][gout];
+          }
+        }
+      }
+    }
+  }
+
+  // If nu-fission is provided, set prompt- and delayed-nu-fission;
+  // if nu-fission is a matrix, set chi-prompt and chi-delayed.
+  if (object_exists(xsdata_grp, "nu-fission")) {
+    hid_t xsdata = open_dataset(xsdata_grp, "nu-fission");
+    int ndims = dataset_ndims(xsdata);
+    // raise ndims to make the isotropic ndims the same as angular
+    if (is_isotropic) ndims += 2;
+
+    if (ndims == 3) {
+      // nu-fission is a 3-d array
+      read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission);
+
+      // set delayed-nu-fission and correct prompt-nu-fission with beta
+      for (int a = 0; a < n_ang; a++) {
+        for (int gin = 0; gin < energy_groups; gin++) {
+          for (int dg = 0; dg < delayed_groups; dg++) {
+            delayed_nu_fission[a][gin][dg] =
+                 temp_beta[a][gin][dg] * prompt_nu_fission[a][gin];
+          }
+
+          // Correct the prompt-nu-fission using the delayed neutron fraction
+          if (delayed_groups > 0) {
+            double beta_sum = std::accumulate(temp_beta[a][gin].begin(),
+                                              temp_beta[a][gin].end(), 0.);
+            prompt_nu_fission[a][gin] *= (1. - beta_sum);
+          }
+        }
+      }
+
+    } else if (ndims == 4) {
+      // nu-fission is a matrix
+      read_nd_vector(xsdata_grp, "nu-fission", chi_prompt);
+
+      // Normalize the chi info so the CDF is 1.
+      for (int a = 0; a < n_ang; a++) {
+        for (int gin = 0; gin < energy_groups; gin++) {
+          double chi_sum = std::accumulate(chi_prompt[a][gin].begin(),
+                                           chi_prompt[a][gin].end(), 0.);
+          // Set the vector nu-fission from the matrix nu-fission
+          prompt_nu_fission[a][gin] = chi_sum;
+
+          if (chi_sum >= 0.) {
+            for (int gout = 0; gout < energy_groups; gout++) {
+              chi_prompt[a][gin][gout] /= chi_sum;
+            }
+          } else {
+            fatal_error("Encountered chi for a group that is <= 0!");
+          }
+        }
+
+        // set chi-delayed to chi-prompt
+        for (int gin = 0; gin < energy_groups; gin++) {
+          for (int gout = 0; gout < energy_groups; gout++) {
+            for (int dg = 0; dg < delayed_groups; dg++) {
+              chi_delayed[a][gin][gout][dg] =
+                   chi_prompt[a][gin][gout];
+            }
+          }
+        }
+
+        // Set the delayed-nu-fission and correct prompt-nu-fission with beta
+        for (int gin = 0; gin < energy_groups; gin++) {
+          for (int dg = 0; dg < delayed_groups; dg++) {
+            delayed_nu_fission[a][gin][dg] =
+                 temp_beta[a][gin][dg] *
+                 prompt_nu_fission[a][gin];
+          }
+
+          // Correct prompt-nu-fission using the delayed neutron fraction
+          if (delayed_groups > 0) {
+            double beta_sum = std::accumulate(temp_beta[a][gin].begin(),
+                                              temp_beta[a][gin].end(), 0.);
+            prompt_nu_fission[a][gin] *= (1. - beta_sum);
+          }
+        }
+      }
+    } else {
+      fatal_error("nu-fission must be provided as a 3D or 4D array!");
+    }
+
+    close_dataset(xsdata);
+  }
+
+  // If chi-prompt is provided, set chi-prompt
+  if (object_exists(xsdata_grp, "chi-prompt")) {
+    double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
+    read_nd_vector(xsdata_grp, "chi-prompt", temp_arr);
+
+    for (int a = 0; a < n_ang; a++) {
+      for (int gin = 0; gin < energy_groups; gin++) {
+        for (int gout = 0; gout < energy_groups; gout++) {
+          chi_prompt[a][gin][gout] = temp_arr[a][gout];
+        }
+
+        // Normalize chi so its CDF goes to 1
+        double chi_sum = std::accumulate(chi_prompt[a][gin].begin(),
+                                         chi_prompt[a][gin].end(), 0.);
+        if (chi_sum >= 0.) {
+          for (int gout = 0; gout < energy_groups; gout++) {
+            chi_prompt[a][gin][gout] /= chi_sum;
+          }
+        } else {
+          fatal_error("Encountered chi-prompt for a group that is <= 0.!");
+        }
+      }
+    }
+  }
+
+  // If chi-delayed is provided, set chi-delayed
+  if (object_exists(xsdata_grp, "chi-delayed")) {
+    hid_t xsdata = open_dataset(xsdata_grp, "chi-delayed");
+    int ndims = dataset_ndims(xsdata);
+    // raise ndims to make the isotropic ndims the same as angular
+    if (is_isotropic) ndims += 2;
+    close_dataset(xsdata);
+
+    if (ndims == 3) {
+      // chi-delayed is a [in group] vector
+      double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
+      read_nd_vector(xsdata_grp, "chi-delayed", temp_arr);
+
+      for (int a = 0; a < n_ang; a++) {
+        // normalize the chi CDF to 1
+        double chi_sum = std::accumulate(temp_arr[a].begin(),
+                                         temp_arr[a].end(), 0.);
+        if (chi_sum <= 0.) {
+          fatal_error("Encountered chi-delayed for a group that is <= 0!");
+        }
+
+        // set chi-delayed
+        for (int gin = 0; gin < energy_groups; gin++) {
+          for (int gout = 0; gout < energy_groups; gout++) {
+            for (int dg = 0; dg < delayed_groups; dg++) {
+              chi_delayed[a][gin][gout][dg] = temp_arr[a][gout] / chi_sum;
+            }
+          }
+        }
+      }
+    } else if (ndims == 4) {
+      // chi_delayed is a matrix
+      read_nd_vector(xsdata_grp, "chi-delayed", chi_delayed);
+
+      // Normalize the chi info so the CDF is 1.
+      for (int a = 0; a < n_ang; a++) {
+        for (int dg = 0; dg < delayed_groups; dg++) {
+          for (int gin = 0; gin < energy_groups; gin++) {
+            double chi_sum = 0.;
+            for (int gout = 0; gout < energy_groups; gout++) {
+              chi_sum += chi_delayed[a][gin][gout][dg];
+            }
+
+            if (chi_sum > 0.) {
+              for (int gout = 0; gout < energy_groups; gout++) {
+                chi_delayed[a][gin][gout][dg] /= chi_sum;
+              }
+            } else {
+              fatal_error("Encountered chi-delayed for a group that is <= 0!");
+            }
+          }
+        }
+      }
+    } else {
+      fatal_error("chi-delayed must be provided as a 3D or 4D array!");
+    }
+  }
+
+  // Get prompt-nu-fission, if present
+  if (object_exists(xsdata_grp, "prompt-nu-fission")) {
+    hid_t xsdata = open_dataset(xsdata_grp, "prompt-nu-fission");
+    int ndims = dataset_ndims(xsdata);
+    // raise ndims to make the isotropic ndims the same as angular
+    if (is_isotropic) ndims += 2;
+    close_dataset(xsdata);
+
+    if (ndims == 3) {
+      // prompt-nu-fission is a [in group] vector
+      read_nd_vector(xsdata_grp, "prompt-nu-fission",
+                     prompt_nu_fission);
+    } else if (ndims == 4) {
+      // prompt nu fission is a matrix,
+      // so set prompt_nu_fiss & chi_prompt
+      double_3dvec temp_arr(n_ang, double_2dvec(energy_groups,
+           double_1dvec(energy_groups)));
+      read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_arr);
+
+      // The prompt_nu_fission vector from the matrix form
+      for (int a = 0; a < n_ang; a++) {
+        for (int gin = 0; gin < energy_groups; gin++) {
+          double prompt_sum = std::accumulate(temp_arr[a][gin].begin(),
+                                              temp_arr[a][gin].end(), 0.);
+          prompt_nu_fission[a][gin] = prompt_sum;
+        }
+
+        // The chi_prompt data is just the normalized fission matrix
+        for (int gin= 0; gin < energy_groups; gin++) {
+          if (prompt_nu_fission[a][gin] > 0.) {
+            for (int gout = 0; gout < energy_groups; gout++) {
+              chi_prompt[a][gin][gout] =
+                   temp_arr[a][gin][gout] / prompt_nu_fission[a][gin];
+            }
+          } else {
+            fatal_error("Encountered chi-prompt for a group that is <= 0!");
+          }
+        }
+      }
+
+    } else {
+      fatal_error("prompt-nu-fission must be provided as a 3D or 4D array!");
+    }
+  }
+
+  // Get delayed-nu-fission, if present
+  if (object_exists(xsdata_grp, "delayed-nu-fission")) {
+    hid_t xsdata = open_dataset(xsdata_grp, "delayed-nu-fission");
+    int ndims = dataset_ndims(xsdata);
+    close_dataset(xsdata);
+    // raise ndims to make the isotropic ndims the same as angular
+    if (is_isotropic) ndims += 2;
+
+    if (ndims == 3) {
+      // delayed-nu-fission is an [in group] vector
+      if (temp_beta[0][0][0] == 0.) {
+        fatal_error("cannot set delayed-nu-fission with a 1D array if "
+                    "beta is not provided");
+      }
+      double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
+      read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
+
+      for (int a = 0; a < n_ang; a++) {
+        for (int gin = 0; gin < energy_groups; gin++) {
+          for (int dg = 0; dg < delayed_groups; dg++) {
+            // Set delayed-nu-fission using beta
+            delayed_nu_fission[a][gin][dg] =
+                 temp_beta[a][gin][dg] * temp_arr[a][gin];
+          }
+        }
+      }
+
+    } else if (ndims == 4) {
+      read_nd_vector(xsdata_grp, "delayed-nu-fission",
+                    delayed_nu_fission);
+
+    } else if (ndims == 5) {
+      // This will contain delayed-nu-fision and chi-delayed data
+      double_4dvec temp_arr(n_ang, double_3dvec(energy_groups,
+           double_2dvec(energy_groups, double_1dvec(delayed_groups))));
+      read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
+
+      // Set the 3D delayed-nu-fission matrix and 4D chi-delayed matrix
+      // from the 4D delayed-nu-fission matrix
+      for (int a = 0; a < n_ang; a++) {
+        for (int dg = 0; dg < delayed_groups; dg++) {
+          for (int gin = 0; gin < energy_groups; gin++) {
+            double gout_sum = 0.;
+            for (int gout = 0; gout < energy_groups; gout++) {
+              gout_sum += temp_arr[a][gin][gout][dg];
+              chi_delayed[a][gin][gout][dg] = temp_arr[a][gin][gout][dg];
+            }
+            delayed_nu_fission[a][gin][dg] = gout_sum;
+            // Normalize chi-delayed
+            if (gout_sum > 0.) {
+              for (int gout = 0; gout < energy_groups; gout++) {
+                chi_delayed[a][gin][gout][dg] /= gout_sum;
+              }
+            } else {
+              fatal_error("Encountered chi-delayed for a group that is <= 0!");
+            }
+          }
+        }
+      }
+
+    } else {
+      fatal_error("prompt-nu-fission must be provided as a 3D, 4D, or 5D "
+                  "array!");
+    }
+  }
+
+  // Combine prompt_nu_fission and delayed_nu_fission into nu_fission
+  for (int a = 0; a < n_ang; a++) {
+    for (int gin = 0; gin < energy_groups; gin++) {
+      nu_fission[a][gin] =
+           std::accumulate(delayed_nu_fission[a][gin].begin(),
+                           delayed_nu_fission[a][gin].end(),
+                           prompt_nu_fission[a][gin]);
+    }
+  }
+}
+
+//==============================================================================
+
+void
+XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
+     int energy_groups, int scatter_format, int final_scatter_format,
+     int order_data, int max_order, int legendre_to_tabular_points)
+{
+  int n_ang = n_pol * n_azi;
+  if (!object_exists(xsdata_grp, "scatter_data")) {
+    fatal_error("Must provide scatter_data group!");
+  }
+  hid_t scatt_grp = open_group(xsdata_grp, "scatter_data");
+
+  // Get the outgoing group boundary indices
+  int_2dvec gmin(n_ang, int_1dvec(energy_groups));
+  read_nd_vector(scatt_grp, "g_min", gmin, true);
+  int_2dvec gmax(n_ang, int_1dvec(energy_groups));
+  read_nd_vector(scatt_grp, "g_max", gmax, true);
+
+  // Make gmin and gmax start from 0 vice 1 as they do in the library
+  for (int a = 0; a < n_ang; a++) {
+    for (int gin = 0; gin < energy_groups; gin++) {
+      gmin[a][gin] -= 1;
+      gmax[a][gin] -= 1;
+    }
+  }
+
+  // Now use this info to find the length of a vector to hold the flattened
+  // data.
+  int length = 0;
+  for (int a = 0; a < n_ang; a++) {
+    for (int gin = 0; gin < energy_groups; gin++) {
+      length += order_data * (gmax[a][gin] - gmin[a][gin] + 1);
+    }
+  }
+  double_1dvec temp_arr(length);
+  read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true);
+
+  // Compare the number of orders given with the max order of the problem;
+  // strip off the superfluous orders if needed
+  int order_dim;
+  if (scatter_format == ANGLE_LEGENDRE) {
+    order_dim = std::min(order_data - 1, max_order) + 1;
+  } else {
+    order_dim = order_data;
+  }
+
+  // convert the flattened temp_arr to a jagged array for passing to
+  // scatt data
+  double_4dvec input_scatt(n_ang, double_3dvec(energy_groups));
+
+  int temp_idx = 0;
+  for (int a = 0; a < n_ang; a++) {
+    for (int gin = 0; gin < energy_groups; gin++) {
+      input_scatt[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1);
+      for (int i_gout = 0; i_gout < input_scatt[a][gin].size(); i_gout++) {
+        input_scatt[a][gin][i_gout].resize(order_dim);
+        for (int l = 0; l < order_dim; l++) {
+          input_scatt[a][gin][i_gout][l] = temp_arr[temp_idx++];
+        }
+        // Adjust index for the orders we didnt take
+        temp_idx += (order_data - order_dim);
+      }
+    }
+  }
+  temp_arr.clear();
+
+  // Get multiplication matrix
+  double_3dvec temp_mult(n_ang, double_2dvec(energy_groups));
+  if (object_exists(scatt_grp, "multiplicity_matrix")) {
+    temp_arr.resize(length / order_data);
+    read_nd_vector(scatt_grp, "multiplicity_matrix", temp_arr);
+
+    // convert the flat temp_arr to a jagged array for passing to scatt data
+    int temp_idx = 0;
+    for (int a = 0; a < n_ang; a++) {
+      for (int gin = 0; gin < energy_groups; gin++) {
+        temp_mult[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1);
+        for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
+          temp_mult[a][gin][i_gout] = temp_arr[temp_idx++];
+        }
+      }
+    }
+  } else {
+    // Use a default: multiplicities are 1.0.
+    for (int a = 0; a < n_ang; a++) {
+      for (int gin = 0; gin < energy_groups; gin++) {
+        temp_mult[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1);
+        for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
+          temp_mult[a][gin][i_gout] = 1.;
+        }
+      }
+    }
+  }
+  temp_arr.clear();
+  close_group(scatt_grp);
+
+  // Finally, convert the Legendre data to tabular, if needed
+  if (scatter_format == ANGLE_LEGENDRE &&
+      final_scatter_format == ANGLE_TABULAR) {
+    for (int a = 0; a < n_ang; a++) {
+      ScattDataLegendre legendre_scatt;
+      legendre_scatt.init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]);
+
+      // Now create a tabular version of legendre_scatt
+      convert_legendre_to_tabular(legendre_scatt,
+           *static_cast<ScattDataTabular*>(scatter[a].get()),
+           legendre_to_tabular_points);
+
+      scatter_format = final_scatter_format;
+    }
+  } else {
+    // We are sticking with the current representation
+    // Initialize the ScattData object with this data
+    for (int a = 0; a < n_ang; a++) {
+      scatter[a]->init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]);
+    }
+  }
+}
+
+//==============================================================================
+
+void
+XsData::combine(const std::vector<XsData*>& those_xs,
+                const double_1dvec& scalars)
+{
+  // Combine the non-scattering data
+  for (int i = 0; i < those_xs.size(); i++) {
+    XsData* that = those_xs[i];
+    if (!equiv(*that)) fatal_error("Cannot combine the XsData objects!");
+    double scalar = scalars[i];
+    for (int a = 0; a < total.size(); a++) {
+      for (int gin = 0; gin < total[a].size(); gin++) {
+        total[a][gin] += scalar * that->total[a][gin];
+        absorption[a][gin] += scalar * that->absorption[a][gin];
+        if (i == 0) {
+          inverse_velocity[a][gin] = that->inverse_velocity[a][gin];
+        }
+        if (that->prompt_nu_fission.size() > 0) {
+          nu_fission[a][gin] += scalar * that->nu_fission[a][gin];
+          prompt_nu_fission[a][gin] +=
+               scalar * that->prompt_nu_fission[a][gin];
+          kappa_fission[a][gin] += scalar * that->kappa_fission[a][gin];
+          fission[a][gin] += scalar * that->fission[a][gin];
+
+          for (int dg = 0; dg < delayed_nu_fission[a][gin].size(); dg++) {
+            delayed_nu_fission[a][gin][dg] +=
+                 scalar * that->delayed_nu_fission[a][gin][dg];
+          }
+
+          for (int gout = 0; gout < chi_prompt[a][gin].size(); gout++) {
+            chi_prompt[a][gin][gout] +=
+                 scalar * that->chi_prompt[a][gin][gout];
+
+            for (int dg = 0; dg < chi_delayed[a][gin][gout].size(); dg++) {
+              chi_delayed[a][gin][gout][dg] +=
+                   scalar * that->chi_delayed[a][gin][gout][dg];
+            }
+          }
+        }
+      }
+
+      for (int dg = 0; dg < decay_rate[a].size(); dg++) {
+        decay_rate[a][dg] += scalar * that->decay_rate[a][dg];
+      }
+
+      // Normalize chi
+      if (chi_prompt.size() > 0) {
+        for (int gin = 0; gin < chi_prompt[a].size(); gin++) {
+          double norm = std::accumulate(chi_prompt[a][gin].begin(),
+                                        chi_prompt[a][gin].end(), 0.);
+          if (norm > 0.) {
+            for (int gout = 0; gout < chi_prompt[a][gin].size(); gout++) {
+              chi_prompt[a][gin][gout] /= norm;
+            }
+          }
+
+          for (int dg = 0; dg < chi_delayed[a][gin][0].size(); dg++) {
+            norm = 0.;
+            for (int gout = 0; gout < chi_delayed[a][gin].size(); gout++) {
+              norm += chi_delayed[a][gin][gout][dg];
+            }
+            if (norm > 0.) {
+              for (int gout = 0; gout < chi_delayed[a][gin].size(); gout++) {
+                chi_delayed[a][gin][gout][dg] /= norm;
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  // Allow the ScattData object to combine itself
+  for (int a = 0; a < total.size(); a++) {
+    // Build vector of the scattering objects to incorporate
+    std::vector<ScattData*> those_scatts(those_xs.size());
+    for (int i = 0; i < those_xs.size(); i++) {
+      those_scatts[i] = those_xs[i]->scatter[a].get();
+    }
+
+    // Now combine these guys
+    scatter[a]->combine(those_scatts, scalars);
+  }
+}
+
+//==============================================================================
+
+bool
+XsData::equiv(const XsData& that)
+{
+  return ((absorption.size() == that.absorption.size()) &&
+      (absorption[0].size() == that.absorption[0].size()));
+}
+
+} //namespace openmc

src/xsdata.h

@@ -0,0 +1,118 @@
+//! \file xsdata.h
+//! A collection of classes for containing the Multi-Group Cross Section data
+
+#ifndef XSDATA_H
+#define XSDATA_H
+
+#include <memory>
+#include <vector>
+
+#include "hdf5_interface.h"
+#include "scattdata.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// XSDATA contains the temperature-independent cross section data for an MGXS
+//==============================================================================
+
+class XsData {
+
+  private:
+    //! \brief Reads scattering data from the HDF5 file
+    void
+    scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups,
+         int scatter_format, int final_scatter_format, int order_data,
+         int max_order, int legendre_to_tabular_points);
+
+    //! \brief Reads fission data from the HDF5 file
+    void
+    fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups,
+         int delayed_groups, bool is_isotropic);
+
+  public:
+
+    // The following quantities have the following dimensions:
+    // [angle][incoming group]
+    double_2dvec total;
+    double_2dvec absorption;
+    double_2dvec nu_fission;
+    double_2dvec prompt_nu_fission;
+    double_2dvec kappa_fission;
+    double_2dvec fission;
+    double_2dvec inverse_velocity;
+
+    // decay_rate has the following dimensions:
+    // [angle][delayed group]
+    double_2dvec decay_rate;
+    // delayed_nu_fission has the following dimensions:
+    // [angle][incoming group][delayed group]
+    double_3dvec delayed_nu_fission;
+    // chi_prompt has the following dimensions:
+    // [angle][incoming group][outgoing group]
+    double_3dvec chi_prompt;
+    // chi_delayed has the following dimensions:
+    // [angle][incoming group][outgoing group][delayed group]
+    double_4dvec chi_delayed;
+    // scatter has the following dimensions: [angle]
+    std::vector<std::shared_ptr<ScattData> > scatter;
+
+    XsData() = default;
+
+    //! \brief Constructs the XsData object metadata.
+    //!
+    //! @param num_groups Number of energy groups.
+    //! @param num_delayed_groups Number of delayed groups.
+    //! @param fissionable Is this a fissionable data set or not.
+    //! @param scatter_format The scattering representation of the file.
+    //! @param n_pol Number of polar angles.
+    //! @param n_azi Number of azimuthal angles.
+    XsData(int num_groups, int num_delayed_groups, bool fissionable,
+           int scatter_format, int n_pol, int n_azi);
+
+    //! \brief Loads the XsData object from the HDF5 file
+    //!
+    //! @param xs_id HDF5 group id for the cross section data.
+    //! @param fissionable Is this a fissionable data set or not.
+    //! @param scatter_format The scattering representation of the file.
+    //! @param final_scatter_format The scattering representation after reading;
+    //!   this is different from scatter_format if converting a Legendre to
+    //!   a tabular representation.
+    //! @param order_data The dimensionality of the scattering data in the file.
+    //! @param max_order Maximum order requested by the user;
+    //!   this is only used for Legendre scattering.
+    //! @param legendre_to_tabular Flag to denote if any Legendre provided
+    //!   should be converted to a Tabular representation.
+    //! @param legendre_to_tabular_points If a conversion is requested, this
+    //!   provides the number of points to use in the tabular representation.
+    //! @param is_isotropic Is this an isotropic or angular with respect to
+    //!   the incoming particle.
+    //! @param n_pol Number of polar angles.
+    //! @param n_azi Number of azimuthal angles.
+    void
+    from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
+         int final_scatter_format, int order_data, int max_order,
+         int legendre_to_tabular_points, bool is_isotropic, int n_pol,
+         int n_azi);
+
+    //! \brief Combines the microscopic data to a macroscopic object.
+    //!
+    //! @param micros Microscopic objects to combine.
+    //! @param scalars Scalars to multiply the microscopic data by.
+    void
+    combine(const std::vector<XsData*>& those_xs, const double_1dvec& scalars);
+
+    //! \brief Checks to see if this and that are able to be combined
+    //!
+    //! This comparison is used when building macroscopic cross sections
+    //! from microscopic cross sections.
+    //! @param that The other XsData to compare to this one.
+    //! @return True if they can be combined.
+    bool
+    equiv(const XsData& that);
+};
+
+
+} //namespace openmc
+#endif // XSDATA_H