From 6bfc2f73e0777093f3f14461cabc023a2bb4eca9 Mon Sep 17 00:00:00 2001
From: Jon Walsh <walshjon@mit.edu>
Date: Mon, 28 Jul 2014 23:22:46 -0600
Subject: [PATCH] account for differences in band probabilities at different
 energies

---
 src/cross_section.F90 | 50 ++++++++++++++++++++++++-------------------
 1 file changed, 28 insertions(+), 22 deletions(-)

diff --git a/src/cross_section.F90 b/src/cross_section.F90
index 7426519cda..cdb4b1e2cc 100644
--- a/src/cross_section.F90
+++ b/src/cross_section.F90
@@ -340,7 +340,8 @@ contains
     real(8), intent(in) :: E         ! energy
 
     integer :: i_energy   ! index for energy
-    integer :: i_table    ! index for table
+    integer :: i_low      ! band index at lower bounding energy
+    integer :: i_up       ! band index at upper bounding energy
     real(8) :: f          ! interpolation factor
     real(8) :: r          ! pseudo-random number
     real(8) :: elastic    ! elastic cross section
@@ -371,21 +372,26 @@ contains
 
     ! sample probability table using the cumulative distribution
     r = prn()
-    i_table = 1
+    i_low = 1
     do
-      if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit
-      i_table = i_table + 1
+      if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit
+      i_low = i_low + 1
+    end do
+    i_up = 1
+    do
+      if (urr % prob(i_energy + 1, URR_CUM_PROB, i_up) > r) exit
+      i_up = i_up + 1
     end do
 
     ! determine elastic, fission, and capture cross sections from probability
     ! table
     if (urr % interp == LINEAR_LINEAR) then
-      elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_table) + &
-           f * urr % prob(i_energy + 1, URR_ELASTIC, i_table)
-      fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_table) + &
-           f * urr % prob(i_energy + 1, URR_FISSION, i_table)
-      capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_table) + &
-           f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table)
+      elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_low) + &
+           f * urr % prob(i_energy + 1, URR_ELASTIC, i_up)
+      fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_low) + &
+           f * urr % prob(i_energy + 1, URR_FISSION, i_up)
+      capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_low) + &
+           f * urr % prob(i_energy + 1, URR_N_GAMMA, i_up)
     elseif (urr % interp == LOG_LOG) then
       ! Get logarithmic interpolation factor
       f = log(E / urr % energy(i_energy)) / &
@@ -393,29 +399,29 @@ contains
 
       ! Calculate elastic cross section/factor
       elastic = ZERO
-      if (urr % prob(i_energy, URR_ELASTIC, i_table) > ZERO .and. &
-          urr % prob(i_energy + 1, URR_ELASTIC, i_table) > ZERO) then
+      if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. &
+          urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
         elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, &
-             i_table)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, &
-             i_table)))
+             i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, &
+             i_up)))
       end if
 
       ! Calculate fission cross section/factor
       fission = ZERO
-      if (urr % prob(i_energy, URR_FISSION, i_table) > ZERO .and. &
-          urr % prob(i_energy + 1, URR_FISSION, i_table) > ZERO) then
+      if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. &
+          urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
         fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, &
-             i_table)) + f * log(urr % prob(i_energy + 1, URR_FISSION, &
-             i_table)))
+             i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, &
+             i_up)))
       end if
 
       ! Calculate capture cross section/factor
       capture = ZERO
-      if (urr % prob(i_energy, URR_N_GAMMA, i_table) > ZERO .and. &
-          urr % prob(i_energy + 1, URR_N_GAMMA, i_table) > ZERO) then
+      if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. &
+          urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
         capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, &
-             i_table)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, &
-             i_table)))
+             i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, &
+             i_up)))
       end if
     end if