diff --git a/src/constants.h b/src/constants.h
index 4129305bf..8290aa8c9 100644
--- a/src/constants.h
+++ b/src/constants.h
@@ -10,7 +10,6 @@
 
 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;
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index 42df45ef0..ec711e130 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -24,6 +24,13 @@ void Mgxs::init(const std::string& in_name, const double in_awr,
   azimuthal = in_azimuthal;
   n_pol = polar.size();
   n_azi = azimuthal.size();
+  index_temp = 0;
+  last_sqrtkT = 0.;
+  index_pol = 0;
+  index_azi = 0;
+  last_uvw[0] = 1.;
+  last_uvw[1] = 0.;
+  last_uvw[2] = 0.;
 }
 
 
@@ -383,7 +390,8 @@ void Mgxs::combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
 }
 
 
-double Mgxs::get_xs(const char* xstype, int gin, int* gout, double* mu, int* dg)
+double Mgxs::get_xs(const char* xstype, const int gin, int* gout, double* mu,
+                    int* dg)
 {
   // This method assumes that the temperature and angle indices are set
   double val;
@@ -469,7 +477,8 @@ double Mgxs::get_xs(const char* xstype, int gin, int* gout, double* mu, int* dg)
 }
 
 
-void Mgxs::sample_fission_energy(int gin, double nu_fission, int& dg, int& gout)
+void Mgxs::sample_fission_energy(const int gin, const double nu_fission,
+                                 int& dg, int& gout)
 {
   // This method assumes that the temperature and angle indices are set
   // Find the probability of having a prompt neutron
@@ -525,8 +534,7 @@ void Mgxs::sample_fission_energy(int gin, double nu_fission, int& dg, int& gout)
 }
 
 
-void Mgxs::sample_scatter(dir_arr& uvw, int gin, int& gout, double& mu,
-                          double& wgt)
+void Mgxs::sample_scatter(const int gin, int& gout, double& mu, double& wgt)
 {
   // This method assumes that the temperature and angle indices are set
   // Sample the data
@@ -534,8 +542,8 @@ void Mgxs::sample_scatter(dir_arr& uvw, int gin, int& gout, double& mu,
 }
 
 
-void Mgxs::calculate_xs(int gin, double sqrtkT, dir_arr& uvw, double& total_xs,
-                        double& abs_xs, double& nu_fiss_xs)
+void Mgxs::calculate_xs(const int gin, const double sqrtkT, const double uvw[3],
+                        double& total_xs, double& abs_xs, double& nu_fiss_xs)
 {
   // Set our indices
   set_temperature_index(sqrtkT);
@@ -543,10 +551,14 @@ void Mgxs::calculate_xs(int gin, double sqrtkT, dir_arr& uvw, double& total_xs,
   total_xs = xs[index_temp].total[index_pol][index_azi][gin];
   abs_xs = xs[index_temp].absorption[index_pol][index_azi][gin];
 
-  // nu-fission is made up of the prompt and all the delayed nu_fission data
-  nu_fiss_xs = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
-  for (auto& val : xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin]) {
-    nu_fiss_xs += val;
+  if (fissionable) {
+    // nu-fission is made up of the prompt and all the delayed nu_fission data
+    nu_fiss_xs = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
+    for (auto& val : xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin]) {
+      nu_fiss_xs += val;
+    }
+  } else {
+    nu_fiss_xs = 0.;
   }
 }
 
@@ -568,7 +580,7 @@ bool Mgxs::equiv(const Mgxs& that)
 }
 
 
