diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp
index 39e0162b5..a81f26e23 100644
--- a/src/hdf5_interface.cpp
+++ b/src/hdf5_interface.cpp
@@ -481,6 +481,28 @@ read_nd_vector(hid_t obj_id, const char* 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();
+    std::vector<int> temp_arr = std::vector<int>(dim1 * dim2);
+    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++) {
+        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,
diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h
index 59a914f52..b00fe7c87 100644
--- a/src/hdf5_interface.h
+++ b/src/hdf5_interface.h
@@ -65,6 +65,10 @@ 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,
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index ab2970b2f..81a510e30 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -34,11 +34,10 @@ Mgxs::init(const std::string& in_name, const double in_awr,
   for (int thread = 0; thread < n_threads; thread++) {
     cache[thread].sqrtkT = 0.;
     cache[thread].t = 0;
-    cache[thread].p = 0;
     cache[thread].a = 0;
-    cache[thread].uvw[0] = 1.;
-    cache[thread].uvw[1] = 0.;
-    cache[thread].uvw[2] = 0.;
+    cache[thread].u = 0.;
+    cache[thread].v = 0.;
+    cache[thread].w = 0.;
   }
 }
 
@@ -48,7 +47,7 @@ void
 Mgxs::_metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
      const int in_num_delayed_groups, double_1dvec& temperature, int& method,
      const double tolerance, int_1dvec& temps_to_read, int& order_dim,
-     bool& is_isotropic, const int n_threads)
+     const int n_threads)
 {
   // get name
   char char_name[MAX_WORD_LEN];
@@ -179,7 +178,7 @@ Mgxs::_metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
       fatal_error("Invalid scatter_shape option!");
     }
   }
-  //TODO: do i even need this flag? - it should be easy to self-determine
+
   bool in_fissionable = false;
   if (attribute_exists(xs_id, "fissionable")) {
     int int_fiss;
@@ -264,9 +263,8 @@ Mgxs::from_hdf5(hid_t xs_id, const int energy_groups,
   // Call generic data gathering routine (will populate the metadata)
   int order_data;
   int_1dvec temps_to_read;
-  bool is_isotropic;
   _metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature,
-       method, tolerance, temps_to_read, order_data, is_isotropic, n_threads);
+       method, tolerance, temps_to_read, order_data, n_threads);
 
   // Set number of energy and delayed groups
   int final_scatter_format = scatter_format;
@@ -284,7 +282,7 @@ Mgxs::from_hdf5(hid_t xs_id, const int energy_groups,
 
     xs[t].from_hdf5(xsdata_grp, fissionable, scatter_format,
                     final_scatter_format, order_data, max_order,
-                    legendre_to_tabular_points, is_isotropic);
+                    legendre_to_tabular_points, is_isotropic, n_pol, n_azi);
     close_group(xsdata_grp);
 
   } // end temperature loop
