diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index a3c9c07cf..8bfa11899 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -22,7 +22,7 @@ namespace openmc {
 //! @return The requested percentile
 //==============================================================================
 
-extern "C" double normal_percentile_c(double p);
+extern "C" double normal_percentile(double p);
 
 //==============================================================================
 //! Calculate the percentile of the Student's t distribution with a specified
@@ -58,7 +58,7 @@ extern "C" void calc_pn_c(int n, double x, double pnx[]);
 //!   evaluated at x
 //==============================================================================
 
-extern "C" double evaluate_legendre_c(int n, const double data[], double x);
+extern "C" double evaluate_legendre(int n, const double data[], double x);
 
 //==============================================================================
 //! Calculate the n-th order real spherical harmonics for a given angle (in
@@ -144,7 +144,7 @@ Direction rotate_angle(Direction u, double mu, double* phi);
 //! @result The sampled outgoing energy
 //==============================================================================
 
-extern "C" double maxwell_spectrum_c(double T);
+extern "C" double maxwell_spectrum(double T);
 
 //==============================================================================
 //! Samples an energy from a Watt energy-dependent fission distribution.
@@ -159,7 +159,7 @@ extern "C" double maxwell_spectrum_c(double T);
 //! @result The sampled outgoing energy
 //==============================================================================
 
-extern "C" double watt_spectrum_c(double a, double b);
+extern "C" double watt_spectrum(double a, double b);
 
 //==============================================================================
 //! Doppler broadens the windowed multipole curvefit.
diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h
index 02a43fece..e0addf5bc 100644
--- a/include/openmc/mgxs_interface.h
+++ b/include/openmc/mgxs_interface.h
@@ -18,6 +18,9 @@ extern std::vector<Mgxs> nuclides_MG;
 extern std::vector<Mgxs> macro_xs;
 extern "C" int num_energy_groups;
 
+//TODO: When more of the Fortran is converted (input_xml, tallies, etc, also
+// bring over energy_bin_avg, energy_bins, etc, as vectors)
+
 //==============================================================================
 // Mgxs data loading interface methods
 //==============================================================================
@@ -44,13 +47,6 @@ extern "C" void
 calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
      double& total_xs, double& abs_xs, double& nu_fiss_xs);
 
-extern "C" void
-sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
-     double uvw[3]);
-
-extern "C" void
-sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout);
-
 extern "C" double
 get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
 
diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h
index fdfebf71e..09a8cc617 100644
--- a/include/openmc/physics_mg.h
+++ b/include/openmc/physics_mg.h
@@ -10,22 +10,21 @@
 
 namespace openmc {
 
-//==============================================================================
-// SCATTER
-//==============================================================================
+//TODO: Remove energy_bin_avg and material_xs parameters when they reside on
+// the C-side this should happen after materials, physics, input, and tallies
+// are brought over
 
 //! \brief samples particle behavior after a collision event.
 extern "C" void
-collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size,
-     const double* energy_bin_avg, const MaterialMacroXS& material_xs);
+collision_mg(Particle* p, const double* energy_bin_avg,
+     const MaterialMacroXS& material_xs);
 
 //! \brief samples a reaction type.
 //!
 //! Note that there is special logic when suvival biasing is turned on since
 //! fission and disappearance are treated implicitly.
 void
-sample_reaction(Particle* p, Bank* fission_bank,
-     const int64_t fission_bank_size, const double* energy_bin_avg,
+sample_reaction(Particle* p, const double* energy_bin_avg,
      const MaterialMacroXS& material_xs);
 
 //! \brief Samples the scattering event
diff --git a/openmc/capi/math.py b/openmc/capi/math.py
index 390ee466f..a8ea849b3 100644
--- a/openmc/capi/math.py
+++ b/openmc/capi/math.py
@@ -12,8 +12,8 @@ _dll.t_percentile_c.argtypes = [c_double, c_int]
 _dll.calc_pn_c.restype = None
 _dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)]
 
-_dll.evaluate_legendre_c.restype = c_double
-_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double]
+_dll.evaluate_legendre.restype = c_double
+_dll.evaluate_legendre.argtypes = [c_int, POINTER(c_double), c_double]
 
 _dll.calc_rn_c.restype = None
 _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
@@ -27,11 +27,11 @@ _dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
 _dll.rotate_angle_c.restype = None
 _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
                                 POINTER(c_double)]
-_dll.maxwell_spectrum_c.restype = c_double
-_dll.maxwell_spectrum_c.argtypes = [c_double]
+_dll.maxwell_spectrum.restype = c_double
+_dll.maxwell_spectrum.argtypes = [c_double]
 