-inline void Mgxs::set_temperature_index(double sqrtkT)
+inline void Mgxs::set_temperature_index(const double sqrtkT)
 {
   // See if we need to find the new index
   if (sqrtkT != last_sqrtkT) {
@@ -588,27 +600,30 @@ inline void Mgxs::set_temperature_index(double sqrtkT)
 }
 
 
-inline void Mgxs::set_angle_index(dir_arr& uvw)
+inline void Mgxs::set_angle_index(const double uvw[3])
 {
   // See if we need to find the new index
-  if (uvw != last_uvw) {
+  if ((uvw[0] != last_uvw[0]) || (uvw[1] != last_uvw[1]) ||
+      (uvw[2] != last_uvw[2])) {
     // 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;
-    index_pol = std::floor(my_pol / delta_angle + 1.);
-    delta_angle = PI / n_azi;
-    index_azi = std::floor((my_azi + PI) / delta_angle + 1.);
+    index_pol = std::floor(my_pol / delta_angle);
+    delta_angle = 2. * PI / n_azi;
+    index_azi = std::floor((my_azi + PI) / delta_angle);
 
     // store this direction as the last one used
-    last_uvw = uvw;
+    last_uvw[0] = uvw[0];
+    last_uvw[1] = uvw[1];
+    last_uvw[2] = uvw[2];
   }
 }
 
 //==============================================================================
-// Mgxs data loading methods
+// Mgxs data loading interface methods
 //==============================================================================
 
 void add_mgxs(hid_t file_id, char* name, int energy_groups,
@@ -680,4 +695,25 @@ void create_macro_xs(char* mat_name, const int n_nuclides,
   macro_xs.push_back(macro);
 }
 
+//==============================================================================
+// Mgxs tracking/transport/tallying interface methods
+//==============================================================================
+
+void calculate_xs(const int i_mat, const int gin, const double sqrtkT,
+     const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs)
+{
+  macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
+       nu_fiss_xs);
+}
+
+void scatter(const int i_mat, const int gin, int& gout, double& mu,
+             double& wgt, double uvw[3])
+{
+  int gout_c = gout - 1;
+  macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
+  gout = gout_c + 1;
+
+  rotate_angle_c(uvw, mu, nullptr);
+}
+
 } // namespace openmc
\ No newline at end of file
diff --git a/src/mgxs.h b/src/mgxs.h
index 90ed520d9..14c4e4e89 100644
--- a/src/mgxs.h
+++ b/src/mgxs.h
@@ -44,7 +44,7 @@ class Mgxs {
     int index_azi;
     double_1dvec polar;
     double_1dvec azimuthal;
-    dir_arr last_uvw;
+    double last_uvw[3];
     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,
@@ -66,16 +66,16 @@ class Mgxs {
          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,
+    double get_xs(const char* xstype, const 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);
+    void sample_fission_energy(const int gin, const double nu_fission,
+         int& dg, int& gout);
+    void sample_scatter(const int gin, int& gout, double& mu, double& wgt);
+    void calculate_xs(const int gin, const double sqrtkT, const double uvw[3],
+         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);
+    inline void set_temperature_index(const double sqrtkT);
+    inline void set_angle_index(const double uvw[3]);
 };
 
 extern "C" void add_mgxs(hid_t file_id, char* name, int energy_groups,
@@ -89,6 +89,13 @@ extern "C" void create_macro_xs(char* mat_name, const int n_nuclides,
      const int i_nuclides[], const int n_temps, const double temps[],
      const double atom_densities[], int& method, const double tolerance);
 
+extern "C" void calculate_xs(const int i_mat, const int gin,
+     const double sqrtkT, const double uvw[3], double& total_xs,
+     double& abs_xs, double& nu_fiss_xs);
+
+extern "C" void scatter(const int i_mat, const int gin, int& gout, double& mu,
+     double& wgt, double uvw[3]);
+
 
 // Storage for the MGXS data
 std::vector<Mgxs> nuclides_MG;
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index ebabe5cdb..94aa1928a 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -166,10 +166,10 @@ module nuclide_header
 !===============================================================================
 
   type MaterialMacroXS
-    real(8) :: total         ! macroscopic total xs
-    real(8) :: absorption    ! macroscopic absorption xs
-    real(8) :: fission       ! macroscopic fission xs
-    real(8) :: nu_fission    ! macroscopic production xs
+    real(C_DOUBLE) :: total         ! macroscopic total xs
+    real(C_DOUBLE) :: absorption    ! macroscopic absorption xs
+    real(C_DOUBLE) :: fission       ! macroscopic fission xs
+    real(C_DOUBLE) :: nu_fission    ! macroscopic production xs
   end type MaterialMacroXS
 
 !===============================================================================
diff --git a/src/physics_mg.F90 b/src/physics_mg.F90
index eb8208594..3030fc099 100644
--- a/src/physics_mg.F90
+++ b/src/physics_mg.F90
@@ -22,6 +22,20 @@ module physics_mg
 
   implicit none
 
+  interface
+    subroutine scatter_c(i_mat, gin, gout, mu, wgt, uvw) &
+         bind(C, name='scatter')
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: i_mat
+      integer(C_INT), value, intent(in) :: gin
+      integer(C_INT),     intent(inout) :: gout
+      real(C_DOUBLE),     intent(inout) :: mu
+      real(C_DOUBLE),     intent(inout) :: wgt
+      real(C_DOUBLE),     intent(inout) :: uvw(1:3)
+    end subroutine scatter_c
+  end interface
+
 contains
 
 !===============================================================================
@@ -143,16 +157,18 @@ contains
 
     type(Particle), intent(inout)  :: p
 
-    call macro_xs(p % material) % obj % sample_scatter(p % coord(1) % uvw, &
-                                                       p % last_g, p % g, &
-                                                       p % mu, p % wgt)
+    ! call macro_xs(p % material) % obj % sample_scatter(p % coord(1) % uvw, &
+    !      p % last_g, p % g, p % mu, p % wgt)
+
+    ! Convert change in angle (mu) to new direction
+    ! p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu)
+
+    call scatter_c(p % material, p % last_g, p % g, p % mu, p % wgt, &
+                   p % coord(1) % uvw)
 
     ! Update energy value for downstream compatability (in tallying)
     p % E = energy_bin_avg(p % g)
 
-    ! Convert change in angle (mu) to new direction
-    p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu)
-
     ! Set event component
     p % event = EVENT_SCATTER
 
diff --git a/src/scattdata.cpp b/src/scattdata.cpp
index d23169fa8..aef620ac5 100644
--- a/src/scattdata.cpp
+++ b/src/scattdata.cpp
@@ -44,7 +44,8 @@ void ScattData::sample_energy(int gin, int& gout, int& i_gout)
 {
   // Sample the outgoing group
   double xi = prn();
-  i_gout = 0; //TODO: + 1?
+
+  i_gout = 0;
   gout = gmin[gin];
   double prob = energy[gin][i_gout];
   while((prob < xi) && (gout < gmax[gin])) {
@@ -247,7 +248,7 @@ void ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt)
 }
 
 
-void ScattDataLegendre::combine(std::vector<ScattData*> those_scatts,
+void ScattDataLegendre::combine(std::vector<ScattData*>& those_scatts,
                                 double_1dvec& scalars)
 {
   // Find the max order in the data set and make sure we can combine the sets
@@ -294,6 +295,7 @@ void ScattDataLegendre::combine(std::vector<ScattData*> those_scatts,
     for (int gin = 0; gin < groups; gin++) {
       // Only spend time adding that's gmin to gmax data since the rest will
       // be zeros
+      int i_gout = 0;
       for (int gout = that->gmin[gin]; gout <= that->gmax[gin]; gout++) {
         // Do the scattering matrix
         for (int l = 0; l < max_order; l++) {
@@ -301,13 +303,14 @@ void ScattDataLegendre::combine(std::vector<ScattData*> those_scatts,
         }
 
         // Incorporate that's contribution to the multiplicity matrix data
-        double nuscatt = that->scattxs[gin] * that->energy[gin][gout];
+        double nuscatt = that->scattxs[gin] * that->energy[gin][i_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];
+        if (that->mult[gin][i_gout] > 0.) {
+          mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
         } else {
           mult_denom[gin][gout] += scalars[i];
         }
+        i_gout++;
       }
     }
   }
@@ -336,14 +339,14 @@ void ScattDataLegendre::combine(std::vector<ScattData*> those_scatts,
     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.;});
+                                  [](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.;});
+                                  [](double val){return val != 0.;});
       if (non_zero) break;
     }
 
