diff --git a/ChangeLog b/ChangeLog index ccadd40f5..7f6e23350 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,10 @@ +2011-03-25 Paul Romano + + * physics.f90: Resample outgoing energy if sampled energy is + greater than 20 MeV. Some cross sections go up to 150 MeV but + others do not. + * source.f90: Resample outgoing energy if >= 20 MeV. + 2011-03-25 Paul Romano * ace.f90: Added ability to read delayed neutron energy diff --git a/src/physics.f90 b/src/physics.f90 index a39252839..385915a6e 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -116,7 +116,7 @@ contains if (E < e_grid(1)) then IE = 1 elseif (E > e_grid(n_grid)) then - IE = n_grid + IE = n_grid - 1 else IE = binary_search(e_grid, n_grid, E) end if @@ -572,11 +572,15 @@ contains ! sample from energy distribution for group j law = table % nu_d_edist(j) % law - if (law == 44 .or. law == 61) then - call sample_energy(table%nu_d_edist(j), E, E_out, mu) - else - call sample_energy(table%nu_d_edist(j), E, E_out) - end if + do + if (law == 44 .or. law == 61) then + call sample_energy(table%nu_d_edist(j), E, E_out, mu) + else + call sample_energy(table%nu_d_edist(j), E, E_out) + end if + ! resample if energy is >= 20 MeV + if (E_out < 20) exit + end do else ! ========================================================== @@ -584,11 +588,15 @@ contains ! sample from prompt neutron energy distribution law = rxn % edist % law - if (law == 44 .or. law == 61) then - call sample_energy(rxn%edist, E, E_out, prob) - else - call sample_energy(rxn%edist, E, E_out) - end if + do + if (law == 44 .or. law == 61) then + call sample_energy(rxn%edist, E, E_out, prob) + else + call sample_energy(rxn%edist, E, E_out) + end if + ! resample if energy is >= 20 MeV + if (E_out < 20) exit + end do end if diff --git a/src/source.f90 b/src/source.f90 index 783fb03d5..80ddbcdd4 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -23,6 +23,7 @@ contains real(8) :: r(3) ! sampled coordinates real(8) :: phi ! azimuthal angle real(8) :: mu ! cosine of polar angle + real(8) :: E ! outgoing energy real(8) :: p_min(3) ! minimum coordinates of source real(8) :: p_max(3) ! maximum coordinates of source character(250) :: msg ! error message @@ -81,7 +82,14 @@ contains p % alive = .true. ! sample energy from Watt fission energy spectrum for U-235 - p % E = watt_spectrum(0.988_8, 2.249_8) + do + E = watt_spectrum(0.988_8, 2.249_8) + ! resample if energy is >= 20 MeV + if (E < 20) exit + end do + + ! set particle energy + p % E = E end do end if end do