diff --git a/include/openmc/physics.h b/include/openmc/physics.h
index 22cb3f30c..87e7080c2 100644
--- a/include/openmc/physics.h
+++ b/include/openmc/physics.h
@@ -47,7 +47,7 @@ void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
 
 Reaction* sample_fission(int i_nuclide, double E);
 
-// void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
+extern "C" void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
 
 void absorption(Particle* p, int i_nuclide);
 
@@ -68,7 +68,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank
 
 // void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle* p);
 
-extern "C" void sample_secondary_photons(Particle* p, int i_nuclide);
+void sample_secondary_photons(Particle* p, int i_nuclide);
 
 extern "C" void thick_target_bremsstrahlung(Particle* p, double* E_lost);
 
diff --git a/src/physics.F90 b/src/physics.F90
index ffc74e78f..bcee5d60a 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -258,11 +258,11 @@ contains
 ! SAMPLE_PHOTON_PRODUCT
 !===============================================================================
 
-  subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product)
-    integer, intent(in)  :: i_nuclide  ! index in nuclides array
-    real(8), intent(in)  :: E          ! energy of neutron
-    integer, intent(out) :: i_reaction ! index in nuc % reactions array
-    integer, intent(out) :: i_product  ! index in reaction % products array
+  subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product) bind(C)
+    integer(C_INT), value  :: i_nuclide  ! index in nuclides array
+    real(C_DOUBLE), value  :: E          ! energy of neutron
+    integer(C_INT), intent(out) :: i_reaction ! index in nuc % reactions array
+    integer(C_INT), intent(out) :: i_product  ! index in reaction % products array
 
     integer :: i_grid
     integer :: i_temp
@@ -885,50 +885,4 @@ contains
 
   end subroutine inelastic_scatter
 
-!===============================================================================
-! SAMPLE_SECONDARY_PHOTONS
-!===============================================================================
-
-  subroutine sample_secondary_photons(p, i_nuclide) bind(C)
-    type(Particle), intent(inout) :: p
-    integer(C_INT), value         :: i_nuclide
-
-    integer :: i_reaction   ! index in nuc % reactions array
-    integer :: i_product    ! index in nuc % reactions % products array
-
-    real(8) :: nu_t
-    real(8) :: mu
-    real(8) :: E
-    real(8) :: uvw(3)
-    integer :: nu
-    integer :: i
-
-    ! Sample the number of photons produced
-    nu_t =  p % wgt * micro_xs(i_nuclide) % photon_prod / &
-         micro_xs(i_nuclide) % total
-    if (prn() > nu_t - int(nu_t)) then
-      nu = int(nu_t)
-    else
-      nu = int(nu_t) + 1
-    end if
-
-    ! Sample each secondary photon
-    do i = 1, nu
-
-      ! Sample the reaction and product
-      call sample_photon_product(i_nuclide, p % E, i_reaction, i_product)
-
-      ! Sample the outgoing energy and angle
-      call nuclides(i_nuclide) % reactions(i_reaction) % &
-           product_sample(i_product, p % E, E, mu)
-
-      ! Sample the new direction
-      uvw = rotate_angle(p % coord(1) % uvw, mu)
-
-      ! Create the secondary photon
-      call particle_create_secondary(p, uvw, E, PHOTON, run_CE=.true._C_BOOL)
-    end do
-
-  end subroutine sample_secondary_photons
-
 end module physics
diff --git a/src/physics.cpp b/src/physics.cpp
index 25e811f30..0abdd29c3 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -4,6 +4,7 @@
 #include "openmc/constants.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
+#include "openmc/math_functions.h"
 #include "openmc/message_passing.h"
 #include "openmc/nuclide.h"
 #include "openmc/photon.h"
@@ -1127,33 +1128,36 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank
 //   }
 // }
 
-// void sample_secondary_photons(Particle* p, int i_nuclide)
-// {
-//   // Sample the number of photons produced
-//   nu_t =  p->wgt * simulation::micro_xs[i_nuclide-1].photon_prod / &
-//         simulation::micro_xs[i_nuclide-1].total
-//   if (prn() > nu_t - int(nu_t)) {
-//     nu = int(nu_t)
-//   } else {
-//     nu = int(nu_t) + 1
-//   }
+void sample_secondary_photons(Particle* p, int i_nuclide)
+{
+  // Sample the number of photons produced
+  double y_t = p->wgt * simulation::micro_xs[i_nuclide-1].photon_prod /
+    simulation::micro_xs[i_nuclide-1].total;
+  int y = static_cast<int>(y_t);
+  if (prn() <= y_t - y) ++y;
 
-//   // Sample each secondary photon
-//   do i = 1, nu
+  // Sample each secondary photon
+  for (int i = 0; i < y; ++i) {
+    // Sample the reaction and product
+    int i_rx;
+    int i_product;
+    sample_photon_product(i_nuclide, p->E, &i_rx, &i_product);
 
-//     // Sample the reaction and product
-//     sample_photon_product(i_nuclide, p->E, i_reaction, i_product)
+    // Sample the outgoing energy and angle
+    auto& rx = data::nuclides[i_nuclide-1]->reactions_[i_rx-1];
+    double E;
+    double mu;
+    rx->products_[i_product-1].sample(p->E, E, mu);
 
-//     // Sample the outgoing energy and angle
-//     nuclides(i_nuclide) % reactions(i_reaction) % &
-//           product_sample(i_product, p->E, E, mu)
+    // Sample the new direction
+    double uvw[3];
+    std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw);
+    rotate_angle_c(uvw, mu, nullptr);
 
-//     // Sample the new direction
-//     uvw = rotate_angle(p->coord(1) % uvw, mu)
-
-//     // Create the secondary photon
-//     particle_create_secondary(p, uvw, E, PHOTON, run_CE=true)
-//   end do
-// }
+    // Create the secondary photon
+    int photon = static_cast<int>(ParticleType::photon);
+    p->create_secondary(uvw, E, photon, true);
+  }
+}
 
 } // namespace openmc