@@ -419,45 +417,41 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
   double val;
   switch(xstype) {
   case MG_GET_XS_TOTAL:
-    val = xs[cache[tid].t].total[cache[tid].p][cache[tid].a][gin];
+    val = xs[cache[tid].t].total[cache[tid].a][gin];
     break;
-  case MG_GET_XS_ABSORPTION:
-    val = xs[cache[tid].t].absorption[cache[tid].p][cache[tid].a][gin];
-    break;
-  case MG_GET_XS_INVERSE_VELOCITY:
-    val = xs[cache[tid].t].inverse_velocity[cache[tid].p][cache[tid].a][gin];
-    break;
-  case MG_GET_XS_DECAY_RATE:
-    if (dg != nullptr) {
-      val = xs[cache[tid].t].decay_rate[cache[tid].p][cache[tid].a][*dg + 1];
+  case MG_GET_XS_NU_FISSION:
+    if (fissionable) {
+      val = xs[cache[tid].t].nu_fission[cache[tid].a][gin];
     } else {
-      val = xs[cache[tid].t].decay_rate[cache[tid].p][cache[tid].a][0];
+      val = 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[cache[tid].t].scatter[cache[tid].p]
-                    [cache[tid].a]->get_xs(xstype, gin, gout, mu);
+  case MG_GET_XS_ABSORPTION:
+    val = xs[cache[tid].t].absorption[cache[tid].a][gin];
     break;
   case MG_GET_XS_FISSION:
     if (fissionable) {
-      val = xs[cache[tid].t].fission[cache[tid].p][cache[tid].a][gin];
+      val = xs[cache[tid].t].fission[cache[tid].a][gin];
     } else {
       val = 0.;
     }
     break;
   case MG_GET_XS_KAPPA_FISSION:
     if (fissionable) {
-      val = xs[cache[tid].t].kappa_fission[cache[tid].p][cache[tid].a][gin];
+      val = xs[cache[tid].t].kappa_fission[cache[tid].a][gin];
     } else {
       val = 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[cache[tid].t].scatter[cache[tid].a]->get_xs(xstype, gin, gout, mu);
+    break;
   case MG_GET_XS_PROMPT_NU_FISSION:
     if (fissionable) {
-      val = xs[cache[tid].t].prompt_nu_fission[cache[tid].p][cache[tid].a][gin];
+      val = xs[cache[tid].t].prompt_nu_fission[cache[tid].a][gin];
     } else {
       val = 0.;
     }
@@ -465,10 +459,10 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
   case MG_GET_XS_DELAYED_NU_FISSION:
     if (fissionable) {
       if (dg != nullptr) {
-        val = xs[cache[tid].t].delayed_nu_fission[cache[tid].p][cache[tid].a][gin][*dg];
+        val = xs[cache[tid].t].delayed_nu_fission[cache[tid].a][gin][*dg];
       } else {
         val = 0.;
-        for (auto& num : xs[cache[tid].t].delayed_nu_fission[cache[tid].p]
+        for (auto& num : xs[cache[tid].t].delayed_nu_fission
              [cache[tid].a][gin]) {
           val += num;
         }
@@ -477,21 +471,14 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
       val = 0.;
     }
     break;
-  case MG_GET_XS_NU_FISSION:
-    if (fissionable) {
-      val = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
-    } else {
-      val = 0.;
-    }
-    break;
   case MG_GET_XS_CHI_PROMPT:
     if (fissionable) {
       if (gout != nullptr) {
-        val = xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][*gout];
+        val = xs[cache[tid].t].chi_prompt[cache[tid].a][gin][*gout];
       } else {
         // provide an outgoing group-wise sum
         val = 0.;
-        for (auto& num : xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin]) {
+        for (auto& num : xs[cache[tid].t].chi_prompt[cache[tid].a][gin]) {
           val += num;
         }
       }
@@ -503,22 +490,22 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
     if (fissionable) {
       if (gout != nullptr) {
         if (dg != nullptr) {
-          val = xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][*gout][*dg];
+          val = xs[cache[tid].t].chi_delayed[cache[tid].a][gin][*gout][*dg];
         } else {
-          val = xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][*gout][0];
+          val = xs[cache[tid].t].chi_delayed[cache[tid].a][gin][*gout][0];
         }
       } else {
         if (dg != nullptr) {
           val = 0.;
-          for (int i = 0; i < xs[cache[tid].t].chi_delayed[cache[tid].p]
+          for (int i = 0; i < xs[cache[tid].t].chi_delayed
                [cache[tid].a][gin].size(); i++) {
-            val += xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][i][*dg];
+            val += xs[cache[tid].t].chi_delayed[cache[tid].a][gin][i][*dg];
           }
         } else {
           val = 0.;
-          for (int i = 0; i < xs[cache[tid].t].chi_delayed[cache[tid].p]
+          for (int i = 0; i < xs[cache[tid].t].chi_delayed
                                             [cache[tid].a][gin].size(); i++) {
-            for (auto& num : xs[cache[tid].t].chi_delayed[cache[tid].p]
+            for (auto& num : xs[cache[tid].t].chi_delayed
                                            [cache[tid].a][gin][i]) {
               val += num;
             }
@@ -529,6 +516,16 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
       val = 0.;
     }
     break;
+  case MG_GET_XS_INVERSE_VELOCITY:
+    val = xs[cache[tid].t].inverse_velocity[cache[tid].a][gin];
+    break;
+  case MG_GET_XS_DECAY_RATE:
+    if (dg != nullptr) {
+      val = xs[cache[tid].t].decay_rate[cache[tid].a][*dg + 1];
+    } else {
+      val = xs[cache[tid].t].decay_rate[cache[tid].a][0];
+    }
+    break;
   default:
     val = 0.;
   }
@@ -541,11 +538,11 @@ void
 Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
 {
   // This method assumes that the temperature and angle indices are set
-  double nu_fission = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
+  double nu_fission = xs[cache[tid].t].nu_fission[cache[tid].a][gin];
 
   // Find the probability of having a prompt neutron
   double prob_prompt =
-       xs[cache[tid].t].prompt_nu_fission[cache[tid].p][cache[tid].a][gin];
+       xs[cache[tid].t].prompt_nu_fission[cache[tid].a][gin];
 
   // sample random numbers
   double xi_pd = prn() * nu_fission;
@@ -561,10 +558,10 @@ Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
     // sample the outgoing energy group
     gout = 0;
     double prob_gout =
-         xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][gout];
+         xs[cache[tid].t].chi_prompt[cache[tid].a][gin][gout];
     while (prob_gout < xi_gout) {
       gout++;
-      prob_gout += xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][gout];
+      prob_gout += xs[cache[tid].t].chi_prompt[cache[tid].a][gin][gout];
     }
 
   } else {
@@ -575,7 +572,7 @@ Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
     while (xi_pd >= prob_prompt) {
       dg++;
       prob_prompt +=
-           xs[cache[tid].t].delayed_nu_fission[cache[tid].p][cache[tid].a][gin][dg];
+           xs[cache[tid].t].delayed_nu_fission[cache[tid].a][gin][dg];
     }
 
     // adjust dg in case of round-off error
@@ -584,11 +581,11 @@ Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
     // sample the outgoing energy group
     gout = 0;
     double prob_gout =
-         xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][gout][dg];
+         xs[cache[tid].t].chi_delayed[cache[tid].a][gin][gout][dg];
     while (prob_gout < xi_gout) {
       gout++;
       prob_gout +=
-           xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][gout][dg];
+           xs[cache[tid].t].chi_delayed[cache[tid].a][gin][gout][dg];
     }
   }
 }
@@ -601,7 +598,7 @@ Mgxs::sample_scatter(const int tid, const int gin, int& gout, double& mu,
 {
   // This method assumes that the temperature and angle indices are set
   // Sample the data
-  xs[cache[tid].t].scatter[cache[tid].p][cache[tid].a]->sample(gin, gout, mu, wgt);
+  xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt);
 }
 
 //==============================================================================
@@ -613,11 +610,11 @@ Mgxs::calculate_xs(const int tid, const int gin, const double sqrtkT,
   // Set our indices
   set_temperature_index(tid, sqrtkT);
   set_angle_index(tid, uvw);
-  total_xs = xs[cache[tid].t].total[cache[tid].p][cache[tid].a][gin];
-  abs_xs = xs[cache[tid].t].absorption[cache[tid].p][cache[tid].a][gin];
+  total_xs = xs[cache[tid].t].total[cache[tid].a][gin];
+  abs_xs = xs[cache[tid].t].absorption[cache[tid].a][gin];
 
   if (fissionable) {
-    nu_fiss_xs = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
+    nu_fiss_xs = xs[cache[tid].t].nu_fission[cache[tid].a][gin];
   } else {
     nu_fiss_xs = 0.;
   }
@@ -670,22 +667,25 @@ void
 Mgxs::set_angle_index(const int tid, const double uvw[3])
 {
   // See if we need to find the new index
-  if ((uvw[0] != cache[tid].uvw[0]) || (uvw[1] != cache[tid].uvw[1]) ||
-      (uvw[2] != cache[tid].uvw[2])) {
+  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;
-    cache[tid].p = std::floor(my_pol / delta_angle);
+    int p = std::floor(my_pol / delta_angle);
     delta_angle = 2. * PI / n_azi;
-    cache[tid].a = std::floor((my_azi + PI) / delta_angle);
+    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].uvw[0] = uvw[0];
-    cache[tid].uvw[1] = uvw[1];
-    cache[tid].uvw[2] = uvw[2];
+    cache[tid].u = uvw[0];
+    cache[tid].v = uvw[1];
+    cache[tid].w = uvw[2];
   }
 }
 
diff --git a/src/mgxs.h b/src/mgxs.h
index d7b9a814a..2a79bfe2a 100644
--- a/src/mgxs.h
+++ b/src/mgxs.h
@@ -29,9 +29,11 @@ namespace openmc {
 struct CacheData {
   double sqrtkT; // last temperature corresponding to t
   int t; // temperature index
-  int p; // polar angle index
-  int a; // azimuthal angle index
-  double uvw[3]; // last angle that corresponds to p and a
+  int a; // angle index
+  // last angle that corresponds to p and a
+  double u;
+  double v;
+  double w;
 };
 
 //==============================================================================
@@ -45,6 +47,8 @@ class Mgxs {
     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;
@@ -52,7 +56,7 @@ class Mgxs {
     void _metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
          const int in_num_delayed_groups, double_1dvec& temperature,
          int& method, const double tolerance, int_1dvec& temps_to_read,
-         int& order_dim, bool& is_isotropic, const int n_threads);
+         int& order_dim, const int n_threads);
     bool equiv(const Mgxs& that);
 
   public:
diff --git a/src/scattdata.cpp b/src/scattdata.cpp
index c0515f1ba..e321bc275 100644
--- a/src/scattdata.cpp
+++ b/src/scattdata.cpp
@@ -876,9 +876,9 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
 {
   // 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 P1 (super fast), or the default if
-    // a higher order
-    if (leg.get_order() <= 1) {
+    // 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;
diff --git a/src/xsdata.cpp b/src/xsdata.cpp
index 991e2a901..05e6ab01b 100644
--- a/src/xsdata.cpp
+++ b/src/xsdata.cpp
@@ -10,59 +10,49 @@ XsData::XsData(const int energy_groups, const int num_delayed_groups,
                const bool fissionable, const int scatter_format,
                const int n_pol, const 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_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.)));
+  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_3dvec(n_pol, double_2dvec(n_azi,
-         double_1dvec(energy_groups, 0.)));
-    nu_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.)));
+    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_3dvec(n_pol, double_2dvec(n_azi,
-       double_1dvec(num_delayed_groups, 0.)));
+  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_4dvec(n_pol, double_3dvec(n_azi,
-         double_2dvec(energy_groups, double_1dvec(num_delayed_groups, 0.))));
+    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_4dvec(n_pol, double_3dvec(n_azi,
-         double_2dvec(energy_groups, double_1dvec(energy_groups, 0.))));
+    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_5dvec(n_pol, double_4dvec(n_azi,
-         double_3dvec(energy_groups, double_2dvec(energy_groups,
-         double_1dvec(num_delayed_groups, 0.)))));
+    chi_delayed = double_4dvec(n_ang, 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;
-      }
+  scatter.resize(n_ang);
+  for (int a = 0; a < n_ang; a++) {
+    if (scatter_format == ANGLE_HISTOGRAM) {
+      scatter[a] = new ScattDataHistogram;
+    } else if (scatter_format == ANGLE_TABULAR) {
+      scatter[a] = new ScattDataTabular;
+    } else if (scatter_format == ANGLE_LEGENDRE) {
+      scatter[a] = new ScattDataLegendre;
     }
   }
 }