@@ -425,6 +428,7 @@ void ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
     for (int i_gout = 0; i_gout < num_groups; i_gout++) {
       double norm = std::accumulate(matrix[gin][i_gout].begin(),
                                     matrix[gin][i_gout].end(), 0.);
+      in_energy[gin][i_gout] = norm;
       if (norm != 0.) {
         for (auto& n : matrix[gin][i_gout]) n /= norm;
       }
@@ -440,7 +444,7 @@ void ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
   dmu = 2. / order;
   mu[0] = -1.;
   for (int imu = 1; imu < order; imu++) {
-    mu[imu] = -1. + (imu - 1) * dmu;
+    mu[imu] = -1. + imu * dmu;
   }
 
   // Calculate f(mu) and integrate it so we can avoid rejection sampling
@@ -507,7 +511,7 @@ void ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
 
   int imu;
   if (xi < dist[gin][i_gout][0]) {
-    imu = 1;
+    imu = 0;
   } else {
     // TODO lower_bound?  + 1?
     imu = std::upper_bound(dist[gin][i_gout].begin(),
@@ -551,7 +555,7 @@ double_3dvec ScattDataHistogram::get_matrix(int max_order)
 }
 
 
-void ScattDataHistogram::combine(std::vector<ScattData*> those_scatts,
+void ScattDataHistogram::combine(std::vector<ScattData*>& those_scatts,
                                 double_1dvec& scalars)
 {
   // Find the max order in the data set and make sure we can combine the sets
@@ -600,6 +604,7 @@ void ScattDataHistogram::combine(std::vector<ScattData*> those_scatts,
     for (int gin = 0; gin < groups; gin++) {
       // Only spend time adding that's gmin to gmax data since the rest will
       // be zeros
+      int i_gout = 0;
       for (int gout = that->gmin[gin]; gout <= that->gmax[gin]; gout++) {
         // Do the scattering matrix
         for (int l = 0; l < max_order; l++) {
@@ -607,13 +612,14 @@ void ScattDataHistogram::combine(std::vector<ScattData*> those_scatts,
         }
 
         // Incorporate that's contribution to the multiplicity matrix data
-        double nuscatt = that->scattxs[gin] * that->energy[gin][gout];
+        double nuscatt = that->scattxs[gin] * that->energy[gin][i_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];
+        if (that->mult[gin][i_gout] > 0.) {
+          mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
         } else {
           mult_denom[gin][gout] += scalars[i];
         }
+        i_gout++;
       }
     }
   }
@@ -642,14 +648,14 @@ void ScattDataHistogram::combine(std::vector<ScattData*> those_scatts,
     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.;});
+                                  [](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.;});
+                                  [](double val){return val != 0.;});
       if (non_zero) break;
     }
 
@@ -695,7 +701,7 @@ void ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
   dmu = 2. / (order - 1);
   mu[0] = -1.;
   for (int imu = 1; imu < order - 1; imu++) {
-    mu[imu] = -1. + (imu - 1) * dmu;
+    mu[imu] = -1. + imu * dmu;
   }
   mu[order - 1] = 1.;
 
@@ -724,9 +730,7 @@ void ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
         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;
-      }
+      in_energy[gin][i_gout] = norm;
     }
   }
 