-_dll.watt_spectrum_c.restype = c_double
-_dll.watt_spectrum_c.argtypes = [c_double, c_double]
+_dll.watt_spectrum.restype = c_double
+_dll.watt_spectrum.argtypes = [c_double, c_double]
 
 _dll.broaden_wmp_polynomials_c.restype = None
 _dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int,
@@ -100,9 +100,8 @@ def evaluate_legendre(data, x):
     """
 
     data_arr = np.array(data, dtype=np.float64)
-    return _dll.evaluate_legendre_c(len(data),
-                                    data_arr.ctypes.data_as(POINTER(c_double)),
-                                    x)
+    return _dll.evaluate_legendre(len(data),
+                                  data_arr.ctypes.data_as(POINTER(c_double)), x)
 
 
 def calc_rn(n, uvw):
@@ -182,7 +181,7 @@ def calc_zn_rad(n, rho):
     zn_rad = np.zeros(num_bins, dtype=np.float64)
     _dll.calc_zn_rad_c(n, rho, zn_rad)
     return zn_rad
-    
+
 
 def rotate_angle(uvw0, mu, phi=None):
     """ Rotates direction cosines through a polar angle whose cosine is
@@ -231,7 +230,7 @@ def maxwell_spectrum(T):
 
     """
 
-    return _dll.maxwell_spectrum_c(T)
+    return _dll.maxwell_spectrum(T)
 
 
 def watt_spectrum(a, b):
@@ -251,7 +250,7 @@ def watt_spectrum(a, b):
 
     """
 
-    return _dll.watt_spectrum_c(a, b)
+    return _dll.watt_spectrum(a, b)
 
 
 def broaden_wmp_polynomials(E, dopp, n):
diff --git a/src/distribution.cpp b/src/distribution.cpp
index 01b09e65d..e546b11dd 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -90,7 +90,7 @@ Maxwell::Maxwell(pugi::xml_node node)
 
 double Maxwell::sample() const
 {
-  return maxwell_spectrum_c(theta_);
+  return maxwell_spectrum(theta_);
 }
 
 //==============================================================================
@@ -110,7 +110,7 @@ Watt::Watt(pugi::xml_node node)
 
 double Watt::sample() const
 {
-  return watt_spectrum_c(a_, b_);
+  return watt_spectrum(a_, b_);
 }
 
 //==============================================================================
diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index 863ffd8b0..e311e322d 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -279,7 +279,7 @@ double MaxwellEnergy::sample(double E) const
 
   while (true) {
     // Sample maxwell fission spectrum
-    double E_out = maxwell_spectrum_c(theta);
+    double E_out = maxwell_spectrum(theta);
 
     // Accept energy based on restriction energy
     if (E_out <= E - u_) return E_out;
@@ -343,7 +343,7 @@ double WattEnergy::sample(double E) const
 
   while (true) {
     // Sample energy-dependent Watt fission spectrum
-    double E_out = watt_spectrum_c(a, b);
+    double E_out = watt_spectrum(a, b);
 
     // Accept energy based on restriction energy
     if (E_out <= E - u_) return E_out;
diff --git a/src/math.F90 b/src/math.F90
index 9964922c1..143e21e03 100644
--- a/src/math.F90
+++ b/src/math.F90
@@ -12,10 +12,7 @@ module math
   public :: calc_rn
   public :: calc_zn
   public :: calc_zn_rad
-  public :: evaluate_legendre
   public :: rotate_angle
-  public :: maxwell_spectrum
-  public :: watt_spectrum
   public :: faddeeva
   public :: w_derivative
   public :: broaden_wmp_polynomials
@@ -42,16 +39,6 @@ module math
       real(C_DOUBLE), intent(out) :: pnx(n + 1)
     end subroutine calc_pn
 
-    pure function evaluate_legendre_c_intfc(n, data, x) &
-         bind(C, name='evaluate_legendre_c') result(val)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: n
-      real(C_DOUBLE), intent(in) :: data(n)
-      real(C_DOUBLE), value, intent(in) :: x
-      real(C_DOUBLE) :: val
-    end function evaluate_legendre_c_intfc
-
     pure subroutine calc_rn(n, uvw, rn) bind(C, name='calc_rn_c')
       use ISO_C_BINDING
       implicit none
@@ -85,23 +72,6 @@ module math
       real(C_DOUBLE), optional, intent(in) :: phi
     end subroutine rotate_angle_c_intfc
 
-    function maxwell_spectrum(T) bind(C, name='maxwell_spectrum_c') &
-         result(E_out)
-      use ISO_C_BINDING
-      implicit none
-      real(C_DOUBLE), value, intent(in) :: T
-      real(C_DOUBLE) :: E_out
-    end function maxwell_spectrum
-
-    function watt_spectrum(a, b) bind(C, name='watt_spectrum_c') &
-         result(E_out)
-      use ISO_C_BINDING
-      implicit none
-      real(C_DOUBLE), value, intent(in) :: a
-      real(C_DOUBLE), value, intent(in) :: b
-      real(C_DOUBLE) :: E_out
-    end function watt_spectrum
-
     subroutine broaden_wmp_polynomials(E, dopp, n, factors) &
            bind(C, name='broaden_wmp_polynomials_c')
       use ISO_C_BINDING
@@ -157,20 +127,6 @@ module math
 
 contains
 
-!===============================================================================
-! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
-! and the value of x
-!===============================================================================
-
-  pure function evaluate_legendre(data, x) result(val) bind(C)
-    real(C_DOUBLE), intent(in) :: data(:)
-    real(C_DOUBLE), intent(in) :: x
-    real(C_DOUBLE)             :: val
-
-    val = evaluate_legendre_c_intfc(size(data) - 1, data, x)
-
-  end function evaluate_legendre
-
 !===============================================================================
 ! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
 ! mu and through an azimuthal angle sampled uniformly. Note that this is done
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 5c1a2f4dc..8d30ea84d 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -6,7 +6,7 @@ namespace openmc {
 // Mathematical methods
 //==============================================================================
 
-double normal_percentile_c(double p) {
+double normal_percentile(double p) {
   constexpr double p_low = 0.02425;
   constexpr double a[6] = {-3.969683028665376e1, 2.209460984245205e2,
                            -2.759285104469687e2, 1.383577518672690e2,
@@ -79,7 +79,7 @@ double t_percentile_c(double p, int df){
     // 16 (4), pp. 1123-1132 (1987).
     double n = df;
     double k = 1. / (n - 2.);
-    double z = normal_percentile_c(p);
+    double z = normal_percentile(p);
     double z2 = z * z;
     t = std::sqrt(n * k) * (z + (z2 - 3.) * z * k / 4. + ((5. * z2 - 56.) * z2 +
          75.) * z * k * k / 96. + (((z2 - 27.) * 3. * z2 + 417.) * z2 - 315.) *
@@ -103,7 +103,7 @@ void calc_pn_c(int n, double x, double pnx[]) {
 }
 
 
-double evaluate_legendre_c(int n, const double data[], double x) {
+double evaluate_legendre(int n, const double data[], double x) {
   double pnx[n + 1];
   double val = 0.0;
   calc_pn_c(n, x, pnx);
@@ -658,7 +658,7 @@ Direction rotate_angle(Direction u, double mu, double* phi)
 }
 
 
-double maxwell_spectrum_c(double T) {
+double maxwell_spectrum(double T) {
   // Set the random numbers
   double r1 = prn();
   double r2 = prn();
@@ -674,8 +674,8 @@ double maxwell_spectrum_c(double T) {
 }
 
 
-double watt_spectrum_c(double a, double b) {
-  double w = maxwell_spectrum_c(a);
+double watt_spectrum(double a, double b) {
+  double w = maxwell_spectrum(a);
   double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w);
 
   return E_out;
diff --git a/src/mgxs_interface.F90 b/src/mgxs_interface.F90
index 55a18afab..77f629b80 100644
--- a/src/mgxs_interface.F90
+++ b/src/mgxs_interface.F90
@@ -62,26 +62,6 @@ module mgxs_interface
       real(C_DOUBLE),     intent(inout) :: nu_fiss_xs
     end subroutine calculate_xs_c
 
-    subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: i_mat
-      integer(C_INT), value, intent(in) :: gin
-      integer(C_INT),     intent(inout) :: gout
-      real(C_DOUBLE),     intent(inout) :: mu
-      real(C_DOUBLE),     intent(inout) :: wgt
-      real(C_DOUBLE),     intent(inout) :: uvw(1:3)
-    end subroutine sample_scatter_c
-
-    subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: i_mat
-      integer(C_INT), value, intent(in) :: gin
-      integer(C_INT),     intent(inout) :: dg
-      integer(C_INT),     intent(inout) :: gout
-    end subroutine sample_fission_energy_c
-
     subroutine get_name_c(index, name_len, name) bind(C)
       use ISO_C_BINDING
       implicit none
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index 56f3392ff..768ffe979 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -97,36 +97,6 @@ calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
 
 //==============================================================================
 
-void
-sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
-                 double uvw[3])
-{
-  int gout_c = gout - 1;
-  macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
-
-  // adjust return value for fortran indexing
-  gout = gout_c + 1;
-
-  // Rotate the angle
-  rotate_angle_c(uvw, mu, nullptr);
-}
-
-//==============================================================================
-
-void
-sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout)
-{
-  int dg_c = 0;
-  int gout_c = 0;
-  macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
-
-  // adjust return values for fortran indexing
-  dg = dg_c + 1;
-  gout = gout_c + 1;
-}
-
-//==============================================================================
-
 double
 get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
 {
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index 8651e29f1..10e6878b2 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -20,15 +20,14 @@
 namespace openmc {
 
 void
-collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size,
-     const double* energy_bin_avg, const MaterialMacroXS& material_xs)
+collision_mg(Particle* p, const double* energy_bin_avg,
+             const MaterialMacroXS& material_xs)
 {
   // Add to the collision counter for the particle
   p->n_collision++;
 
   // Sample the reaction type
-  sample_reaction(p, fission_bank, fission_bank_size, energy_bin_avg,
-                  material_xs);
+  sample_reaction(p, energy_bin_avg, material_xs);
 
   // Display information about collision
   if ((settings::verbosity >= 10) || (openmc_trace)) {
@@ -39,8 +38,7 @@ collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size,
 }
 
 void
-sample_reaction(Particle* p, Bank* fission_bank,
-                const int64_t fission_bank_size, const double* energy_bin_avg,
+sample_reaction(Particle* p, const double* energy_bin_avg,
                 const MaterialMacroXS& material_xs)
 {
   // Create fission bank sites. Note that while a fission reaction is sampled,
@@ -50,7 +48,11 @@ sample_reaction(Particle* p, Bank* fission_bank,
 
   if (materials[p->material - 1]->fissionable) {
     if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-      create_fission_sites(p, fission_bank, n_bank, fission_bank_size,
+      Bank* result_bank;
+      int64_t result_bank_size;
+      // Get pointer to fission bank from Fortran side
+      openmc_fission_bank(&result_bank, &result_bank_size);
+      create_fission_sites(p, result_bank, n_bank, result_bank_size,
                            material_xs);
     } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
                (settings::create_fission_neutrons)) {
diff --git a/src/scattdata.cpp b/src/scattdata.cpp
index 8d7138624..21e7f9df7 100644
--- a/src/scattdata.cpp
+++ b/src/scattdata.cpp
@@ -316,8 +316,8 @@ ScattDataLegendre::update_max_val()
         }
 
         // Calculate probability
-        double f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
-                                       dist[gin][i_gout].data(), mu);
+        double f = evaluate_legendre(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;
@@ -339,8 +339,8 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu)
     f = 0.;
   } else {
     int i_gout = gout - gmin[gin];
-    f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
-                            dist[gin][i_gout].data(), mu);
+    f = evaluate_legendre(dist[gin][i_gout].size() - 1,
+                          dist[gin][i_gout].data(), mu);
   }
   return f;
 }
@@ -881,8 +881,8 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
       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]);
+             evaluate_legendre(leg.dist[gin][i_gout].size() - 1,
+                               leg.dist[gin][i_gout].data(), tab.mu[imu]);
       }
 
       // Ensure positivity
diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp
index 45f6d3cb8..dabd9de3a 100644
--- a/src/secondary_nbody.cpp
+++ b/src/secondary_nbody.cpp
@@ -32,13 +32,13 @@ void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu) const
   double E_max = (Ap - 1.0)/Ap * (A_/(A_ + 1.0)*E_in + Q_);
 
   // x is essentially a Maxwellian distribution
-  double x = maxwell_spectrum_c(1.0);
+  double x = maxwell_spectrum(1.0);
 
   double y;
   double r1, r2, r3, r4, r5, r6;
   switch (n_bodies_) {
   case 3:
-    y = maxwell_spectrum_c(1.0);
+    y = maxwell_spectrum(1.0);
     break;
   case 4:
     r1 = prn();
diff --git a/src/tracking.F90 b/src/tracking.F90
index 02f10f7cc..982daeccc 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -34,12 +34,9 @@ module tracking
   implicit none
 
   interface
-    subroutine collision_mg(p, fission_bank, fission_bank_size, &
-         energy_bin_avg, material_xs) bind(C)
-      import Particle, Bank, C_INT64_T, C_DOUBLE, MaterialMacroXS
+    subroutine collision_mg(p, energy_bin_avg, material_xs) bind(C)
+      import Particle, C_DOUBLE, MaterialMacroXS
       type(Particle),            intent(inout) :: p
-      type(Bank),                intent(inout) :: fission_bank(*)
-      integer(C_INT64_T), value, intent(in)    :: fission_bank_size
       real(C_DOUBLE),            intent(in)    :: energy_bin_avg(*)
       type(MaterialMacroXS),     intent(in)    :: material_xs
     end subroutine collision_mg
@@ -239,8 +236,7 @@ contains
         if (run_CE) then
           call collision(p)
         else
-          call collision_mg(p, fission_bank, &
-               size(fission_bank, kind=C_INT64_T), energy_bin_avg, material_xs)
+          call collision_mg(p, energy_bin_avg, material_xs)
         end if
 
         ! Score collision estimator tallies -- this is done after a collision