@@ -73,13 +63,13 @@ void
 XsData::from_hdf5(const hid_t xsdata_grp, const bool fissionable,
                   const int scatter_format, const int final_scatter_format,
                   const int order_data, const int max_order,
-                  const int legendre_to_tabular_points, const bool is_isotropic)
+                  const int legendre_to_tabular_points, const bool is_isotropic,
+                  const int n_pol, const int n_azi)
 {
   // 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();
+  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) {
@@ -98,11 +88,9 @@ XsData::from_hdf5(const hid_t xsdata_grp, const bool fissionable,
 
   // 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[p][a][gin] == 0.) absorption[p][a][gin] = 1.e-10;
-      }
+  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;
     }
   }
 
@@ -110,23 +98,17 @@ XsData::from_hdf5(const hid_t xsdata_grp, const bool fissionable,
   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];
-        }
+    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];
       }
     }
   }
 
-  // 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;
-      }
+  // 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;
     }
   }
 }
@@ -138,15 +120,14 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
      const int n_azi, const int energy_groups, const int delayed_groups,
      const 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_4dvec temp_beta =
-       double_4dvec(n_pol, double_3dvec(n_azi,
-       double_2dvec(energy_groups, double_1dvec(delayed_groups, 0.))));
+  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);
@@ -160,14 +141,12 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
 
       // 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];