@@ -806,14 +810,14 @@ void ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
 
   double c_k = dist[gin][i_gout][0];
   int k;
-  for (k = 0; k < NP - 2; k++) {
+  for (k = 0; k < NP - 1; k++) {
     double c_k1 = dist[gin][i_gout][k + 1];
     if (xi < c_k1) break;
     c_k = c_k1;
   }
 
   // Check to make sure k is <= NP - 1
-  k = std::min(k, NP - 1);
+  k = std::min(k, NP - 2);
 
   // Find the pdf values we want
   double p0 = fmu[gin][i_gout][k];
@@ -861,7 +865,7 @@ double_3dvec ScattDataTabular::get_matrix(int max_order)
   return matrix;
 }
 
-void ScattDataTabular::combine(std::vector<ScattData*> those_scatts,
+void ScattDataTabular::combine(std::vector<ScattData*>& those_scatts,
                                double_1dvec& scalars)
 {
   // Find the max order in the data set and make sure we can combine the sets
@@ -910,6 +914,7 @@ void ScattDataTabular::combine(std::vector<ScattData*> those_scatts,
     for (int gin = 0; gin < groups; gin++) {
       // Only spend time adding that's gmin to gmax data since the rest will
       // be zeros
+      int i_gout = 0;
       for (int gout = that->gmin[gin]; gout <= that->gmax[gin]; gout++) {
         // Do the scattering matrix
         for (int l = 0; l < max_order; l++) {
@@ -917,13 +922,14 @@ void ScattDataTabular::combine(std::vector<ScattData*> those_scatts,
         }
 
         // Incorporate that's contribution to the multiplicity matrix data
-        double nuscatt = that->scattxs[gin] * that->energy[gin][gout];
+        double nuscatt = that->scattxs[gin] * that->energy[gin][i_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];
+        if (that->mult[gin][i_gout] > 0.) {
+          mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
         } else {
           mult_denom[gin][gout] += scalars[i];
         }
+        i_gout++;
       }
     }
   }
@@ -952,14 +958,14 @@ void ScattDataTabular::combine(std::vector<ScattData*> those_scatts,
     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.;});
+                                  [](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.;});
+                                  [](double val){return val != 0.;});
       if (non_zero) break;
     }
 
