From 58a5f871c443a66506db0f47a0ea3d987d196788 Mon Sep 17 00:00:00 2001
From: Paul Romano <romano7@gmail.com>
Date: Fri, 11 Nov 2011 13:10:41 -0500
Subject: [PATCH] Ability to sample multiple secondary energy distributions for
 single reaction. Closes gh-5.

---
 src/endf_header.f90   |  2 +-
 src/interpolation.f90 |  2 +-
 src/physics.f90       | 29 +++++++++++++++++++++--------
 3 files changed, 23 insertions(+), 10 deletions(-)

diff --git a/src/endf_header.f90 b/src/endf_header.f90
index bee61aa592..ce0a2929b2 100644
--- a/src/endf_header.f90
+++ b/src/endf_header.f90
@@ -7,7 +7,7 @@ module endf_header
 !===============================================================================
 
   type Tab1
-     integer :: n_regions           ! # of interpolation regions
+     integer :: n_regions = 0       ! # of interpolation regions
      integer, allocatable :: nbt(:) ! values separating interpolation regions
      integer, allocatable :: int(:) ! interpolation scheme
      integer :: n_pairs             ! # of pairs of (x,y) values
diff --git a/src/interpolation.f90 b/src/interpolation.f90
index e9d1f5337a..f32e1367d5 100644
--- a/src/interpolation.f90
+++ b/src/interpolation.f90
@@ -138,7 +138,7 @@ contains
 
     ! determine number of interpolation regions and pairs
     n_regions = obj % n_regions
-    n_pairs = obj % n_pairs
+    n_pairs   = obj % n_pairs
 
     ! find which bin the abscissa is in -- if the abscissa is outside the
     ! tabulated range, the first or last point is chosen, i.e. no interpolation
diff --git a/src/physics.f90 b/src/physics.f90
index 7086fd3bbc..3998c8e8c4 100644
--- a/src/physics.f90
+++ b/src/physics.f90
@@ -1425,7 +1425,7 @@ contains
 ! SAMPLE_ENERGY
 !===============================================================================
 
-  subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
+  recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
 
     type(DistEnergy),  pointer       :: edist
     real(8), intent(in)              :: E_in
@@ -1447,6 +1447,8 @@ contains
     integer :: ND          ! number of discrete lines
     integer :: NP          ! number of points in distribution
 
+    real(8) :: p_valid     ! probability of law validity
+
     real(8) :: E_i_1, E_i_K   ! endpoints on outgoing grid i
     real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
     real(8) :: E_1, E_K       ! endpoints interpolated between i and i+1
@@ -1479,15 +1481,26 @@ contains
     real(8) :: x, y, v     ! intermediate variables for n-body dist
     real(8) :: r1, r2, r3, r4, r5, r6
 
-    ! TODO: If there are multiple scattering laws, sample scattering law
+    ! ==========================================================================
+    ! SAMPLE ENERGY DISTRIBUTION IF THERE ARE MULTIPLE
 
-    ! Check for multiple interpolation regions
-    if (edist % p_valid % n_regions > 0) then
-       message = "Multiple interpolation regions not supported while &
-            &attempting to sample secondary energy distribution."
-       call fatal_error()
+    if (associated(edist % next)) then
+       if (edist % p_valid % n_regions > 0) then
+          p_valid = interpolate_tab1(edist % p_valid, E_in)
+
+          if (rang() > p_valid) then
+             if (edist % law == 44 .or. edist % law == 61) then
+                call sample_energy(edist%next, E_in, E_out, mu_out)
+             elseif (edist % law == 66) then
+                call sample_energy(edist%next, E_in, E_out, A=A, Q=Q)
+             else
+                call sample_energy(edist%next, E_in, E_out)
+             end if
+             return
+          end if
+       end if
     end if
-       
+
     ! Determine which secondary energy distribution law to use
     select case (edist % law)
     case (1)