-            }
+      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];
           }
         }
       }
@@ -181,41 +160,38 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
 
   // 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)));
+    double_2dvec temp_arr = double_2dvec(n_ang, double_1dvec(energy_groups));
     read_nd_vector(xsdata_grp, "chi", temp_arr);
 
-    for (int p = 0; p < n_pol; p++) {
-      for (int a = 0; a < n_azi; a++) {
-        // First set the first group
+    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++) {
-          chi_prompt[p][a][0][gout] = temp_arr[p][a][gout];
-        }
-
-        // 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];
-            }
+          for (int dg = 0; dg < delayed_groups; dg++) {
+            chi_delayed[a][gin][gout][dg] =
+                 chi_prompt[a][gin][gout];
           }
         }
       }
@@ -234,21 +210,18 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
       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];
-            }
+      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[p][a][gin].begin(),
-                                                temp_beta[p][a][gin].end(), 0.);
-              prompt_nu_fission[p][a][gin] *= (1. - beta_sum);
-            }
+          // 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);
           }
         }
       }
@@ -258,47 +231,45 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
       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.);
-            // Set the vector nu-fission from the matrix nu-fission
-            prompt_nu_fission[p][a][gin] = chi_sum;
+      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[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++) {
+          if (chi_sum >= 0.) {
             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];
-              }
+              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[p][a][gin][dg] =
-                   temp_beta[p][a][gin][dg] *
-                   prompt_nu_fission[p][a][gin];
-            }
+        // 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[p][a][gin].begin(),
-                                                temp_beta[p][a][gin].end(), 0.);
-              prompt_nu_fission[p][a][gin] *= (1. - beta_sum);
-            }
+          // 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);
           }
         }
       }
