Initial implementation of MGXS C++ code

src/constants.h

@@ -0,0 +1,60 @@
+
+#ifndef CONSTANTS_H
+#define CONSTANTS_H
+
+#include <array>
+#include <vector>
+
+namespace openmc {
+
+typedef std::array<double, 3> dir_arr;
+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;
+
+int constexpr 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};
+
+
+} // namespace openmc
+
+#endif // CONSTANTS_H

src/hdf5_interface.cpp

@@ -447,6 +447,152 @@ 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[0], 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();
+    std::vector<double> temp_arr = std::vector<double>(dim1 * dim2);
+    read_double(obj_id, name, &temp_arr[0], 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)) {
+    dim1 = result.size();
+    dim2 = result[0].size();
+    dim3 = result[0][0].size();
+    std::vector<double> temp_arr = std::vector<double>(dim1 * dim2 * dim3);
+    read_double(obj_id, name, &temp_arr[0], 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)) {
+    dim1 = result.size();
+    dim2 = result[0].size();
+    dim3 = result[0][0].size();
+    std::vector<int> temp_arr = std::vector<int>(dim1 * dim2 * dim3);
+    read_int(obj_id, name, &temp_arr[0], 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)) {
+    dim1 = result.size();
+    dim2 = result[0].size();
+    dim3 = result[0][0].size();
+    dim4 = result[0][0][0].size();
+    std::vector<double> temp_arr = std::vector<double>(
+         dim1 * dim2 * dim3 * dim4);
+    read_double(obj_id, name, &temp_arr[0], 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)) {
+    dim1 = result.size();
+    dim2 = result[0].size();
+    dim3 = result[0][0].size();
+    dim4 = result[0][0][0].size();
+    dim5 = result[0][0][0][0].size();
+    std::vector<double> temp_arr = std::vector<double>(
+         dim1 * dim2 * dim3 * dim4 * dim5);
+    read_double(obj_id, name, &temp_arr[0], 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

@@ -55,6 +55,35 @@ 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<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/mgxs.h

@@ -0,0 +1,95 @@
+//! \file mgxs.h
+//! A collection of classes for Multi-Group Cross Section data
+
+#ifndef MGXS_H
+#define MGXS_H
+
+#include <cmath>
+#include <cstdlib>
+#include <algorithm>
+#include <string>
+#include <valarray>
+#include <vector>
+
+#include "constants.h"
+#include "hdf5_interface.h"
+#include "math_functions.h"
+#include "random_lcg.h"
+#include "scattdata.h"
+#include "string_functions.h"
+#include "xsdata.h"
+
+
+namespace openmc {
+
+
+//==============================================================================
+// MGXS contains the mgxs data for a nuclide/material
+//==============================================================================
+
+class Mgxs {
+  private:
+    std::string name;   // name of dataset, e.g., UO2
+    double awr;         // atomic weight ratio
+    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
+    int index_temp;     // cache of temperature index
+    double last_sqrtkT; // cache of the temperature corresponding to index_temp
+    std::vector<XsData> xs; // Cross section data
+    int n_pol;
+    int n_azi;
+    int index_pol; // cache fof the angle indices
+    int index_azi;
+    double_1dvec polar;
+    double_1dvec azimuthal;
+    dir_arr last_uvw;
+    void _metadata_from_hdf5(hid_t xs_id, int in_num_groups,
+         int in_num_delayed_groups, double_1dvec temperature, int& method,
+         double tolerance, double_1dvec& temps_to_read, int& order_dim);
+
+  public:
+    bool fissionable;   // Is this fissionable
+    void init(const std::string& in_name, double in_awr, double_1dvec& in_kTs,
+              bool in_fissionable, int in_scatter_format, int in_num_groups,
+              int in_num_delayed_groups, double_1dvec& in_polar,
+              double_1dvec& in_azimuthal);
+    void build_macro(const std::string& in_name, double_1dvec& mat_kTs,
+                     std::vector<Mgxs>& micros, double_1dvec& atom_densities,
+                     int& method, double tolerance);
+    void combine(std::vector<Mgxs>& micros, double_1dvec& scalars,
+                 int_1dvec& micro_ts, int this_t);
+    void from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups,
+                           double_1dvec temperature, int& method,
+                           double tolerance, int max_order,
+                           bool legendre_to_tabular,
+                           int legendre_to_tabular_points);
+    double get_xs(const char* xstype, int gin, int* gout, double* mu,
+                  int* dg);
+    void sample_fission_energy(int gin, double nu_fission, int& dg, int& gout);
+    void sample_scatter(dir_arr& uvw, int gin, int& gout, double& mu,
+                        double& wgt);
+    void calculate_xs(int gin, double sqrtkT, dir_arr& uvw, double& total_xs,
+                      double& abs_xs, double& nu_fiss_xs);
+    bool equiv(const Mgxs& that);
+    inline void set_temperature_index(double sqrtkT);
+    inline void set_angle_index(dir_arr& uvw);
+};
+
+extern "C" void read_mgxs_library(hid_t file_id, int n_nuclides, char** names,
+     int energy_groups, int delayed_groups, int n_temps, double temps[],
+     int& method, double tolerance, int max_order, bool legendre_to_tabular,
+     int legendre_to_tabular_points);
+extern "C" bool query_fissionable(const int i_nuclides[], const int n_nuclides);
+void create_macro_xs(int n_materials, double_2dvec& mat_kTs,
+     std::vector<std::string>& mat_names, double_1dvec& atom_densities,
+     int& method, double tolerance);
+
+
+// Storage for the MGXS data
+std::vector<Mgxs> nuclides_MG;
+std::vector<Mgxs> macro_xs;
+
+} // namespace openmc
+#endif // MGXS_H

src/scattdata.cpp

@@ -0,0 +1,831 @@
+#include "scattdata.h"
+
+namespace openmc {
+
+//==============================================================================
+// Methods for use by all extended types
+//==============================================================================
+
+
+
+//==============================================================================
+// ScattData base-class methods
+//==============================================================================
+
+void ScattData::generic_init(int order, int_1dvec in_gmin,
+     int_1dvec in_gmax, double_2dvec in_energy, 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++) {
+    // Make sure the energy is normalized
+    double norm = std::accumulate(in_energy[gin].begin(),
+                                  in_energy[gin].end(), 0.);
+
+    if (norm != 0.) {
+      for (auto& n : in_energy[gin]) n /= norm;
+    }
+
+    // Store the inputted data
+    energy[gin] = in_energy[gin];
+    mult[gin] = in_mult[gin];
+
+    // Initialize the distribution data
+    dist[gin].resize(in_gmax[gin] - in_gmin[gin] + 1);
+    for (auto& v : dist[gin]) {
+      v.resize(order);
+      for (auto& n : v) n = 0.;
+    }
+  }
+}
+
+
+void ScattData::sample_energy(int gin, int& gout, int& i_gout)
+{
+  // Sample the outgoing group
+  double xi = prn();
+  i_gout = 0; //TODO: + 1?
+  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(const char* xstype, int gin, int* gout, 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 = 0.;
+  if (std::strcmp(xstype, "scatter")) {
+    if (gout != nullptr) {
+      val = scattxs[gin] * energy[gin][i_gout];
+    } else {
+      val = scattxs[gin];
+    }
+  } else if (std::strcmp(xstype, "scatter/mult")) {
+    if (gout != nullptr) {
+      val = scattxs[gin] * energy[gin][i_gout] / mult[gin][i_gout];
+    } else {
+      val = scattxs[gin] / std::inner_product(mult[gin].begin(),
+                                              mult[gin].end(),
+                                              energy[gin].begin(), 0.0);
+    }
+  } else if (std::strcmp(xstype, "scatter*f_mu/mult")) {
+    if ((gout != nullptr) && (mu != nullptr)) {
+      val = scattxs[gin] * 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");
+    }
+  } else if (std::strcmp(xstype, "scatter*f_mu")) {
+    if ((gout != nullptr) && (mu != nullptr)) {
+      val = scattxs[gin] * 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");
+    }
+  }
+  return val;
+}
+
+
+//==============================================================================
+// ScattDataLegendre methods
+//==============================================================================
+
+void ScattDataLegendre::init(int_1dvec in_gmin, int_1dvec in_gmax,
+                             double_2dvec in_mult, double_3dvec coeffs)
+{
+  int groups = coeffs.size();
+  int order = coeffs[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 = gmax[gin] - gmin[gin] + 1;
+    scattxs[gin] = 0.;
+    for (int i_gout = 0; i_gout < num_groups; i_gout++) {
+      scattxs[gin] = std::accumulate(matrix[gin][i_gout].begin(),
+                                     matrix[gin][i_gout].end(),
+                                     scattxs[gin]);
+    }
+  }
+
+  // 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 = gmax[gin] - 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::generic_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)
+{
+  // TODO: gout >= or gout >?
+  double f;
+  if ((gout < gmin[gin]) || (gout >= gmax[gin])) {
+    f = 0.;
+  } else {
+      // TODO: size() -1 or just size?
+    int i_gout = gout - gmin[gin]; //TODO: + 1?
+    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) {
+    double 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(std::vector<ScattData*> those_scatts,
+                                double_1dvec& scalars)
+{
+  int groups = energy.size();
+  // Find the maximum order in the data set
+  int max_order = get_order();
+  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 (!equiv(*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
+
+  // Now allocate and zero our storage spaces
+  double_3dvec this_matrix = get_matrix(max_order);
+  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++) {
+    ScattDataLegendre* that = dynamic_cast<ScattDataLegendre*>(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
+      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][gout];
+        mult_numer[gin][gout] += scalars[i] * nuscatt;
+        if (that->mult[gin][gout] > 0.) {
+          mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][gout];
+        } else {
+          mult_denom[gin][gout] += scalars[i];
+        }
+      }
+    }
+  }
+
+  // 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.
+  int_1dvec in_gmin(groups);
+  int_1dvec in_gmax(groups);
+  double_3dvec sparse_scatter(groups);
+  double_2dvec sparse_mult(groups);
+  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 = std::all_of(this_matrix[gin][gmin_].begin(),
+                                  this_matrix[gin][gmin_].end(),
+                                  [](double val){return val > 0.;});
+      if (non_zero) break;
+    }
+    int gmax_;
+    for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) {
+      bool non_zero = std::all_of(this_matrix[gin][gmax_].begin(),
+                                  this_matrix[gin][gmax_].end(),
+                                  [](double val){return val > 0.;});
+      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];
+    }
+  }
+
+  // Got everything we need, store it.
+  init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
+}
+
+
+bool ScattDataLegendre::equiv(const ScattDataLegendre& that)
+{
+  // ensure that the number of groups match
+  return (this->energy.size() == that.energy.size());
+}
+
+
+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(order_dim,
+       double_1dvec(order_dim, 0.)));
+
+  for (int gin = 0; gin < groups; gin++) {
+    for (int i_gout = 0; i_gout < energy[0].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(int_1dvec in_gmin, int_1dvec in_gmax,
+                            double_2dvec in_mult, double_3dvec coeffs)
+{
+  int groups = coeffs.size();
+  int order = coeffs[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 = gmax[gin] - 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.);
+      if (norm != 0.) {
+        for (auto& n : matrix[gin][i_gout]) n /= norm;
+      }
+    }
+  }
+
+  // Initialize the base class attributes
+  ScattData::generic_init(order, in_gmin, in_gmax, in_energy,
+                                    in_mult);
+
+  // Build the angular distributio mu values
+  mu = double_1dvec(order);
+  dmu = 2. / order;
+  mu[0] = -1.;
+  for (int imu = 1; imu < order; imu++) {
+    mu[imu] = -1. + (imu - 1) * 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)
+{
+  // TODO: gout >= or gout >?
+  double f;
+  if ((gout < gmin[gin]) || (gout >= gmax[gin])) {
+    f = 0.;
+  } else {
+    // Find mu bin
+    int i_gout = gout - gmin[gin]; //TODO: + 1?
+    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.);
+    }
+
+    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 = 1;
+  } 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];
+}
+
+
+bool ScattDataHistogram::equiv(const ScattDataHistogram& that)
+{
+  bool match = false;
+  if (this->energy.size() == that.energy.size() &&
+      this->dmu == that.dmu &&
+      std::equal(this->mu.begin(), this->mu.end(), that.mu.begin())) {
+    match = true;
+  }
+  return match;
+}
+
+
+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(order_dim,
+       double_1dvec(order_dim, 0.)));
+
+  for (int gin = 0; gin < groups; gin++) {
+    for (int i_gout = 0; i_gout < energy[0].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;
+}
+
+//==============================================================================
+// ScattDataTabular methods
+//==============================================================================
+
+void ScattDataTabular::init(int_1dvec in_gmin, int_1dvec in_gmax,
+                            double_2dvec in_mult, double_3dvec coeffs)
+{
+  int groups = coeffs.size();
+  int order = coeffs[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 - 1) * 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;
+  in_energy.resize(groups);
+  for (int gin = 0; gin < groups; gin++) {
+    int num_groups = gmax[gin] - 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]);
+      }
+      if (norm != 0.) {
+        for (auto& n : matrix[gin][i_gout]) n /= norm;
+      }
+    }
+  }
+
+  // Initialize the base class attributes
+  ScattData::generic_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)
+{
+  // TODO: gout >= or gout >?
+  double f;
+  if ((gout < gmin[gin]) || (gout >= gmax[gin])) {
+    f = 0.;
+  } else {
+    // Find mu bin
+    int i_gout = gout - gmin[gin]; //TODO: + 1?
+    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.);
+    }
+
+    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 - 2; 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 - 1);
+
+  // 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];
+}
+
+
+bool ScattDataTabular::equiv(const ScattDataTabular& that)
+{
+  bool match = false;
+  if (this->energy.size() == that.energy.size() &&
+      this->dmu == that.dmu &&
+      std::equal(this->mu.begin(), this->mu.end(), that.mu.begin())) {
+    match = true;
+  }
+  return match;
+}
+
+
+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(order_dim,
+       double_1dvec(order_dim, 0.)));
+
+  for (int gin = 0; gin < groups; gin++) {
+    for (int i_gout = 0; i_gout < energy[0].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 convert_legendre_to_tabular(ScattDataLegendre& leg,
+                                 ScattDataTabular& tab, int n_mu)
+{
+  // Copy the obvious data
+  tab.energy = leg.energy;
+  tab.mult = leg.mult;
+  tab.gmin = leg.gmin;
+  tab.gmax = leg.gmax;
+
+  // 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 - 1) * 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.;
+      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,106 @@
+//! \file scattdata.h
+//! A collection of multi-group scattering data classes
+
+#ifndef SCATTDATA_H
+#define SCATTDATA_H
+
+#include <algorithm>
+#include <numeric>
+#include <vector>
+#include <cmath>
+
+#include "constants.h"
+#include "math_functions.h"
+#include "random_lcg.h"
+#include "error.h"
+
+namespace openmc {
+
+//==============================================================================
+// SCATTDATA contains all the data needed to describe the scattering energy and
+// angular distribution data
+//==============================================================================
+// temporary declaations so we can name our friend functions
+class ScattDataLegendre;
+class ScattDataTabular;
+
+class ScattData {
+  protected:
+    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
+  public:
+    double_1dvec scattxs; // Isotropic Sigma_{s,g_{in}}
+    virtual double calc_f(int gin, int gout, double mu) = 0;
+    virtual void sample(int gin, int& gout, double& mu, double& wgt) = 0;
+    virtual void init(int_1dvec in_gmin, int_1dvec in_gmax,
+                      double_2dvec in_mult, double_3dvec coeffs) = 0;
+    void sample_energy(int gin, int& gout, int& i_gout);
+    double get_xs(const char* xstype, int gin, int* gout, double* mu);
+    void generic_init(int order, int_1dvec in_gmin, int_1dvec in_gmax,
+                      double_2dvec in_energy, double_2dvec in_mult);
+    virtual void combine(std::vector<ScattData*> those_scatts,
+                         double_1dvec& scalars) = 0;
+    virtual int get_order() = 0;
+    virtual double_3dvec get_matrix(int max_order) = 0;
+};
+
+class ScattDataLegendre: public ScattData {
+  protected:
+    // Maximal value for rejection sampling from a rectangle
+    double_2dvec max_val;
+    friend void convert_legendre_to_tabular(ScattDataLegendre& leg,
+                                            ScattDataTabular& tab, int n_mu);
+  public:
+    void init(int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_mult,
+              double_3dvec coeffs);
+    void update_max_val();
+    double calc_f(int gin, int gout, double mu);
+    void sample(int gin, int& gout, double& mu, double& wgt);
+    bool equiv(const ScattDataLegendre& that);
+    void combine(std::vector<ScattData*> those_scatts, double_1dvec& scalars);
+    int get_order() {return dist[0][0].size() - 1;};
+    double_3dvec get_matrix(int max_order);
+};
+
+class ScattDataHistogram: public ScattData {
+  protected:
+    double_1dvec mu;
+    double dmu;
+    double_3dvec fmu;
+  public:
+    void init(int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_mult,
+              double_3dvec coeffs);
+    double calc_f(int gin, int gout, double mu);
+    void sample(int gin, int& gout, double& mu, double& wgt);
+    void combine(std::vector<ScattData*> those_scatts, double_1dvec& scalars);
+    bool equiv(const ScattDataHistogram& that);
+    int get_order() {return dist[0][0].size();};
+    double_3dvec get_matrix(int max_order);
+};
+
+class ScattDataTabular: public ScattData {
+  protected:
+    double_1dvec mu;
+    double dmu;
+    double_3dvec fmu;
+    friend void convert_legendre_to_tabular(ScattDataLegendre& leg,
+                                            ScattDataTabular& tab, int n_mu);
+  public:
+    void init(int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_mult,
+              double_3dvec coeffs);
+    double calc_f(int gin, int gout, double mu);
+    void sample(int gin, int& gout, double& mu, double& wgt);
+    void combine(std::vector<ScattData*> those_scatts, double_1dvec& scalars);
+    bool equiv(const ScattDataTabular& that);
+    int get_order() {return dist[0][0].size();};
+    double_3dvec get_matrix(int max_order);
+};
+
+void convert_legendre_to_tabular(ScattDataLegendre& leg,
+                                 ScattDataTabular& tab, int n_mu);
+
+} // namespace openmc
+#endif // SCATTDATA_H

src/string_functions.h

@@ -0,0 +1,54 @@
+//! \file string_functions.h
+//! A collection of helper routines for C-strings and STL strings
+
+#ifndef STRING_FUNCTIONS_H
+#define STRING_FUNCTIONS_H
+
+// for string functions
+#include <functional>
+#include <cctype>
+#include <locale>
+#include <string>
+
+namespace openmc {
+
+void strtrim(char* str)
+{
+    int start = 0; // number of leading spaces
+    char* buffer = str;
+
+    while (*str && *str++ == ' ') ++start;
+
+    while (*str++); // move to end of string
+
+    // backup over trailing spaces
+    int end = str - buffer - 1;
+    while (end > 0 && buffer[end - 1] == ' ') --end;
+    buffer[end] = 0; // remove trailing spaces
+
+    // exit if no leading spaces or string is now empty
+    if (end <= start || start == 0) return;
+    str = buffer + start;
+
+    while ((*buffer++ = *str++));  // remove leading spaces: K&R
+}
+
+std::string strtrim(std::string in_str)
+{
+  std::string str = in_str;
+  // perform the left trim
+  str.erase(str.begin(), std::find_if(str.begin(), str.end(),
+            std::not1(std::ptr_fun<int, int>(std::isspace))));
+  // perform the right trim
+  str.erase(std::find_if(str.rbegin(), str.rend(),
+            std::not1(std::ptr_fun<int, int>(std::isspace))).base(),
+                      str.end());
+}
+
+void to_lower(std::string& str)
+{
+  for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
+}
+
+} // namespace openmc
+#endif // STRING_FUNCTIONS_H

src/xsdata.cpp

@@ -0,0 +1,782 @@
+#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)
+{
+  // 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_3dvec(n_pol, double_2dvec(n_azi,
+       double_1dvec(energy_groups, 0.)));
+  absorption = double_3dvec(n_pol, double_2dvec(n_azi,
+       double_1dvec(energy_groups, 0.)));
+  inverse_velocity = double_3dvec(n_pol,
+       double_2dvec(n_azi, double_1dvec(energy_groups, 0.)));
+  if (fissionable) {
+    fission = double_3dvec(n_pol, double_2dvec(n_azi,
+         double_1dvec(energy_groups, 0.)));
+    prompt_nu_fission = double_3dvec(n_pol, double_2dvec(n_azi,
+         double_1dvec(energy_groups, 0.)));
+    kappa_fission = double_3dvec(n_pol, double_2dvec(n_azi,
+         double_1dvec(energy_groups, 0.)));
+  }
+
+  // allocate decay_rate; [temperature][phi][theta][delayed group]
+  decay_rate = double_3dvec(n_pol, double_2dvec(n_azi,
+       double_1dvec(num_delayed_groups, 0.)));
+
+  if (fissionable) {
+    // allocate delayed_nu_fission; [temperature][phi][theta][in group][out group]
+    delayed_nu_fission = double_4dvec(n_pol, double_3dvec(n_azi,
+         double_2dvec(energy_groups, double_1dvec(energy_groups, 0.))));
+
+    // chi_prompt; [temperature][phi][theta][in group][delayed group]
+    chi_prompt = double_4dvec(n_pol, double_3dvec(n_azi,
+         double_2dvec(energy_groups, double_1dvec(num_delayed_groups, 0.))));
+
+    // chi_delayed; [temperature][phi][theta][in group][out group][delay group]
+    chi_delayed = double_5dvec(n_pol, double_4dvec(n_azi,
+         double_3dvec(energy_groups, double_2dvec(energy_groups,
+         double_1dvec(num_delayed_groups, 0.)))));
+  }
+
+  scatter.resize(n_pol);
+  for (int p = 0; p < n_pol; p++) {
+    scatter[p].resize(n_azi);
+    for (int a = 0; a < n_azi; a++) {
+      if (scatter_format == ANGLE_HISTOGRAM) {
+        scatter[p][a] = new ScattDataHistogram;
+      } else if (scatter_format == ANGLE_TABULAR) {
+        scatter[p][a] = new ScattDataTabular;
+      } else if (scatter_format == ANGLE_LEGENDRE) {
+        scatter[p][a] = new ScattDataLegendre;
+      }
+    }
+  }
+}
+
+XsData::~XsData()
+{
+  for (int p = 0; p < scatter.size(); p++) {
+    for (int a = 0; a < scatter[p].size(); a++) delete scatter[p][a];
+    scatter[p].clear();
+  }
+  scatter.clear();
+}
+
+
+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)
+{
+  // Reconstruct the dimension information so it doesn't need to be passed
+  int n_pol = total.size();
+  int n_azi = total[0].size();
+  int energy_groups = total[0][0].size();
+  int delayed_groups = decay_rate[0][0].size();
+
+  // Set the fissionable-specific data
+  if (fissionable) {
+    _fissionable_from_hdf5(xsdata_grp, n_pol, n_azi, energy_groups,
+                           delayed_groups);
+  }
+  // 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 p = 0; p < n_pol; p++) {
+    for (int a = 0; a < n_azi; a++) {
+      for (int gin = 0; gin < energy_groups; gin++) {
+        if (absorption[gin][p][a] == 0.) {
+          absorption[p][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 p = 0; p < n_pol; p++) {
+      for (int a = 0; a < n_azi; a++) {
+        for (int gin = 0; gin < energy_groups; gin++) {
+          total[p][a][gin] = absorption[p][a][gin] +
+               scatter[p][a]->scattxs[gin];
+        }
+      }
+    }
+  }
+
+  // Check total to ensure it is not 0 since it is often the denominator in
+  // tally methods
+  for (int p = 0; p < n_pol; p++) {
+    for (int a = 0; a < n_azi; a++) {
+      for (int gin = 0; gin < energy_groups; gin++) {
+        if (total[p][a][gin] == 0.) total[p][a][gin] = 1.e-10;
+      }
+    }
+  }
+}
+
+
+void XsData::_fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
+     int energy_groups, int delayed_groups)
+{
+  double_4dvec temp_beta =
+       double_4dvec(n_pol, double_3dvec(n_azi,
+       double_2dvec(energy_groups, double_1dvec(delayed_groups, 0.))));
+
+  // Set/get beta
+  if (object_exists(xsdata_grp, "beta")) {
+    hid_t xsdata = open_dataset(xsdata_grp, "beta");
+    int ndims = dataset_ndims(xsdata);
+
+    if (ndims == 3) {
+      // Beta is input as [delayed group]
+      double_1dvec temp_arr = double_1dvec(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 p = 0; p < n_pol; p++) {
+        for (int a = 0; a < n_azi; a++) {
+          for (int dg = 0; dg < delayed_groups; dg++) {
+            // Set the first group index and copy the rest
+            temp_beta[p][a][0][dg] = temp_arr[temp_idx++];
+            for (int gin = 1; gin < energy_groups; gin++) {
+              temp_beta[p][a][gin] = temp_beta[p][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 1D or 2D array!");
+    }
+  }
+
+  // If chi is provided, set chi-prompt and chi-delayed
+  if (object_exists(xsdata_grp, "chi")) {
+    double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi,
+         double_1dvec(energy_groups)));
+    read_nd_vector(xsdata_grp, "chi", temp_arr);
+
+    int temp_idx = 0;
+    for (int p = 0; p < n_pol; p++) {
+      for (int a = 0; a < n_azi; a++) {
+        // First set the first group
+        for (int gout = 0; gout < energy_groups; gout++) {
+          chi_prompt[p][a][0][gout] = temp_arr[p][a][temp_idx++];
+        }
+
+        // Now normalize this data
+        double chi_sum = std::accumulate(chi_prompt[p][a][0].begin(),
+                                         chi_prompt[p][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[p][a][0][gout] /= chi_sum;
+        }
+
+        // And extend to the remaining incoming groups
+        for (int gin = 1; gin < energy_groups; gin++) {
+          chi_prompt[p][a][gin] = chi_prompt[p][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[p][a][gin][gout][dg] =
+                   chi_prompt[p][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);
+
+    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 p = 0; p < n_pol; p++) {
+        for (int a = 0; a < n_azi; a++) {
+          for (int gin = 0; gin < energy_groups; gin++) {
+            for (int dg = 0; dg < delayed_groups; dg++) {
+              delayed_nu_fission[p][a][gin][dg] =
+                   temp_beta[p][a][gin][dg] *
+                   prompt_nu_fission[p][a][gin];
+            }
+
+            // Correct the prompt-nu-fission using the delayed neutron fraction
+            if (delayed_groups > 0) {
+              double beta_sum = std::accumulate(temp_beta[p][a][gin].begin(),
+                                                temp_beta[p][a][gin].end(), 0.);
+              prompt_nu_fission[p][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 p = 0; p < n_pol; p++) {
+        for (int a = 0; a < n_azi; a++) {
+          for (int gin = 0; gin < energy_groups; gin++) {
+            double chi_sum = std::accumulate(chi_prompt[p][a][gin].begin(),
+                                             chi_prompt[p][a][gin].end(), 0.);
+            if (chi_sum >= 0.) {
+              for (int gout = 0; gout < energy_groups; gout++) {
+                chi_prompt[p][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[p][a][gin][gout][dg] =
+                     chi_prompt[p][a][gin][gout];
+              }
+            }
+          }
+
+          // Set the vector nu-fission from the matrix nu-fission
+          for (int gin = 0; gin < energy_groups; gin++) {
+            double sum = std::accumulate(chi_prompt[p][a][gin].begin(),
+                                         chi_prompt[p][a][gin].end(), 0.);
+            prompt_nu_fission[p][a][gin] = sum;
+          }
+
+          // 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[p][a][gin][dg] =
+                   temp_beta[p][a][gin][dg] *
+                   prompt_nu_fission[p][a][gin];
+            }
+
+            // Correct prompt-nu-fission using the delayed neutron fraction
+            if (delayed_groups > 0) {
+              double beta_sum = std::accumulate(temp_beta[p][a][gin].begin(),
+                                                temp_beta[p][a][gin].end(), 0.);
+              prompt_nu_fission[p][a][gin] *= (1. - beta_sum);
+            }
+          }
+        }
+      }
+    } else {
+      fatal_error("beta 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_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi,
+         double_1dvec(energy_groups)));
+    read_nd_vector(xsdata_grp, "chi-prompt", temp_arr);
+
+    for (int a = 0; a < n_azi; a++) {
+      for (int p = 0; p < n_pol; p++) {
+        for (int gin = 0; gin < energy_groups; gin++) {
+          for (int gout = 0; gout < energy_groups; gout++) {
+            chi_prompt[p][a][gin][gout] = temp_arr[p][a][gout];
+          }
+
+          // Normalize chi so its CDF goes to 1
+          double chi_sum = std::accumulate(chi_prompt[p][a][gin].begin(),
+                                           chi_prompt[p][a][gin].end(), 0.);
+          if (chi_sum >= 0.) {
+            for (int gout = 0; gout < energy_groups; gout++) {
+              chi_prompt[p][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);
+    close_dataset(xsdata);
+
+    if (ndims == 3) {
+      // chi-delayed is a [in group] vector
+      double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi,
+           double_1dvec(energy_groups)));
+      read_nd_vector(xsdata_grp, "chi-delayed", temp_arr);
+
+      for (int a = 0; a < n_azi; a++) {
+        for (int p = 0; p < n_pol; p++) {
+          // normalize the chi CDF to 1
+          double chi_sum = std::accumulate(temp_arr[p][a].begin(),
+                                           temp_arr[p][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[p][a][gin][gout][dg] =
+                     temp_arr[p][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_azi; a++) {
+        for (int p = 0; p < n_pol; p++) {
+          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[p][a][gin][gout][dg];
+              }
+
+              if (chi_sum > 0.) {
+                for (int gout = 0; gout < energy_groups; gout++) {
+                  chi_delayed[p][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);
+    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_4dvec temp_arr = double_4dvec(n_pol, double_3dvec(n_azi,
+           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_azi; a++) {
+        for (int p = 0; p < n_pol; p++) {
+          for (int gin = 0; gin < energy_groups; gin++) {
+            double prompt_sum = std::accumulate(temp_arr[p][a][gin].begin(),
+                                                temp_arr[p][a][gin].end(), 0.);
+            prompt_nu_fission[p][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[p][a][gin] > 0.) {
+              for (int gout = 0; gout < energy_groups; gout++) {
+                chi_prompt[p][a][gin][gout] =
+                     temp_arr[p][a][gin][gout] /
+                     prompt_nu_fission[p][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);
+
+    if (ndims == 3) {
+      // delayed-nu-fission is a [in group] vector
+      if (temp_beta[0][0][0][0] == 0.) {
+        fatal_error("cannot set delayed-nu-fission with a 1D array if "
+                    "beta is not provided");
+      }
+      double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi,
+         double_1dvec(energy_groups)));
+      read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
+
+      for (int p = 0; p < n_pol; p++) {
+        for (int a = 0; a < n_azi; 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[p][a][gin][dg] =
+                   temp_beta[p][a][gin][dg] * temp_arr[p][a][gin];
+            }
+          }
+        }
+      }
+
+    } else if (ndims == 4) {
+      // delayed nu fission is a [pol][azi][energy_group][delayed_group] matrix;
+      // matrix use this to set delayed-nu-fission separately for each
+      // delayed group
+      std::vector<hsize_t> dims(ndims);
+      get_shape(xsdata, &dims[0]);
+
+      if (dims[2] != delayed_groups) {
+        fatal_error("The delayed-nu-fission matrix was input with a 1st "
+                    "dimension not equal to the number of delayed groups");
+      }
+      if (dims[3] != energy_groups) {
+        fatal_error("The delayed-nu-fission matrix was input with a 2nd "
+                    "dimension not equal to the number of energy groups");
+      }
+      if (delayed_groups == energy_groups) {
+        warning("delayed-nu-fission was input as a dimension-4 matrix "
+                "with the same number of delayed groups and energy "
+                "groups. OpenMC assumes the dimensions in the matrix "
+                "are [delayed_groups][energy_groups].  Currently, "
+                "delayed-nu-fission cannot be set as a group-by-group "
+                "matrix");
+      }
+
+      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_5dvec temp_arr = double_5dvec(n_pol, double_4dvec(n_azi,
+         double_3dvec(energy_groups, double_2dvec(energy_groups,
+         double_1dvec(delayed_groups)))));
+      read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
+
+      // Set the 4D delayed-nu-fission matrix and 5D chi-delayed matrix
+      // from the 5D delayed-nu-fission matrix
+      for (int p = 0; p < n_pol; p++) {
+        for (int a = 0; a < n_azi; 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[p][a][gin][gout][dg];
+                chi_delayed[p][a][gin][gout][dg] =
+                     temp_arr[p][a][gin][gout][dg];
+              }
+              delayed_nu_fission[p][a][gin][dg] = gout_sum;
+              // Normalize chi-delayed
+              if (gout_sum > 0.) {
+                for (int gout = 0; gout < energy_groups; gout++) {
+                  chi_delayed[p][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!");
+    }
+    close_dataset(xsdata);
+  }
+
+  // Get the fission and kappa_fission data xs
+  read_nd_vector(xsdata_grp, "fission", fission);
+  read_nd_vector(xsdata_grp, "kappa-fission", kappa_fission);
+}
+
+
+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)
+{
+  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_3dvec gmin = int_3dvec(n_pol, int_2dvec(n_azi,
+       int_1dvec(energy_groups)));
+  read_nd_vector(scatt_grp, "g_min", gmin, true);
+  int_3dvec gmax = int_3dvec(n_pol, int_2dvec(n_azi,
+       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 p = 0; p < n_pol; p++) {
+    for (int a = 0; a < n_azi; a++) {
+      for (int gin = 0; gin < energy_groups; gin++) {
+        gmin[p][a][gin] -= 1;
+        gmax[p][a][gin] -= 1;
+      }
+    }
+  }
+
+  // Now use this info to find the length of a vector to hold the flattened
+  // data.
+  int length = 0;
+  for (int p = 0; p < n_pol; p++) {
+    for (int a = 0; a < n_azi; a++) {
+      for (int gin = 0; gin < energy_groups; gin++) {
+        length += order_data * (gmax[p][a][gin] - gmin[p][a][gin] + 1);
+      }
+    }
+  }
+  double_1dvec temp_arr = double_1dvec(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_5dvec input_scatt =
+       double_5dvec(n_pol, double_4dvec(n_azi, double_3dvec(energy_groups)));
+
+  int temp_idx = 0;
+  for (int p = 0; p < n_pol; p++) {
+    for (int a = 0; a < n_azi; a++) {
+      for (int gin = 0; gin < energy_groups; gin++) {
+        input_scatt[p][a][gin].resize(gmax[p][a][gin] - gmin[p][a][gin] + 1);
+        for (int i_gout = 0; i_gout < input_scatt[p][a][gin].size(); i_gout++) {
+          input_scatt[p][a][gin][i_gout].resize(order_dim);
+          for (int l = 0; l < order_dim; l++) {
+            input_scatt[p][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_4dvec temp_mult = double_4dvec(n_pol, double_3dvec(n_azi,
+       double_2dvec(energy_groups)));
+  if (object_exists(scatt_grp, "multiplicity_matrix")) {
+    temp_arr.resize(length);
+    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 p = 0; p < n_pol; p++) {
+      for (int a = 0; a < n_azi; a++) {
+        for (int gin = 0; gin < energy_groups; gin++) {
+          temp_mult[p][a][gin].resize(gmax[p][a][gin] - gmin[p][a][gin] + 1);
+          for (int i_gout = 0; i_gout < temp_mult[p][a][gin].size(); i_gout++) {
+            temp_mult[p][a][gin][i_gout] = temp_arr[temp_idx++];
+          }
+        }
+      }
+    }
+  } else {
+    // Use a default: multiplicities are 1.0.
+    for (int p = 0; p < n_pol; p++) {
+      for (int a = 0; a < n_azi; a++) {
+        for (int gin = 0; gin < energy_groups; gin++) {
+          temp_mult[p][a][gin].resize(gmax[p][a][gin] - gmin[p][a][gin] + 1);
+          for (int i_gout = 0; i_gout < temp_mult[p][a][gin].size(); i_gout++) {
+            temp_mult[p][a][gin][i_gout] = 1.;
+          }
+        }
+      }
+    }
+  }
+  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 p = 0; p < n_pol; p++) {
+      for (int a = 0; a < n_azi; a++) {
+        ScattDataLegendre legendre_scatt;
+        legendre_scatt.init(gmin[p][a], gmax[p][a], temp_mult[p][a],
+                            input_scatt[p][a]);
+
+        // Now create a tabular version of legendre_scatt
+        convert_legendre_to_tabular(legendre_scatt,
+             *static_cast<ScattDataTabular*>(scatter[p][a]),
+             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 p = 0; p < n_pol; p++) {
+      for (int a = 0; a < n_azi; a++) {
+        scatter[p][a]->init(gmin[p][a], gmax[p][a], temp_mult[p][a],
+                                  input_scatt[p][a]);
+      }
+    }
+  }
+}
+
+
+void XsData::combine(std::vector<XsData*> those_xs, 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 p = 0; p < total.size(); p++) {
+      for (int a = 0; a < total[p].size(); a++) {
+        for (int gin = 0; gin < total[p][a].size(); gin++) {
+          total[p][a][gin] += scalar * that->total[p][a][gin];
+          absorption[p][a][gin] += scalar * that->absorption[p][a][gin];
+          inverse_velocity[p][a][gin] +=
+                 scalar * that->inverse_velocity[p][a][gin];
+
+          prompt_nu_fission[p][a][gin] +=
+               scalar * that->prompt_nu_fission[p][a][gin];
+          kappa_fission[p][a][gin] +=
+               scalar * that->kappa_fission[p][a][gin];
+          fission[p][a][gin] +=
+               scalar * that->fission[p][a][gin];
+
+          for (int dg = 0; dg < delayed_nu_fission[p][a][gin].size(); dg++) {
+            delayed_nu_fission[p][a][gin][dg] +=
+                 scalar * that->delayed_nu_fission[p][a][gin][dg];
+          }
+
+          for (int gout = 0; gout < chi_prompt[p][a][gin].size(); gout++) {
+            chi_prompt[p][a][gin][gout] +=
+                 scalar * that->chi_prompt[p][a][gin][gout];
+
+            for (int dg = 0; dg < chi_delayed[p][a][gin][gout].size(); dg++) {
+              chi_delayed[p][a][gin][gout][dg] +=
+                   scalar * that->chi_delayed[p][a][gin][gout][dg];
+            }
+          }
+        }
+
+        for (int dg = 0; dg < decay_rate[p][a].size(); dg++) {
+          decay_rate[p][a][dg] += scalar * that->decay_rate[p][a][dg];
+        }
+
+        // Normalize chi
+        for (int gin = 0; gin < chi_prompt[p][a].size(); gin++) {
+          double norm = std::accumulate(chi_prompt[p][a][gin].begin(),
+                                        chi_prompt[p][a][gin].end(), 0.);
+          if (norm > 0.) {
+            for (int gout = 0; gout < chi_prompt[p][a][gin].size(); gout++) {
+              chi_prompt[p][a][gin][gout] /= norm;
+            }
+          }
+
+          for (int dg = 0; dg < chi_delayed[p][a][gin][0].size(); dg++) {
+            norm = 0.;
+            for (int gout = 0; gout < chi_delayed[p][a][gin].size(); gout++) {
+              norm += chi_delayed[p][a][gin][gout][dg];
+            }
+            if (norm > 0.) {
+              for (int gout = 0; gout < chi_delayed[p][a][gin].size(); gout++) {
+                chi_delayed[p][a][gin][gout][dg] /= norm;
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  // Allow the ScattData object to combine itself
+  for (int p = 0; p < total.size(); p++) {
+    for (int a = 0; a < total[p].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[p][a];
+      }
+
+      // Now combine these guys
+      scatter[p][a]->combine(those_scatts, scalars);
+    }
+  }
+}
+
+
+bool XsData::equiv(const XsData& that)
+{
+  bool match = false;
+  // check n_pol (total.size()), n_azi (total[0].size()), and
+  // groups (total[0][0].size())
+  // This assumes correct initializatino of the remaining cross sections
+  if ((total.size() == that.total.size()) &&
+      (total[0].size() == that.total[0].size()) &&
+      (total[0][0].size() == that.total[0][0].size())) {
+    match = true;
+  }
+  return match;
+}
+
+} //namespace openmc

src/xsdata.h

@@ -0,0 +1,72 @@
+//! \file xsdata.h
+//! A collection of classes for containing the Multi-Group Cross Section data
+
+#ifndef XSDATA_H
+#define XSDATA_H
+
+#include <cmath>
+#include <cstdlib>
+#include <algorithm>
+#include <memory>
+#include <string>
+#include <valarray>
+#include <vector>
+
+#include "constants.h"
+#include "hdf5_interface.h"
+#include "math_functions.h"
+#include "random_lcg.h"
+#include "scattdata.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// XSDATA contains the temperature-independent cross section data for an MGXS
+//==============================================================================
+
+class XsData {
+  private:
+    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);
+    void _fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
+         int energy_groups, int delayed_groups);
+  public:
+    // The following quantities have the following dimensions:
+    // [phi][theta][incoming group]
+    double_3dvec total;
+    double_3dvec absorption;
+    double_3dvec prompt_nu_fission;
+    double_3dvec kappa_fission;
+    double_3dvec fission;
+    double_3dvec inverse_velocity;
+    // decay_rate has the following dimensions:
+    // [phi][theta][delayed group]
+    double_3dvec decay_rate;
+    // delayed_nu_fission has the following dimensions:
+    // [phi][theta][incoming group][delayed group]
+    double_4dvec delayed_nu_fission;
+    // chi_prompt has the following dimensions:
+    // [phi][theta][incoming group][outgoing group]
+    double_4dvec chi_prompt;
+    // chi_delayed has the following dimensions:
+    // [phi][theta][incoming group][outgoing group][delayed group]
+    double_5dvec chi_delayed;
+    // scatter has the following dimensions: [phi][theta]
+    std::vector<std::vector<ScattData*> > scatter;
+
+    XsData() = default;
+    XsData(int num_groups, int num_delayed_groups, bool fissionable,
+           int scatter_format, int n_pol, int n_azi);
+    ~XsData();
+    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);
+    void combine(std::vector<XsData*> those_xs, double_1dvec& scalars);
+    bool equiv(const XsData& that);
+};
+
+
+} //namespace openmc
+#endif // XSDATA_H