diff --git a/src/scattdata.h b/src/scattdata.h
index 6e24a50bc..4553156e0 100644
--- a/src/scattdata.h
+++ b/src/scattdata.h
@@ -41,7 +41,7 @@ class ScattData {
     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,
+    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;
@@ -59,7 +59,7 @@ class ScattDataLegendre: public ScattData {
     void update_max_val();
     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);
+    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);
 };
@@ -74,7 +74,7 @@ class ScattDataHistogram: public ScattData {
               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);
+    void combine(std::vector<ScattData*>& those_scatts, double_1dvec& scalars);
     int get_order() {return dist[0][0].size();};
     double_3dvec get_matrix(int max_order);
 };
@@ -91,7 +91,7 @@ class ScattDataTabular: public ScattData {
               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);
+    void combine(std::vector<ScattData*>& those_scatts, double_1dvec& scalars);
     int get_order() {return dist[0][0].size();};
     double_3dvec get_matrix(int max_order);
 };
diff --git a/src/tracking.F90 b/src/tracking.F90
index 3fe6a065b..2ba62a8dd 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -1,5 +1,7 @@
 module tracking
 
+  use, intrinsic :: ISO_C_BINDING
+
   use constants
   use error,              only: warning, write_message
   use geometry_header,    only: cells
@@ -27,6 +29,21 @@ module tracking
 
   implicit none
 
+  interface
+    subroutine calculate_xs_c(i_mat, gin, sqrtkT, uvw, total_xs, abs_xs, &
+         nu_fiss_xs) bind(C, name='calculate_xs')
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: i_mat
+      integer(C_INT), value, intent(in) :: gin
+      real(C_DOUBLE), value, intent(in) :: sqrtkT
+      real(C_DOUBLE),        intent(in) :: uvw(1:3)
+      real(C_DOUBLE),     intent(inout) :: total_xs
+      real(C_DOUBLE),     intent(inout) :: abs_xs
+      real(C_DOUBLE),     intent(inout) :: nu_fiss_xs
+    end subroutine calculate_xs_c
+  end interface
+
 contains
 
 !===============================================================================
@@ -112,8 +129,14 @@ contains
           end if
         else
           ! Get the MG data
+          !!TODO: Remove Fortran call - needed until I'm done replacing Fortran
+          !!with C++ code because it sets index_temp
           call macro_xs(p % material) % obj % calculate_xs(p % g, p % sqrtkT, &
                p % coord(p % n_coord) % uvw, material_xs)
+          call calculate_xs_c(p % material, p % g, p % sqrtkT, &
+               p % coord(p % n_coord) % uvw, material_xs % total, &
+               material_xs % absorption, material_xs % nu_fission)
+
 
           ! Finally, update the particle group while we have already checked
           ! for if multi-group