@@ -311,27 +282,24 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
 
   // 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)));
+    double_2dvec temp_arr = double_2dvec(n_ang, 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];
-          }
+    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[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.!");
+        // 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.!");
         }
       }
     }
@@ -346,26 +314,22 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
 
     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)));
+      double_2dvec temp_arr = double_2dvec(n_ang, 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!");
-          }
+      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[p][a][gin][gout][dg] =
-                     temp_arr[p][a][gout] / chi_sum;
-              }
+        // 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;
             }
           }
         }
@@ -375,22 +339,20 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
       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];
-              }
+      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[p][a][gin][gout][dg] /= chi_sum;
-                }
-              } else {
-                fatal_error("Encountered chi-delayed for a group that is <= 0!");
+            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!");
             }
           }
         }
@@ -414,30 +376,27 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
     } 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))));
+      double_3dvec temp_arr = double_3dvec(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_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;
-          }
+      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[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!");
+        // 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!");
           }
         }
       }
@@ -455,23 +414,20 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
     if (is_isotropic) ndims += 2;
 
     if (ndims == 3) {
-      // delayed-nu-fission is a [in group] vector
-      if (temp_beta[0][0][0][0] == 0.) {
+      // 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_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi,
-         double_1dvec(energy_groups)));
+      double_2dvec temp_arr = double_2dvec(n_ang, 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];
-            }
+      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];
           }
         }
       }
@@ -482,32 +438,28 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
 
     } 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)))));
+      double_4dvec temp_arr = double_4dvec(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 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.;
+      // 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++) {
-                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!");
+                chi_delayed[a][gin][gout][dg] /= gout_sum;
               }
+            } else {
+              fatal_error("Encountered chi-delayed for a group that is <= 0!");
             }
           }
         }
@@ -520,14 +472,12 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
   }
 
   // Combine prompt_nu_fission and delayed_nu_fission into nu_fission
-  for (int p = 0; p < n_pol; p++) {
-    for (int a = 0; a < n_azi; a++) {
-      for (int gin = 0; gin < energy_groups; gin++) {
-        nu_fission[p][a][gin] =
-             std::accumulate(delayed_nu_fission[p][a][gin].begin(),
-                             delayed_nu_fission[p][a][gin].end(),
-                             prompt_nu_fission[p][a][gin]);
-      }
+  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]);
     }
   }
 }
@@ -540,37 +490,32 @@ XsData::_scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
      const int final_scatter_format, const int order_data, const int max_order,
      const 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_3dvec gmin = int_3dvec(n_pol, int_2dvec(n_azi,
-       int_1dvec(energy_groups)));
+  int_2dvec gmin = int_2dvec(n_ang, 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)));
+  int_2dvec gmax = int_2dvec(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 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;
-      }
+  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 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);
-      }
+  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 = double_1dvec(length);
@@ -587,55 +532,48 @@ XsData::_scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
 
   // 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)));
+  double_4dvec input_scatt =
+       double_4dvec(n_ang, 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);
+  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_4dvec temp_mult = double_4dvec(n_pol, double_3dvec(n_azi,
-       double_2dvec(energy_groups)));
+  double_3dvec temp_mult = double_3dvec(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 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++];
-          }
+    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 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.;
-          }
+    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.;
         }
       }
     }
@@ -646,28 +584,22 @@ XsData::_scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
   // 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]);
+    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[p][a]),
-             legendre_to_tabular_points);
+      // Now create a tabular version of legendre_scatt
+      convert_legendre_to_tabular(legendre_scatt,
+           *static_cast<ScattDataTabular*>(scatter[a]),
+           legendre_to_tabular_points);
 
-        scatter_format = final_scatter_format;
-      }
+      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]);
-      }
+    for (int a = 0; a < n_ang; a++) {
+      scatter[a]->init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]);
     }
   }
 }
@@ -675,7 +607,7 @@ XsData::_scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
 //==============================================================================
 
 void
-XsData::combine(const std::vector<XsData*> those_xs,
+XsData::combine(const std::vector<XsData*>& those_xs,
                 const double_1dvec& scalars)
 {
   // Combine the non-scattering data
@@ -683,64 +615,60 @@ XsData::combine(const std::vector<XsData*> those_xs,
     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];
-          if (that->prompt_nu_fission.size() > 0) {
-            nu_fission[p][a][gin] +=
-                 scalar * that->nu_fission[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 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[p][a][gin].size(); dg++) {
-              delayed_nu_fission[p][a][gin][dg] +=
-                   scalar * that->delayed_nu_fission[p][a][gin][dg];
-            }
+          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[p][a][gin].size(); gout++) {
-              chi_prompt[p][a][gin][gout] +=
-                   scalar * that->chi_prompt[p][a][gin][gout];
+          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[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 < 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[p][a].size(); dg++) {
-          decay_rate[p][a][dg] += scalar * that->decay_rate[p][a][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[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;
-              }
+      // 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[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;
-                }
+          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;
               }
             }
           }
@@ -750,17 +678,15 @@ XsData::combine(const std::vector<XsData*> those_xs,
   }
 
   // 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);
+  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];
     }
+
+    // Now combine these guys
+    scatter[a]->combine(those_scatts, scalars);
   }
 }
 
@@ -770,12 +696,8 @@ 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())) {
+  if ((absorption.size() == that.absorption.size()) &&
+      (absorption[0].size() == that.absorption[0].size())) {
     match = true;
   }
   return match;
diff --git a/src/xsdata.h b/src/xsdata.h
index 4d3a5d5d4..dd18a8ce9 100644
--- a/src/xsdata.h
+++ b/src/xsdata.h
@@ -28,32 +28,34 @@ class XsData {
          const int n_azi, const int energy_groups, int scatter_format,
          const int final_scatter_format, const int order_data,
          const int max_order, const 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, bool is_isotropic);
+    void _fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
+         const int n_azi, const int energy_groups, const int delayed_groups,
+         const bool is_isotropic);
   public:
     // The following quantities have the following dimensions:
-    // [phi][theta][incoming group]
-    double_3dvec total;
-    double_3dvec absorption;
-    double_3dvec nu_fission;
-    double_3dvec prompt_nu_fission;
-    double_3dvec kappa_fission;
-    double_3dvec fission;
-    double_3dvec inverse_velocity;
+    // [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:
-    // [phi][theta][delayed group]
-    double_3dvec decay_rate;
+    // [angle][delayed group]
+    double_2dvec decay_rate;
     // delayed_nu_fission has the following dimensions:
-    // [phi][theta][incoming group][delayed group]
-    double_4dvec delayed_nu_fission;
+    // [angle][incoming group][delayed group]
+    double_3dvec delayed_nu_fission;
     // chi_prompt has the following dimensions:
-    // [phi][theta][incoming group][outgoing group]
-    double_4dvec chi_prompt;
+    // [angle][incoming group][outgoing group]
+    double_3dvec 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;
+    // [angle][incoming group][outgoing group][delayed group]
+    double_4dvec chi_delayed;
+    // scatter has the following dimensions: [angle]
+    std::vector<ScattData*> scatter;
 
     XsData() = default;
     XsData(const int num_groups, const int num_delayed_groups,
@@ -63,8 +65,8 @@ class XsData {
                    const int scatter_format, const int final_scatter_format,
                    const int order_data, const int max_order,
                    const int legendre_to_tabular_points,
-                   const bool is_isotropic);
-    void combine(const std::vector<XsData*> those_xs,
+                   const bool is_isotropic, const int n_pol, const int n_azi);
+    void combine(const std::vector<XsData*>& those_xs,
                  const double_1dvec& scalars);
     bool equiv(const XsData& that);
 };