diff --git a/src/ace.F90 b/src/ace.F90
index e97338b558..1c49262985 100644
--- a/src/ace.F90
+++ b/src/ace.F90
@@ -1,17 +1,23 @@
module ace
- use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing, &
- DistEnergy
+ use ace_header, only: Nuclide, Reaction, SAlphaBeta, XsListing
use constants
- use endf, only: reaction_name, is_fission, is_disappearance
- use error, only: fatal_error, warning
- use fission, only: nu_total
+ use distribution_univariate, only: Uniform, Equiprobable, Tabular
+ use endf, only: reaction_name, is_fission, is_disappearance
+ use energy_distribution, only: TabularEquiprobable, LevelInelastic, &
+ ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, NBodyPhaseSpace
+ use error, only: fatal_error, warning
+ use fission, only: nu_total
use global
- use list_header, only: ListInt
- use material_header, only: Material
- use output, only: write_message
- use set_header, only: SetChar
- use string, only: to_str, to_lower
+ use list_header, only: ListInt
+ use material_header, only: Material
+ use output, only: write_message
+ use set_header, only: SetChar
+ use secondary_header, only: AngleEnergy
+ use secondary_correlated, only: CorrelatedAngleEnergy
+ use secondary_kalbach, only: KalbachMann
+ use secondary_uncorrelated, only: UncorrelatedAngleEnergy
+ use string, only: to_str, to_lower
implicit none
@@ -372,8 +378,8 @@ contains
else
call read_nu_data(nuc)
call read_reactions(nuc)
- call read_angular_dist(nuc)
call read_energy_dist(nuc)
+ call read_angular_dist(nuc)
call read_unr_res(nuc)
end if
@@ -516,9 +522,10 @@ contains
integer :: LED ! location of energy distribution locators
integer :: LDIS ! location of all energy distributions
integer :: LOCC ! location of energy distributions for given MT
+ integer :: LAW
+ integer :: IDAT
integer :: lc ! locator
integer :: length ! length of data to allocate
- type(DistEnergy), pointer :: edist
JXS2 = JXS(2)
JXS24 = JXS(24)
@@ -661,11 +668,15 @@ contains
! Loop over all delayed neutron precursor groups
do i = 1, NPCR
! find location of energy distribution data
- LOCC = int(XSS(LED + i - 1))
+ LOCC = nint(XSS(LED + i - 1))
+
+ ! Determine law and location of data
+ LAW = nint(XSS(LDIS + LOCC))
+ IDAT = nint(XSS(LDIS + LOCC + 1))
! read energy distribution data
- edist => nuc % nu_d_edist(i)
- call get_energy_dist(edist, LOCC, .true.)
+ call get_energy_dist(nuc%nu_d_edist(i)%obj, LAW, LDIS, IDAT, &
+ ZERO, ZERO)
end do
! =======================================================================
@@ -738,8 +749,8 @@ contains
rxn%multiplicity = 1
rxn%threshold = 1
rxn%scatter_in_cm = .true.
- rxn%has_angle_dist = .false.
- rxn%has_energy_dist = .false.
+ allocate(rxn%secondary%distribution(1))
+ allocate(UncorrelatedAngleEnergy :: rxn%secondary%distribution(1)%obj)
end associate
! Add contribution of elastic scattering to total cross section
@@ -754,10 +765,6 @@ contains
do i = 1, NMT
associate (rxn => nuc % reactions(i+1))
- ! set defaults
- rxn % has_angle_dist = .false.
- rxn % has_energy_dist = .false.
-
! read MT number, Q-value, and neutrons produced
rxn % MT = int(XSS(LMT + i - 1))
rxn % Q_value = XSS(JXS4 + i - 1)
@@ -887,86 +894,96 @@ contains
subroutine read_angular_dist(nuc)
type(Nuclide), intent(inout) :: nuc
- integer :: JXS8 ! location of angular distribution locators
- integer :: JXS9 ! location of angular distributions
integer :: LOCB ! location of angular distribution for given MT
integer :: NE ! number of incoming energies
integer :: NP ! number of points for cosine distribution
- integer :: LC ! locator
integer :: i ! index in reactions array
integer :: j ! index over incoming energies
- integer :: length ! length of data array to allocate
-
- JXS8 = JXS(8)
- JXS9 = JXS(9)
+ integer :: k ! index over energy distributions
+ integer :: interp
+ integer, allocatable :: LC(:) ! locator
! loop over all reactions with secondary neutrons -- NXS(5) does not include
! elastic scattering
do i = 1, NXS(5) + 1
associate (rxn => nuc%reactions(i))
-
! find location of angular distribution
- LOCB = int(XSS(JXS8 + i - 1))
- if (LOCB == -1) then
- ! Angular distribution data are specified through LAWi = 44 in the DLW
- ! block
- cycle
- elseif (LOCB == 0) then
- ! No angular distribution data are given for this reaction, isotropic
- ! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0)
- cycle
- end if
- rxn % has_angle_dist = .true.
+ LOCB = int(XSS(JXS(8) + i - 1))
- ! allocate space for incoming energies and locations
- NE = int(XSS(JXS9 + LOCB - 1))
- rxn % adist % n_energy = NE
- allocate(rxn % adist % energy(NE))
- allocate(rxn % adist % type(NE))
- allocate(rxn % adist % location(NE))
+ ! Angular distribution given as part of a correlated angle-energy distribution
+ if (LOCB == -1) cycle
- ! read incoming energy grid and location of nucs
- XSS_index = JXS9 + LOCB
- rxn % adist % energy = get_real(NE)
- rxn % adist % location = get_int(NE)
+ ! No angular distribution data are given for this reaction, isotropic
+ ! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0)
+ if (LOCB == 0) cycle
- ! determine dize of data block
- length = 0
- do j = 1, NE
- LC = rxn % adist % location(j)
- if (LC == 0) then
- ! isotropic
- rxn % adist % type(j) = ANGLE_ISOTROPIC
- elseif (LC > 0) then
- ! 32 equiprobable bins
- rxn % adist % type(j) = ANGLE_32_EQUI
- length = length + 33
- elseif (LC < 0) then
- ! tabular distribution
- rxn % adist % type(j) = ANGLE_TABULAR
- NP = int(XSS(JXS9 + abs(LC)))
- length = length + 2 + 3*NP
- end if
- end do
+ ! Loop over each separate energy distribution. Even though there is only
+ ! "one" angular distribution, it is repeated as many times as there are
+ ! energy distributions for this reaction since the
+ ! UncorrelatedAngleEnergy type holds one angle and energy distribution.
+ do k = 1, size(rxn%secondary%distribution)
+ select type (aedist => rxn%secondary%distribution(k)%obj)
+ type is (UncorrelatedAngleEnergy)
+ ! allocate space for incoming energies and locations
+ NE = int(XSS(JXS(9) + LOCB - 1))
+ allocate(aedist%angle%energy(NE))
+ allocate(aedist%angle%distribution(NE))
+ allocate(LC(NE))
- ! allocate angular distribution data and read
- allocate(rxn % adist % data(length))
+ ! read incoming energy grid and location of nucs
+ XSS_index = JXS(9) + LOCB
+ aedist%angle%energy(:) = get_real(NE)
+ LC(:) = get_int(NE)
- ! read angular distribution -- currently this does not actually parse the
- ! angular distribution tables for each incoming energy, that must be done
- ! on-the-fly
- XSS_index = JXS9 + LOCB + 2 * NE
- rxn % adist % data = get_real(length)
+ ! determine dize of data block
+ do j = 1, NE
+ if (LC(j) == 0) then
+ ! isotropic
+ allocate(Uniform :: aedist%angle%distribution(j)%obj)
+ select type (adist => aedist%angle%distribution(j)%obj)
+ type is (Uniform)
+ adist%a = -ONE
+ adist%b = ONE
+ end select
- ! change location pointers since they are currently relative to JXS(9)
- LC = LOCB + 2 * NE + 1
- do j = 1, NE
- ! For consistency, leave location as 0 if type is isotropic.
- ! This is not necessary for current correctness, but can avoid
- ! future issues
- if (rxn % adist % location(j) /= 0) then
- rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC
- end if
+ elseif (LC(j) > 0) then
+ ! 32 equiprobable bins
+ allocate(Equiprobable :: aedist%angle%distribution(j)%obj)
+ select type (adist => aedist%angle%distribution(j)%obj)
+ type is (Equiprobable)
+ allocate(adist%x(33))
+ end select
+
+ elseif (LC(j) < 0) then
+ ! tabular distribution
+ allocate(Tabular :: aedist%angle%distribution(j)%obj)
+ end if
+ end do
+
+ ! read angular distribution -- currently this does not actually parse the
+ ! angular distribution tables for each incoming energy, that must be done
+ ! on-the-fly
+ do j = 1, NE
+ XSS_index = JXS(9) + abs(LC(j)) - 1
+ select type(adist => aedist%angle%distribution(j)%obj)
+ type is (Equiprobable)
+ adist%x(:) = get_real(33)
+ type is (Tabular)
+ ! determine interpolation and number of points
+ interp = nint(XSS(XSS_index))
+ NP = nint(XSS(XSS_index + 1))
+
+ ! Get probability density data
+ XSS_index = XSS_index + 2
+ allocate(adist%x(NP), adist%p(NP), adist%c(NP))
+ adist%x(:) = get_real(NP)
+ adist%p(:) = get_real(NP)
+ adist%c(:) = get_real(NP)
+ end select
+ end do
+ deallocate(LC)
+
+ end select
end do
end associate
end do
@@ -981,25 +998,62 @@ contains
subroutine read_energy_dist(nuc)
type(Nuclide), intent(inout) :: nuc
- integer :: LED ! location of energy distribution locators
- integer :: LOCC ! location of energy distributions for given MT
integer :: i ! loop index
-
- LED = JXS(10)
+ integer :: n
+ integer :: IDAT ! locator for distribution data
+ integer :: LNW ! location of next energy law
+ integer :: LAW ! Type of energy law
! Loop over all reactions
do i = 1, NXS(5)
- associate (rxn => nuc % reactions(i+1)) ! skip over elastic scattering
- rxn % has_energy_dist = .true.
+ ! Determine how many energy distributions are present for this reaction
+ LNW = nint(XSS(JXS(10) + i - 1))
+ n = 0
+ do while (LNW > 0)
+ n = n + 1
+ LNW = nint(XSS(JXS(11) + LNW - 1))
+ end do
- ! find location of energy distribution data
- LOCC = int(XSS(LED + i - 1))
+ ! Allocate space for distributions and probability of validity
+ associate (secondary => nuc%reactions(i + 1)%secondary)
+ allocate(secondary%applicability(n))
+ allocate(secondary%distribution(n))
- ! allocate energy distribution
- allocate(rxn % edist)
+ LNW = nint(XSS(JXS(10) + i - 1))
+ n = 0
+ do while (LNW > 0)
+ n = n + 1
- ! read data for energy distribution
- call get_energy_dist(rxn % edist, LOCC)
+ ! Determine energy law and location of data
+ LAW = nint(XSS(JXS(11) + LNW))
+ IDAT = nint(XSS(JXS(11) + LNW + 1))
+
+ ! Read probability of law validity
+ call secondary%applicability(n)%from_ace(XSS, JXS(11) + LNW + 2)
+
+ ! Read energy law data
+ call get_energy_dist(secondary%distribution(n)%obj, LAW, &
+ JXS(11), IDAT, nuc%awr, nuc%reactions(i + 1)%Q_value)
+
+ ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<
+ ! Before the secondary distribution refactor, when the angle/energy
+ ! distribution was uncorrelated, no angle was actually sampled. With
+ ! the refactor, an angle is always sampled for an uncorrelated
+ ! distribution even when no angle distribution exists in the ACE file
+ ! (isotropic is assumed). To preserve the RNG stream, we explicitly
+ ! mark fission reactions so that we avoid the angle sampling.
+ if (any(nuc%reactions(i + 1)%MT == &
+ [N_FISSION, N_F, N_NF, N_2NF, N_3NF])) then
+ select type (aedist => secondary%distribution(n)%obj)
+ type is (UncorrelatedAngleEnergy)
+ aedist%fission = .true.
+ end select
+ end if
+ ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ ! Get locator for next distribution
+ LNW = nint(XSS(JXS(11) + LNW - 1))
+ end do
end associate
end do
@@ -1011,281 +1065,319 @@ contains
! single reaction
!===============================================================================
- recursive subroutine get_energy_dist(edist, loc_law, delayed_n)
- type(DistEnergy), intent(inout) :: edist ! energy distribution
- integer, intent(in) :: loc_law ! locator for data
- logical, intent(in), optional :: delayed_n ! is this for delayed neutrons?
+ recursive subroutine get_energy_dist(aedist, law, LDIS, IDAT, awr, Q_value)
+ class(AngleEnergy), allocatable, intent(inout) :: aedist
+ integer, intent(in) :: law
+ integer, intent(in) :: LDIS
+ integer, intent(in) :: IDAT
+ real(8), intent(in) :: awr
+ real(8), intent(in) :: Q_value
- integer :: LDIS ! location of all energy distributions
- integer :: LNW ! location of next energy distribution if multiple
- integer :: LAW ! secondary energy distribution law
+ integer :: i, j
integer :: NR ! number of interpolation regions
integer :: NE ! number of incoming energies
- integer :: IDAT ! location of first energy distribution for given MT
- integer :: lc ! locator
- integer :: length ! length of data to allocate
- integer :: length_interp_data ! length of interpolation data
+ integer :: NP ! number of outgoing energies/angles
+ integer :: interp
+ integer, allocatable :: L(:) ! locations of distributions for each Ein
+ integer, allocatable :: LC(:) ! locations of distributions for each Ein
- ! determine location of energy distribution
- if (present(delayed_n)) then
- LDIS = JXS(27)
+ XSS_index = LDIS + IDAT - 1
+
+ if (law == 44) then
+ allocate(KalbachMann :: aedist)
+ elseif (law == 61) then
+ allocate(CorrelatedAngleEnergy :: aedist)
else
- LDIS = JXS(11)
+ allocate(UncorrelatedAngleEnergy :: aedist)
end if
- ! locator for next law and information on this law
- LNW = int(XSS(LDIS + loc_law - 1))
- LAW = int(XSS(LDIS + loc_law))
- IDAT = int(XSS(LDIS + loc_law + 1))
- NR = int(XSS(LDIS + loc_law + 2))
- edist % law = LAW
- edist % p_valid % n_regions = NR
+ select type (aedist)
+ type is (UncorrelatedAngleEnergy)
+ ! ========================================================================
+ ! UNCORRELATED ENERGY DISTRIBUTIONS
- ! allocate space for ENDF interpolation parameters
- if (NR > 0) then
- allocate(edist % p_valid % nbt(NR))
- allocate(edist % p_valid % int(NR))
- end if
-
- ! read ENDF interpolation parameters
- XSS_index = LDIS + loc_law + 3
- if (NR > 0) then
- edist % p_valid % nbt = int(get_real(NR))
- edist % p_valid % int = int(get_real(NR))
- end if
-
- ! allocate space for law validity data
- NE = int(XSS(LDIS + loc_law + 3 + 2*NR))
- edist % p_valid % n_pairs = NE
- allocate(edist % p_valid % x(NE))
- allocate(edist % p_valid % y(NE))
-
- length_interp_data = 5 + 2*(NR + NE)
-
- ! read law validity data
- XSS_index = LDIS + loc_law + 4 + 2*NR
- edist % p_valid % x = get_real(NE)
- edist % p_valid % y = get_real(NE)
-
- ! Set index to beginning of IDAT array
- lc = LDIS + IDAT - 2
-
- ! determine length of energy distribution
- length = length_energy_dist(lc, LAW, loc_law, length_interp_data)
-
- ! allocate secondary energy distribution array
- allocate(edist % data(length))
-
- ! read secondary energy distribution
- XSS_index = lc + 1
- edist % data = get_real(length)
-
- ! read next energy distribution if present
- if (LNW > 0) then
- allocate(edist % next)
- call get_energy_dist(edist % next, LNW)
- end if
-
- end subroutine get_energy_dist
-
-!===============================================================================
-! LENGTH_ENERGY_DIST determines how many values are contained in an LDAT energy
-! distribution array based on the secondary energy law and location in XSS
-!===============================================================================
-
- function length_energy_dist(lc, law, LOCC, lid) result(length)
- integer, intent(in) :: lc ! location in XSS array
- integer, intent(in) :: law ! energy distribution law
- integer, intent(in) :: LOCC ! location of energy distribution
- integer, intent(in) :: lid ! length of interpolation data
- integer :: length ! length of energy distribution (LDAT)
-
- integer :: i ! loop index for incoming energies
- integer :: j ! loop index for outgoing energies
- integer :: k ! dummy index in XSS
- integer :: NR ! number of interpolation regions
- integer :: NE ! number of incoming energies
- integer :: NP ! number of points in outgoing energy distribution
- integer :: NMU ! number of points in outgoing cosine distribution
- integer :: NRa ! number of interpolation regions for Watt 'a'
- integer :: NEa ! number of energies for Watt 'a'
- integer :: NRb ! number of interpolation regions for Watt 'b'
- integer :: NEb ! number of energies for Watt 'b'
- real(8), allocatable :: L(:) ! locations of distributions for each Ein
-
- ! initialize length
- length = 0
-
- select case (law)
- case (1)
- ! Tabular equiprobable energy bins
- NR = int(XSS(lc + 1))
- NE = int(XSS(lc + 2 + 2*NR))
- NP = int(XSS(lc + 3 + 2*NR + NE))
- length = 3 + 2*NR + NE + 3*NP*NE
-
- case (2)
- ! Discrete photon energy
- length = 2
-
- case (3)
- ! Level scattering
- length = 2
-
- case (4)
- ! Continuous tabular distribution
- NR = int(XSS(lc + 1))
- NE = int(XSS(lc + 2 + 2*NR))
- allocate(L(NE))
- L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
-
- ! Continue with finding data length
- length = length + 2 + 2*NR + 2*NE
- do i = 1,NE
- ! Some older data sets use the same LDAT for multiple Ein tables.
- ! If this is the case, we should skip incrementing length when it is
- ! not needed.
- if (i < NE) then
- if (any(L(i) == L(i + 1: NE))) then
- ! adjust location for this block
- j = lc + 2 + 2*NR + NE + i
- XSS(j) = XSS(j) - LOCC - lid
- cycle
+ select case (law)
+ case (1)
+ allocate(TabularEquiprobable :: aedist%energy)
+ select type (edist => aedist%energy)
+ type is (TabularEquiprobable)
+ NR = nint(XSS(XSS_index))
+ NE = nint(XSS(XSS_index + 1 + 2*NR))
+ if (NR > 0) then
+ call fatal_error("Multiple interpolation regions not yet supported &
+ &for tabular equiprobable energy distributions.")
end if
- end if
- ! determine length
- NP = int(XSS(lc + length + 2))
- length = length + 2 + 3*NP
+ edist%n_region = NR
- ! adjust location for this block
- j = lc + 2 + 2*NR + NE + i
- XSS(j) = XSS(j) - LOCC - lid
+ ! Read incoming energies for which outgoing energies are tabulated
+ allocate(edist%energy_in(NE))
+ XSS_index = XSS_index + 2 + 2*NR
+ edist%energy_in(:) = get_real(NE)
+
+ ! Read outgoing energy tables
+ NP = nint(XSS(XSS_index))
+ allocate(edist%energy_out(NP, NE))
+ XSS_index = XSS_index + 1
+ do i = 1, NE
+ edist%energy_out(:, i) = get_real(NP)
+ end do
+ end select
+
+ case (3)
+ allocate(LevelInelastic :: aedist%energy)
+ select type (edist => aedist%energy)
+ type is (LevelInelastic)
+ edist%threshold = XSS(XSS_index)
+ edist%mass_ratio = XSS(XSS_index + 1)
+ end select
+
+ case (4)
+ allocate(ContinuousTabular :: aedist%energy)
+ select type (edist => aedist%energy)
+ type is (ContinuousTabular)
+ NR = nint(XSS(XSS_index))
+ XSS_index = XSS_index + 1
+ if (NR > 1) then
+ call fatal_error("Multiple interpolation regions not yet supported &
+ &for continuous tabular energy distributions.")
+ end if
+ edist%n_region = NR
+
+ ! Read breakpoints and interpolation parameters
+ if (NR > 0) then
+ allocate(edist%breakpoints(NR))
+ allocate(edist%interpolation(NR))
+ edist%breakpoints(:) = get_int(NR)
+ edist%interpolation(:) = get_int(NR)
+ end if
+
+ ! Read incoming energies for which outgoing energies are tabulated and
+ ! locators
+ NE = nint(XSS(XSS_index))
+ XSS_index = XSS_index + 1
+ allocate(edist%energy_in(NE))
+ allocate(L(NE))
+ edist%energy_in(:) = get_real(NE)
+ L(:) = get_int(NE)
+
+ ! Read outgoing energy tables
+ allocate(edist%energy_out(NE))
+ do i = 1, NE
+ ! Determine interpolation and number of discrete points
+ XSS_index = LDIS + L(i) - 1
+ interp = nint(XSS(XSS_index))
+ edist%energy_out(i)%interpolation = mod(interp, 10)
+ edist%energy_out(i)%n_discrete = (interp - &
+ edist%energy_out(i)%interpolation)/10
+
+ ! check for discrete lines present
+ if (edist%energy_out(i)%n_discrete > 0) then
+ call fatal_error("Discrete lines in continuous tabular &
+ &distribution not yet supported")
+ end if
+
+ ! Determine number of points and allocate space
+ NP = nint(XSS(XSS_index + 1))
+ allocate(edist%energy_out(i)%e_out(NP))
+ allocate(edist%energy_out(i)%p(NP))
+ allocate(edist%energy_out(i)%c(NP))
+
+ ! Read tabular PDF for outgoing energy
+ XSS_index = XSS_index + 2
+ edist%energy_out(i)%e_out(:) = get_real(NP)
+ edist%energy_out(i)%p(:) = get_real(NP)
+ edist%energy_out(i)%c(:) = get_real(NP)
+ end do
+
+ deallocate(L)
+ end select
+
+ case (7)
+ allocate(MaxwellEnergy :: aedist%energy)
+ select type (edist => aedist%energy)
+ type is (MaxwellEnergy)
+ call edist%theta%from_ace(XSS, XSS_index)
+ edist%u = XSS(XSS_index + 2 + 2*edist%theta%n_regions + &
+ 2*edist%theta%n_pairs)
+ end select
+
+ case (9)
+ allocate(Evaporation :: aedist%energy)
+ select type(edist => aedist%energy)
+ type is (Evaporation)
+ call edist%theta%from_ace(XSS, XSS_index)
+ edist%u = XSS(XSS_index + 2 + 2*edist%theta%n_regions + &
+ 2*edist%theta%n_pairs)
+ end select
+
+ case (11)
+ allocate(WattEnergy :: aedist%energy)
+ select type(edist => aedist%energy)
+ type is (WattEnergy)
+ call edist%a%from_ace(XSS, XSS_index)
+ XSS_index = XSS_index + 2 + 2*edist%a%n_regions + 2*edist%a%n_pairs
+ call edist%b%from_ace(XSS, XSS_index)
+ XSS_index = XSS_index + 2 + 2*edist%b%n_regions + 2*edist%b%n_pairs
+ edist%u = XSS(XSS_index)
+ end select
+
+ case (66)
+ allocate(NBodyPhaseSpace :: aedist%energy)
+ select type(edist => aedist%energy)
+ type is (NBodyPhaseSpace)
+ edist%n_bodies = int(XSS(XSS_index))
+ edist%mass_ratio = XSS(XSS_index + 1)
+ edist%A = awr
+ edist%Q = Q_value
+ end select
+
+ end select
+
+ type is (KalbachMann)
+ ! ========================================================================
+ ! CORRELATED KALBACH-MANN DISTRIBUTION
+
+ NR = int(XSS(XSS_index))
+ NE = int(XSS(XSS_index + 1 + 2*NR))
+ if (NR > 0) then
+ call fatal_error("Multiple interpolation regions not yet supported &
+ &for Kalbach-Mann energy distributions.")
+ end if
+ aedist%n_region = NR
+
+ ! Read incoming energies for which outgoing energies are tabulated and locators
+ allocate(aedist%energy_in(NE))
+ allocate(L(NE))
+ XSS_index = XSS_index + 2 + 2*NR
+ aedist%energy_in(:) = get_real(NE)
+ L(:) = get_int(NE)
+
+ ! Read outgoing energy tables
+ allocate(aedist%table(NE))
+ do i = 1, NE
+ ! Determine interpolation and number of discrete points
+ XSS_index = LDIS + L(i) - 1
+ interp = nint(XSS(XSS_index))
+ aedist%table(i)%interpolation = mod(interp, 10)
+ aedist%table(i)%n_discrete = (interp - aedist%table(i)%interpolation)/10
+
+ ! check for discrete lines present
+ if (aedist%table(i)%n_discrete > 0) then
+ call fatal_error("Discrete lines in Kalbach-Mann distribution not &
+ &yet supported")
+ end if
+
+ ! Determine number of points and allocate space
+ NP = nint(XSS(XSS_index + 1))
+ allocate(aedist%table(i)%e_out(NP))
+ allocate(aedist%table(i)%p(NP))
+ allocate(aedist%table(i)%c(NP))
+ allocate(aedist%table(i)%r(NP))
+ allocate(aedist%table(i)%a(NP))
+
+ ! Read tabular PDF for outgoing energy
+ XSS_index = XSS_index + 2
+ aedist%table(i)%e_out(:) = get_real(NP)
+ aedist%table(i)%p(:) = get_real(NP)
+ aedist%table(i)%c(:) = get_real(NP)
+ aedist%table(i)%r(:) = get_real(NP)
+ aedist%table(i)%a(:) = get_real(NP)
end do
+
deallocate(L)
- case (5)
- ! General evaporation spectrum
- NR = int(XSS(lc + 1))
- NE = int(XSS(lc + 2 + 2*NR))
- NP = int(XSS(lc + 3 + 2*NR + 2*NE))
- length = 3 + 2*NR + 2*NE + NP
+ type is (CorrelatedAngleEnergy)
+ ! ========================================================================
+ ! CORRELATED ANGLE-ENERGY DISTRIBUTION
- case (7)
- ! Maxwell fission spectrum
- NR = int(XSS(lc + 1))
- NE = int(XSS(lc + 2 + 2*NR))
- length = 3 + 2*NR + 2*NE
+ NR = int(XSS(XSS_index))
+ NE = int(XSS(XSS_index + 1 + 2*NR))
+ if (NR > 0) then
+ call fatal_error("Multiple interpolation regions not yet supported &
+ &for correlated angle-energy distributions.")
+ end if
+ aedist%n_region = NR
- case (9)
- ! Evaporation spectrum
- NR = int(XSS(lc + 1))
- NE = int(XSS(lc + 2 + 2*NR))
- length = 3 + 2*NR + 2*NE
-
- case (11)
- ! Watt spectrum
- NRa = int(XSS(lc + 1))
- NEa = int(XSS(lc + 2 + 2*NRa))
- NRb = int(XSS(lc + 3 + 2*(NRa+NEa)))
- NEb = int(XSS(lc + 4 + 2*(NRa+NEa+NRb)))
- length = 5 + 2*(NRa + NEa + NRb + NEb)
-
- case (44)
- ! Kalbach-Mann correlated scattering
- NR = int(XSS(lc + 1))
- NE = int(XSS(lc + 2 + 2*NR))
+ ! Read incoming energies for which outgoing energies are tabulated and
+ ! locators
+ allocate(aedist%energy_in(NE))
allocate(L(NE))
- L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
+ XSS_index = XSS_index + 2 + 2*NR
+ aedist%energy_in(:) = get_real(NE)
+ L(:) = get_int(NE)
- ! Continue with finding data length
- length = length + 2 + 2*NR + 2*NE
- do i = 1,NE
- ! Some older data sets use the same LDAT for multiple Ein tables.
- ! If this is the case, we should skip incrementing length when it is
- ! not needed.
- if (i < NE) then
- if (any(L(i) == L(i + 1: NE))) then
- ! adjust location for this block
- j = lc + 2 + 2*NR + NE + i
- XSS(j) = XSS(j) - LOCC - lid
- cycle
- end if
- end if
- NP = int(XSS(lc + length + 2))
- length = length + 2 + 5*NP
+ ! Read outgoing energy tables
+ allocate(aedist%table(NE))
+ do i = 1, NE
+ ! Determine interpolation and number of discrete points
+ XSS_index = LDIS + L(i) - 1
+ interp = nint(XSS(XSS_index))
+ aedist%table(i)%interpolation = mod(interp, 10)
+ aedist%table(i)%n_discrete = (interp - aedist%table(i)%interpolation)/10
- ! adjust location for this block
- j = lc + 2 + 2*NR + NE + i
- XSS(j) = XSS(j) - LOCC - lid
- end do
- deallocate(L)
-
- case (61)
- ! Correlated energy and angle distribution
- NR = int(XSS(lc + 1))
- NE = int(XSS(lc + 2 + 2*NR))
- allocate(L(NE))
- L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
-
- ! Continue with finding data length
- length = length + 2 + 2*NR + 2*NE
- do i = 1,NE
- ! Some older data sets use the same LDAT for multiple Ein tables.
- ! If this is the case, we should skip incrementing length when it is
- ! not needed.
- if (i < NE) then
- if (any(L(i) == L(i + 1: NE))) then
- ! adjust locators for energy distribution
- j = lc + 2 + 2*NR + NE + i
- XSS(j) = XSS(j) - LOCC - lid
- cycle
- end if
+ ! check for discrete lines present
+ if (aedist%table(i)%n_discrete > 0) then
+ call fatal_error("Discrete lines in correlated angle-energy &
+ &distribution not yet supported")
end if
- ! outgoing energy distribution
- NP = int(XSS(lc + length + 2))
+ ! Determine number of points and allocate space
+ NP = nint(XSS(XSS_index + 1))
+ allocate(aedist%table(i)%e_out(NP))
+ allocate(aedist%table(i)%p(NP))
+ allocate(aedist%table(i)%c(NP))
+ allocate(LC(NP))
- ! adjust locators for angular distribution
+ ! Read tabular PDF for outgoing energy
+ XSS_index = XSS_index + 2
+ aedist%table(i)%e_out(:) = get_real(NP)
+ aedist%table(i)%p(:) = get_real(NP)
+ aedist%table(i)%c(:) = get_real(NP)
+ LC(:) = get_int(NP)
+
+ ! allocate angular distributions for each incoming/outgoing energy
+ allocate(aedist%table(i)%angle(NP))
do j = 1, NP
- k = lc + length + 2 + 3*NP + j
- if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid
+ if (LC(j) == 0) then
+ ! isotropic
+ allocate(Uniform :: aedist%table(i)%angle(j)%obj)
+ select type (adist => aedist%table(i)%angle(j)%obj)
+ type is (Uniform)
+ adist%a = -ONE
+ adist%b = ONE
+ end select
+
+ elseif (LC(j) > 0) then
+ ! tabular distribution
+ allocate(Tabular :: aedist%table(i)%angle(j)%obj)
+ end if
end do
- length = length + 2 + 4*NP
+ ! read angular distributions
do j = 1, NP
- ! outgoing angle distribution -- NMU here is actually
- ! referred to as NP in the MCNP documentation
- NMU = int(XSS(lc + length + 2))
- length = length + 2 + 3*NMU
+ XSS_index = LDIS + abs(LC(j)) - 1
+ select type(adist => aedist%table(i)%angle(j)%obj)
+ type is (Tabular)
+ ! determine interpolation and number of points
+ interp = nint(XSS(XSS_index))
+ NP = nint(XSS(XSS_index + 1))
+
+ ! Get probability density data
+ XSS_index = XSS_index + 2
+ allocate(adist%x(NP), adist%p(NP), adist%c(NP))
+ adist%x(:) = get_real(NP)
+ adist%p(:) = get_real(NP)
+ adist%c(:) = get_real(NP)
+ end select
end do
+ deallocate(LC)
- ! adjust locators for energy distribution
- j = lc + 2 + 2*NR + NE + i
- XSS(j) = XSS(j) - LOCC - lid
end do
- deallocate(L)
- case (66)
- ! N-body phase space distribution
- length = 2
- case (67)
- ! Laboratory energy-angle law
- NR = int(XSS(lc + 1))
- NE = int(XSS(lc + 2 + 2*NR))
- ! Before progressing, check to see if data set uses L(I) values
- ! in a way inconsistent with the current form of the ACE Format Guide
- ! (MCNP5 Manual, Vol 3)
- allocate(L(NE))
- L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
- ! Don't currently do anything with L
deallocate(L)
- ! Continue with finding data length
- NMU = int(XSS(lc + 4 + 2*NR + 2*NE))
- length = 4 + 2*(NR + NE + NMU)
-
end select
- end function length_energy_dist
+ end subroutine get_energy_dist
!===============================================================================
! READ_UNR_RES reads in unresolved resonance probability tables if present.
diff --git a/src/ace_header.F90 b/src/ace_header.F90
index d7c298979c..2331f9117c 100644
--- a/src/ace_header.F90
+++ b/src/ace_header.F90
@@ -3,42 +3,11 @@ module ace_header
use constants, only: MAX_FILE_LEN, ZERO
use dict_header, only: DictIntInt
use endf_header, only: Tab1
+ use secondary_header, only: SecondaryDistribution, AngleEnergyContainer
use stl_vector, only: VectorInt
implicit none
-!===============================================================================
-! DISTANGLE contains data for a tabular secondary angle distribution whether it
-! be tabular or 32 equiprobable cosine bins
-!===============================================================================
-
- type DistAngle
- integer :: n_energy ! # of incoming energies
- real(8), allocatable :: energy(:) ! incoming energy grid
- integer, allocatable :: type(:) ! type of distribution
- integer, allocatable :: location(:) ! location of each table
- real(8), allocatable :: data(:) ! angular distribution data
- end type DistAngle
-
-!===============================================================================
-! DISTENERGY contains data for a secondary energy distribution for all
-! scattering laws
-!===============================================================================
-
- type DistEnergy
- integer :: law ! secondary distribution law
- type(Tab1) :: p_valid ! probability of law validity
- real(8), allocatable :: data(:) ! energy distribution data
-
- ! For reactions that may have multiple energy distributions such as (n,2n),
- ! this pointer allows multiple laws to be stored
- type(DistEnergy), pointer :: next => null()
-
- ! Type-Bound procedures
- contains
- procedure :: clear => distenergy_clear ! Deallocates DistEnergy
- end type DistEnergy
-
!===============================================================================
! REACTION contains the cross-section and secondary energy and angle
! distributions for a single reaction in a continuous-energy ACE-format table
@@ -53,10 +22,7 @@ module ace_header
logical :: scatter_in_cm ! scattering system in center-of-mass?
logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity
real(8), allocatable :: sigma(:) ! Cross section values
- logical :: has_angle_dist ! Angle distribution present?
- logical :: has_energy_dist ! Energy distribution present?
- type(DistAngle) :: adist ! Secondary angular distribution
- type(DistEnergy), pointer :: edist => null() ! Secondary energy distribution
+ type(SecondaryDistribution) :: secondary
! Type-Bound procedures
contains
@@ -137,7 +103,7 @@ module ace_header
integer :: n_precursor ! # of delayed neutron precursors
real(8), allocatable :: nu_d_data(:)
real(8), allocatable :: nu_d_precursor_data(:)
- type(DistEnergy), pointer :: nu_d_edist(:) => null()
+ type(AngleEnergyContainer), allocatable :: nu_d_edist(:)
! Unresolved resonance data
logical :: urr_present
@@ -284,37 +250,14 @@ module ace_header
contains
-!===============================================================================
-! DISTENERGY_CLEAR resets and deallocates data in DistEnergy.
-!===============================================================================
-
- recursive subroutine distenergy_clear(this)
-
- class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear
-
- if (associated(this % next)) then
- ! recursively clear this item
- call this % next % clear()
- deallocate(this % next)
- end if
-
- end subroutine distenergy_clear
-
!===============================================================================
! REACTION_CLEAR resets and deallocates data in Reaction.
!===============================================================================
subroutine reaction_clear(this)
-
class(Reaction), intent(inout) :: this ! The Reaction object to clear
if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E)
-
- if (associated(this % edist)) then
- call this % edist % clear()
- deallocate(this % edist)
- end if
-
end subroutine reaction_clear
!===============================================================================
@@ -322,21 +265,11 @@ module ace_header
!===============================================================================
subroutine nuclide_clear(this)
-
- class(Nuclide), intent(inout) :: this ! The Nuclide object to clear
+ class(Nuclide), intent(inout) :: this
integer :: i ! Loop counter
- if (associated(this % nu_d_edist)) then
- do i = 1, size(this % nu_d_edist)
- call this % nu_d_edist(i) % clear()
- end do
- deallocate(this % nu_d_edist)
- end if
-
- if (associated(this % urr_data)) then
- deallocate(this % urr_data)
- end if
+ if (associated(this % urr_data)) deallocate(this % urr_data)
if (allocated(this % reactions)) then
do i = 1, size(this % reactions)
diff --git a/src/angle_distribution.F90 b/src/angle_distribution.F90
new file mode 100644
index 0000000000..282d51b73b
--- /dev/null
+++ b/src/angle_distribution.F90
@@ -0,0 +1,63 @@
+module angle_distribution
+
+ use constants, only: ZERO, ONE
+ use distribution_univariate, only: DistributionContainer
+ use random_lcg, only: prn
+ use search, only: binary_search
+
+ implicit none
+ private
+
+!===============================================================================
+! ANGLEDISTRIBUTION represents an angular distribution that is to be used in an
+! uncorrelated angle-energy distribution. This occurs whenever the angle
+! distrbution is given in File 4 in an ENDF file. The distribution of angles
+! depends on the incoming energy of the neutron, so this type stores a
+! distribution for each of a set of incoming energies.
+!===============================================================================
+
+ type, public :: AngleDistribution
+ real(8), allocatable :: energy(:)
+ type(DistributionContainer), allocatable :: distribution(:)
+ contains
+ procedure :: sample => angle_sample
+ end type AngleDistribution
+
+contains
+
+ function angle_sample(this, E) result(mu)
+ class(AngleDistribution), intent(in) :: this
+ real(8), intent(in) :: E ! incoming energy
+ real(8) :: mu ! sampled cosine of scattering angle
+
+ integer :: i ! index on incoming energy grid
+ integer :: n ! number of incoming energies
+ real(8) :: r ! interpolation factor on incoming energy grid
+
+ ! Determine number of incoming energies
+ n = size(this%energy)
+
+ ! Find energy bin and calculate interpolation factor -- if the energy is
+ ! outside the range of the tabulated energies, choose the first or last bins
+ if (E < this%energy(1)) then
+ i = 1
+ r = ZERO
+ elseif (E > this%energy(n)) then
+ i = n - 1
+ r = ONE
+ else
+ i = binary_search(this%energy, n, E)
+ r = (E - this%energy(i))/(this%energy(i+1) - this%energy(i))
+ end if
+
+ ! Sample between the ith and (i+1)th bin
+ if (r > prn()) i = i + 1
+
+ ! Sample i-th distribution
+ mu = this%distribution(i)%obj%sample()
+
+ ! Make sure mu is in range [-1,1]
+ if (abs(mu) > ONE) mu = sign(ONE, mu)
+ end function angle_sample
+
+end module angle_distribution
diff --git a/src/distribution_univariate.F90 b/src/distribution_univariate.F90
index b3d77e2f82..f596536d0b 100644
--- a/src/distribution_univariate.F90
+++ b/src/distribution_univariate.F90
@@ -1,7 +1,7 @@
module distribution_univariate
- use constants, only: ZERO, HALF, HISTOGRAM, LINEAR_LINEAR, MAX_LINE_LEN, &
- MAX_WORD_LEN
+ use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, &
+ MAX_LINE_LEN, MAX_WORD_LEN
use error, only: fatal_error
use math, only: maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
@@ -78,6 +78,12 @@ module distribution_univariate
procedure :: initialize => tabular_initialize
end type Tabular
+ type, extends(Distribution) :: Equiprobable
+ real(8), allocatable :: x(:)
+ contains
+ procedure :: sample => equiprobable_sample
+ end type Equiprobable
+
contains
function discrete_sample(this) result(x)
@@ -238,6 +244,25 @@ contains
this%c(:) = this%c(:)/this%c(n)
end subroutine tabular_initialize
+ function equiprobable_sample(this) result(x)
+ class(Equiprobable), intent(in) :: this
+ real(8) :: x
+
+ integer :: i
+ integer :: n
+ real(8) :: r
+ real(8) :: xl, xr
+
+ n = size(this%x)
+
+ r = prn()
+ i = 1 + int((n - 1)*r)
+
+ xl = this%x(i)
+ xr = this%x(i+1)
+ x = xl + ((n - 1)*r - i + ONE) * (xr - xl)
+ end function equiprobable_sample
+
subroutine distribution_from_xml(dist, node_dist)
class(Distribution), allocatable, intent(inout) :: dist
type(Node), pointer :: node_dist
diff --git a/src/endf_header.F90 b/src/endf_header.F90
index af62231a5d..7388ea2f54 100644
--- a/src/endf_header.F90
+++ b/src/endf_header.F90
@@ -13,6 +13,40 @@ module endf_header
integer :: n_pairs ! # of pairs of (x,y) values
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
+ contains
+ procedure :: from_ace
end type Tab1
+contains
+
+ subroutine from_ace(this, xss, idx)
+ class(Tab1), intent(inout) :: this
+ real(8), intent(in) :: xss(:)
+ integer, intent(in) :: idx
+
+ integer :: nr, ne
+
+ ! Determine number of regions
+ nr = nint(xss(idx))
+ this%n_regions = nr
+
+ ! Read interpolation region data
+ if (nr > 0) then
+ allocate(this%nbt(nr))
+ allocate(this%int(nr))
+ this%nbt(:) = nint(xss(idx + 1 : idx + nr))
+ this%int(:) = nint(xss(idx + nr + 1 : idx + 2*nr))
+ end if
+
+ ! Determine number of pairs
+ ne = int(XSS(idx + 2*nr + 1))
+ this%n_pairs = ne
+
+ ! Read (x,y) pairs
+ allocate(this%x(ne))
+ allocate(this%y(ne))
+ this%x(:) = xss(idx + 2*nr + 2 : idx + 2*nr + 1 + ne)
+ this%y(:) = xss(idx + 2*nr + 2 + ne : idx + 2*nr + 1 + 2*ne)
+ end subroutine from_ace
+
end module endf_header
diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90
new file mode 100644
index 0000000000..db3dcc4411
--- /dev/null
+++ b/src/energy_distribution.F90
@@ -0,0 +1,429 @@
+module energy_distribution
+
+ use constants, only: ZERO, ONE, TWO, PI, HISTOGRAM, LINEAR_LINEAR
+ use endf_header, only: Tab1
+ use interpolation, only: interpolate_tab1
+ use math, only: maxwell_spectrum, watt_spectrum
+ use random_lcg, only: prn
+ use search, only: binary_search
+
+!===============================================================================
+! ENERGYDISTRIBUTION (abstract) defines an energy distribution that is a
+! function of the incident energy of a projectile. Each derived type must
+! implement a sample() function that returns a sampled outgoing energy given an
+! incoming energy
+!===============================================================================
+
+ type, abstract :: EnergyDistribution
+ contains
+ procedure(iSampleEnergy), deferred :: sample
+ end type EnergyDistribution
+
+ abstract interface
+ function iSampleEnergy(this, E_in) result(E_out)
+ import EnergyDistribution
+ class(EnergyDistribution), intent(in) :: this
+ real(8), intent(in) :: E_in
+ real(8) :: E_out
+ end function iSampleEnergy
+ end interface
+
+ type :: EnergyDistributionContainer
+ class(EnergyDistribution), allocatable :: obj
+ end type EnergyDistributionContainer
+
+!===============================================================================
+! Derived classes
+!===============================================================================
+
+!===============================================================================
+! TABULAREQUIPROBABLE represents an energy distribution with tabular
+! equiprobable energy bins as given in ACE law 1. This is an older
+! representation that has largely been replaced with ACE laws 4, 44, and 61.
+!===============================================================================
+
+ type, extends(EnergyDistribution) :: TabularEquiprobable
+ integer :: n_region ! number of interpolation regions
+ integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
+ integer, allocatable :: interpolation(:) ! interpolation region codes
+ real(8), allocatable :: energy_in(:) ! incoming energies
+ real(8), allocatable :: energy_out(:,:) ! table of outgoing energies for
+ ! each incoming energy
+ contains
+ procedure :: sample => equiprobable_sample
+ end type TabularEquiprobable
+
+!===============================================================================
+! LEVELINELASTIC gives the energy distribution for level inelastic scattering by
+! neutrons as in ENDF MT=51--90.
+!===============================================================================
+
+ type, extends(EnergyDistribution) :: LevelInelastic
+ real(8) :: threshold
+ real(8) :: mass_ratio
+ contains
+ procedure :: sample => level_inelastic_sample
+ end type LevelInelastic
+
+!===============================================================================
+! CONTINUOUSTABULAR gives an energy distribution represented as a tabular
+! distribution with histogram or linear-linear interpolation. This corresponds
+! to ACE law 4, which NJOY produces for a number of ENDF energy distributions.
+!===============================================================================
+
+ type CTTable
+ integer :: interpolation
+ integer :: n_discrete
+ real(8), allocatable :: e_out(:)
+ real(8), allocatable :: p(:)
+ real(8), allocatable :: c(:)
+ end type CTTable
+
+ type, extends(EnergyDistribution) :: ContinuousTabular
+ integer :: n_region
+ integer, allocatable :: breakpoints(:)
+ integer, allocatable :: interpolation(:)
+ real(8), allocatable :: energy_in(:)
+ type(CTTable), allocatable :: energy_out(:)
+ contains
+ procedure :: sample => continuous_sample
+ end type ContinuousTabular
+
+!===============================================================================
+! MAXWELLENERGY gives the energy distribution of neutrons emitted from a Maxwell
+! fission spectrum. This corresponds to ACE law 7 and ENDF File 5, LF=7.
+!===============================================================================
+
+ type, extends(EnergyDistribution) :: MaxwellEnergy
+ type(Tab1) :: theta ! incoming-energy-dependent parameter
+ real(8) :: u ! restriction energy
+ contains
+ procedure :: sample => maxwellenergy_sample
+ end type MaxwellEnergy
+
+!===============================================================================
+! EVAPORATION represents an evaporation spectrum corresponding to ACE law 9 and
+! ENDF File 5, LF=9.
+!===============================================================================
+
+ type, extends(EnergyDistribution) :: Evaporation
+ type(Tab1) :: theta
+ real(8) :: u
+ contains
+ procedure :: sample => evaporation_sample
+ end type Evaporation
+
+!===============================================================================
+! WATTENERGY gives the energy distribution of neutrons emitted from a Watt
+! fission spectrum. This corresponds to ACE law 11 and ENDF File 5, LF=11.
+!===============================================================================
+
+ type, extends(EnergyDistribution) :: WattEnergy
+ type(Tab1) :: a
+ type(Tab1) :: b
+ real(8) :: u
+ contains
+ procedure :: sample => watt_sample
+ end type WattEnergy
+
+!===============================================================================
+! NBODYPHASESPACE gives the energy distribution for particles emitted from
+! neutron and charged-particle reactions. This corresponds to ACE law 66 and
+! ENDF File 6, LAW=6.
+!===============================================================================
+
+ type, extends(EnergyDistribution) :: NBodyPhaseSpace
+ integer :: n_bodies
+ real(8) :: mass_ratio
+ real(8) :: A
+ real(8) :: Q
+ contains
+ procedure :: sample => nbody_sample
+ end type NBodyPhaseSpace
+
+contains
+
+ function equiprobable_sample(this, E_in) result(E_out)
+ class(TabularEquiprobable), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8) :: E_out ! sampled outgoing energy
+
+ integer :: i, k, l ! indices
+ integer :: n_energy_in ! number of incoming energies
+ integer :: n_energy_out ! number of outgoing energies
+ real(8) :: r ! interpolation factor on incoming energy
+ 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
+ real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
+
+ ! Determine number of incoming/outgoing energies
+ n_energy_in = size(this%energy_in)
+ n_energy_out = size(this%energy_out, 1)
+
+ ! Determine index on incoming energy grid and interpolation factor
+ i = binary_search(this%energy_in, size(this%energy_in), E_in)
+ r = (E_in - this%energy_in(i)) / &
+ (this%energy_in(i+1) - this%energy_in(i))
+
+ ! Sample outgoing energy bin
+ k = 1 + int(n_energy_out * prn())
+
+ ! Determine E_1 and E_K
+ E_i_1 = this%energy_out(1, i)
+ E_i_K = this%energy_out(n_energy_out, i)
+
+ E_i1_1 = this%energy_out(1, i+1)
+ E_i1_K = this%energy_out(n_energy_out, i+1)
+
+ E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
+ E_K = E_i_K + r*(E_i1_K - E_i_K)
+
+ ! Randomly select between the outgoing table for incoming energy E_i and
+ ! E_(i+1)
+ if (prn() < r) then
+ l = i + 1
+ else
+ l = i
+ end if
+
+ ! Determine E_l_k and E_l_k+1
+ E_l_k = this%energy_out(k, l)
+ E_l_k1 = this%energy_out(k+1, l)
+
+ ! Determine E' (denoted here as E_out)
+ E_out = E_l_k + prn()*(E_l_k1 - E_l_k)
+
+ ! Now interpolate between incident energy bins i and i + 1
+ if (l == i) then
+ E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
+ else
+ E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
+ end if
+ end function equiprobable_sample
+
+ function level_inelastic_sample(this, E_in) result(E_out)
+ class(LevelInelastic), intent(in) :: this
+ real(8), intent(in) :: E_in
+ real(8) :: E_out
+
+ E_out = this%mass_ratio*(E_in - this%threshold)
+ end function level_inelastic_sample
+
+ function continuous_sample(this, E_in) result(E_out)
+ class(ContinuousTabular), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8) :: E_out ! sampled outgoing energy
+
+ integer :: i, k, l ! indices
+ integer :: n_energy_in ! number of incoming energies
+ integer :: n_energy_out ! number of outgoing energies
+ real(8) :: r ! interpolation factor on incoming energy
+ real(8) :: r1 ! random number on [0,1)
+ real(8) :: frac ! interpolation factor on outgoing energy
+ 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
+ real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
+ real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
+ real(8) :: c_k, c_k1 ! cumulative probability
+ logical :: histogram_interp ! whether histogram interpolation is used
+
+ ! Read number of interpolation regions and incoming energies
+ if (this%n_region == 1) then
+ histogram_interp = (this%interpolation(1) == 1)
+ else
+ histogram_interp = .false.
+ end if
+
+ ! Find energy bin and calculate interpolation factor -- if the energy is
+ ! outside the range of the tabulated energies, choose the first or last bins
+ n_energy_in = size(this%energy_in)
+ if (E_in < this%energy_in(1)) then
+ i = 1
+ r = ZERO
+ elseif (E_in > this%energy_in(n_energy_in)) then
+ i = n_energy_in - 1
+ r = ONE
+ else
+ i = binary_search(this%energy_in, n_energy_in, E_in)
+ r = (E_in - this%energy_in(i)) / &
+ (this%energy_in(i+1) - this%energy_in(i))
+ end if
+
+ ! Sample between the ith and (i+1)th bin
+ if (histogram_interp) then
+ l = i
+ else
+ if (r > prn()) then
+ l = i + 1
+ else
+ l = i
+ end if
+ end if
+
+ ! Interpolation for energy E1 and EK
+ n_energy_out = size(this%energy_out(i)%e_out)
+ E_i_1 = this%energy_out(i)%e_out(1)
+ E_i_K = this%energy_out(i)%e_out(n_energy_out)
+
+ n_energy_out = size(this%energy_out(i+1)%e_out)
+ E_i1_1 = this%energy_out(i+1)%e_out(1)
+ E_i1_K = this%energy_out(i+1)%e_out(n_energy_out)
+
+ E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
+ E_K = E_i_K + r*(E_i1_K - E_i_K)
+
+ ! Determine outgoing energy bin
+ n_energy_out = size(this%energy_out(l)%e_out)
+ r1 = prn()
+ c_k = this%energy_out(l)%c(1)
+ do k = 1, n_energy_out - 1
+ c_k1 = this%energy_out(l)%c(k+1)
+ if (r1 < c_k1) exit
+ c_k = c_k1
+ end do
+
+ ! Check to make sure k is <= NP - 1
+ k = min(k, n_energy_out - 1)
+
+ E_l_k = this%energy_out(l)%e_out(k)
+ p_l_k = this%energy_out(l)%p(k)
+ if (this%energy_out(l)%interpolation == HISTOGRAM) then
+ ! Histogram interpolation
+ if (p_l_k > ZERO) then
+ E_out = E_l_k + (r1 - c_k)/p_l_k
+ else
+ E_out = E_l_k
+ end if
+
+ elseif (this%energy_out(l)%interpolation == LINEAR_LINEAR) then
+ ! Linear-linear interpolation
+ E_l_k1 = this%energy_out(l)%e_out(k+1)
+ p_l_k1 = this%energy_out(l)%p(k+1)
+
+ frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
+ if (frac == ZERO) then
+ E_out = E_l_k + (r1 - c_k)/p_l_k
+ else
+ E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
+ TWO*frac*(r1 - c_k))) - p_l_k)/frac
+ end if
+ end if
+
+ ! Now interpolate between incident energy bins i and i + 1
+ if (.not. histogram_interp) then
+ if (l == i) then
+ E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
+ else
+ E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
+ end if
+ end if
+ end function continuous_sample
+
+ function maxwellenergy_sample(this, E_in) result(E_out)
+ class(MaxwellEnergy), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8) :: E_out ! sampled outgoing energy
+
+ real(8) :: theta ! Maxwell distribution parameter
+
+ ! Get temperature corresponding to incoming energy
+ theta = interpolate_tab1(this%theta, E_in)
+
+ do
+ ! Sample maxwell fission spectrum
+ E_out = maxwell_spectrum(theta)
+
+ ! Accept energy based on restriction energy
+ if (E_out <= E_in - this%u) exit
+ end do
+ end function maxwellenergy_sample
+
+ function evaporation_sample(this, E_in) result(E_out)
+ class(Evaporation), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8) :: E_out ! sampled outgoing energy
+
+ real(8) :: theta ! evaporation spectrum parameter
+ real(8) :: x, y, v
+
+ ! Get temperature corresponding to incoming energy
+ theta = interpolate_tab1(this%theta, E_in)
+
+ y = (E_in - this%U)/theta
+ v = 1 - exp(-y)
+
+ ! Sample outgoing energy based on evaporation spectrum probability
+ ! density function
+ do
+ x = -log((ONE - v*prn())*(ONE - v*prn()))
+ if (x <= y) exit
+ end do
+
+ E_out = x*theta
+ end function evaporation_sample
+
+ function watt_sample(this, E_in) result(E_out)
+ class(WattEnergy), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8) :: E_out ! sampled outgoing energy
+
+ real(8) :: a, b ! Watt spectrum parameters
+
+ ! Determine Watt parameter 'a' from tabulated function
+ a = interpolate_tab1(this%a, E_in)
+
+ ! Determine Watt parameter 'b' from tabulated function
+ b = interpolate_tab1(this%b, E_in)
+
+ do
+ ! Sample energy-dependent Watt fission spectrum
+ E_out = watt_spectrum(a, b)
+
+ ! Accept energy based on restriction energy
+ if (E_out <= E_in - this%u) exit
+ end do
+ end function watt_sample
+
+ function nbody_sample(this, E_in) result(E_out)
+ class(NBodyPhaseSpace), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8) :: E_out ! sampled outgoing energy
+
+ real(8) :: Ap ! total mass of particles in neutron masses
+ real(8) :: E_max ! maximum possible COM energy
+ real(8) :: x, y, v
+ real(8) :: r1, r2, r3, r4, r5, r6
+
+ ! Determine E_max parameter
+ Ap = this%mass_ratio
+ E_max = (Ap - ONE)/Ap * (this%A/(this%A + ONE)*E_in + this%Q)
+
+ ! x is essentially a Maxwellian distribution
+ x = maxwell_spectrum(ONE)
+
+ select case (this%n_bodies)
+ case (3)
+ y = maxwell_spectrum(ONE)
+ case (4)
+ r1 = prn()
+ r2 = prn()
+ r3 = prn()
+ y = -log(r1*r2*r3)
+ case (5)
+ r1 = prn()
+ r2 = prn()
+ r3 = prn()
+ r4 = prn()
+ r5 = prn()
+ r6 = prn()
+ y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/TWO*r6)**2
+ end select
+
+ ! Now determine v and E_out
+ v = x/(x+y)
+ E_out = E_max * v
+ end function nbody_sample
+
+end module energy_distribution
diff --git a/src/global.F90 b/src/global.F90
index a5943b6b74..4c243ee41f 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -22,7 +22,6 @@ module global
#endif
implicit none
- save
! ============================================================================
! GEOMETRY-RELATED VARIABLES
diff --git a/src/interpolation.F90 b/src/interpolation.F90
index 9e28bc086e..5f87870677 100644
--- a/src/interpolation.F90
+++ b/src/interpolation.F90
@@ -2,7 +2,6 @@ module interpolation
use constants
use endf_header, only: Tab1
- use error, only: fatal_error
use search, only: binary_search
use string, only: to_str
diff --git a/src/output.F90 b/src/output.F90
index 66e6831371..f133bad4bb 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -322,14 +322,8 @@ contains
integer :: i ! loop index over nuclides
integer :: unit_ ! unit to write to
- integer :: size_total ! memory used by nuclide (bytes)
- integer :: size_angle_total ! total memory used for angle dist. (bytes)
- integer :: size_energy_total ! total memory used for energy dist. (bytes)
integer :: size_xs ! memory used for cross-sections (bytes)
- integer :: size_angle ! memory used for an angle distribution (bytes)
- integer :: size_energy ! memory used for a energy distributions (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
- character(11) :: law ! secondary energy distribution law
type(UrrData), pointer :: urr
! set default unit for writing information
@@ -340,8 +334,6 @@ contains
end if
! Initialize totals
- size_angle_total = 0
- size_energy_total = 0
size_urr = 0
size_xs = 0
@@ -356,33 +348,15 @@ contains
write(unit_,*) ' # of reactions = ' // trim(to_str(nuc % n_reaction))
! Information on each reaction
- write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)'
+ write(unit_,*) ' Reaction Q-value COM IE'
do i = 1, nuc % n_reaction
associate (rxn => nuc % reactions(i))
- ! Determine size of angle distribution
- if (rxn % has_angle_dist) then
- size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8
- else
- size_angle = 0
- end if
-
- ! Determine size of energy distribution and law
- if (rxn % has_energy_dist) then
- size_energy = size(rxn % edist % data) * 8
- law = to_str(rxn % edist % law)
- else
- size_energy = 0
- law = 'None'
- end if
-
- write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') &
+ write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
- law(1:4), rxn % threshold, size_angle, size_energy
+ rxn % threshold
! Accumulate data size
size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8
- size_angle_total = size_angle_total + size_angle
- size_energy_total = size_energy_total + size_energy
end associate
end do
@@ -408,19 +382,11 @@ contains
size_urr = urr % n_energy * (urr % n_prob * 6 + 1) * 8
end if
- ! Calculate total memory
- size_total = size_xs + size_angle_total + size_energy_total + size_urr
-
! Write memory used
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
- write(unit_,*) ' Secondary angle distributions = ' // &
- trim(to_str(size_angle_total)) // ' bytes'
- write(unit_,*) ' Secondary energy distributions = ' // &
- trim(to_str(size_energy_total)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
- write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes'
! Blank line at end of nuclide
write(unit_,*)
diff --git a/src/physics.F90 b/src/physics.F90
index 217d232182..31eb981abb 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -1,6 +1,6 @@
module physics
- use ace_header, only: Nuclide, Reaction, DistEnergy
+ use ace_header, only: Nuclide, Reaction
use constants
use cross_section, only: elastic_xs_0K
use endf, only: reaction_name
@@ -17,12 +17,10 @@ module physics
use random_lcg, only: prn
use search, only: binary_search
use string, only: to_str
+ use secondary_uncorrelated, only: UncorrelatedAngleEnergy
implicit none
-! TODO: Figure out how to write particle restart files in sample_angle,
-! sample_energy, etc.
-
contains
!===============================================================================
@@ -458,7 +456,10 @@ contains
vel = sqrt(dot_product(v_n, v_n))
! Sample scattering angle
- mu_cm = sample_angle(rxn, E)
+ select type (dist => rxn%secondary%distribution(1)%obj)
+ type is (UncorrelatedAngleEnergy)
+ mu_cm = dist%angle%sample(E)
+ end select
! Determine direction cosines in CM
uvw_cm = v_n/vel
@@ -1195,16 +1196,13 @@ contains
integer :: lc ! index before start of energies/nu values
integer :: NR ! number of interpolation regions
integer :: NE ! number of energies tabulated
- integer :: law ! energy distribution law
integer :: n_sample ! number of times resampling
real(8) :: nu_t ! total nu
real(8) :: nu_d ! delayed nu
- real(8) :: mu ! fission neutron angular cosine
real(8) :: beta ! delayed neutron fraction
real(8) :: xi ! random number
real(8) :: yield ! delayed neutron precursor yield
real(8) :: prob ! cumulative probability
- type(DistEnergy), pointer :: edist
! Determine total nu
nu_t = nu_total(nuc, p % E)
@@ -1248,18 +1246,13 @@ contains
! set the delayed group for the particle born from fission
p % delayed_group = j
- ! select energy distribution for group j
- law = nuc % nu_d_edist(j) % law
- edist => nuc % nu_d_edist(j)
-
! sample from energy distribution
n_sample = 0
do
- if (law == 44 .or. law == 61) then
- call sample_energy(edist, p % E, E_out, mu)
- else
- call sample_energy(edist, p % E, E_out)
- end if
+ select type (aedist => nuc%nu_d_edist(j)%obj)
+ type is (UncorrelatedAngleEnergy)
+ E_out = aedist%energy%sample(p%E)
+ end select
! resample if energy is greater than maximum neutron energy
if (E_out < energy_max_neutron) exit
@@ -1281,14 +1274,9 @@ contains
p % delayed_group = 0
! sample from prompt neutron energy distribution
- law = rxn % edist % law
n_sample = 0
do
- if (law == 44 .or. law == 61) then
- call sample_energy(rxn%edist, p % E, E_out, prob)
- else
- call sample_energy(rxn%edist, p % E, E_out)
- end if
+ call rxn%secondary%sample(p%E, E_out, prob)
! resample if energy is greater than maximum neutron energy
if (E_out < energy_max_neutron) exit
@@ -1317,41 +1305,18 @@ contains
type(Particle), intent(inout) :: p
integer :: i ! loop index
- integer :: law ! secondary energy distribution law
real(8) :: E ! energy in lab (incoming/outgoing)
real(8) :: mu ! cosine of scattering angle in lab
real(8) :: A ! atomic weight ratio of nuclide
real(8) :: E_in ! incoming energy
real(8) :: E_cm ! outgoing energy in center-of-mass
- real(8) :: Q ! Q-value of reaction
real(8) :: yield ! neutron yield
! copy energy of neutron
E_in = p % E
- ! determine A and Q
- A = nuc % awr
- Q = rxn % Q_value
-
- ! determine secondary energy distribution law
- law = rxn % edist % law
-
- ! sample scattering angle
- mu = sample_angle(rxn, E_in)
-
- ! sample outgoing energy
- if (law == 44 .or. law == 61) then
- call sample_energy(rxn%edist, E_in, E, mu)
- ! Because of floating-point roundoff, it may be possible for mu to be
- ! outside of the range [-1,1). In these cases, we just set mu to exactly
- ! -1 or 1
-
- if (abs(mu) > ONE) mu = sign(ONE,mu)
- elseif (law == 66) then
- call sample_energy(rxn%edist, E_in, E, A=A, Q=Q)
- else
- call sample_energy(rxn%edist, E_in, E)
- end if
+ ! sample outgoing energy and scattering cosine
+ call rxn%secondary%sample(E_in, E, mu)
! if scattering system is in center-of-mass, transfer cosine of scattering
! angle and outgoing energy from CM to LAB
@@ -1359,19 +1324,19 @@ contains
E_cm = E
! determine outgoing energy in lab
+ A = nuc%awr
E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) &
/ ((A+ONE)*(A+ONE))
! determine outgoing angle in lab
mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E)
-
- ! Because of floating-point roundoff, it may be possible for mu to be
- ! outside of the range [-1,1). In these cases, we just set mu to exactly
- ! -1 or 1
-
- if (abs(mu) > ONE) mu = sign(ONE,mu)
end if
+ ! Because of floating-point roundoff, it may be possible for mu to be
+ ! outside of the range [-1,1). In these cases, we just set mu to exactly -1
+ ! or 1
+ if (abs(mu) > ONE) mu = sign(ONE,mu)
+
! Set outgoing energy and scattering angle
p % E = E
p % mu = mu
@@ -1391,893 +1356,4 @@ contains
end subroutine inelastic_scatter
-!===============================================================================
-! SAMPLE_ANGLE samples the cosine of the angle between incident and exiting
-! particle directions either from 32 equiprobable bins or from a tabular
-! distribution.
-!===============================================================================
-
- function sample_angle(rxn, E) result(mu)
- type(Reaction), intent(in) :: rxn ! reaction
- real(8), intent(in) :: E ! incoming energy
-
- real(8) :: xi ! random number on [0,1)
- integer :: interp ! type of interpolation
- integer :: type ! angular distribution type
- integer :: i ! incoming energy bin
- integer :: n ! number of incoming energy bins
- integer :: lc ! location in data array
- integer :: NP ! number of points in cos distribution
- integer :: k ! index on cosine grid
- real(8) :: r ! interpolation factor on incoming energy
- real(8) :: frac ! interpolation fraction on cosine
- real(8) :: mu0 ! cosine in bin k
- real(8) :: mu1 ! cosine in bin k+1
- real(8) :: mu ! final cosine sampled
- real(8) :: c_k ! cumulative frequency at k
- real(8) :: c_k1 ! cumulative frequency at k+1
- real(8) :: p0,p1 ! probability distribution
-
- ! check if reaction has angular distribution -- if not, sample outgoing
- ! angle isotropically
- if (.not. rxn % has_angle_dist) then
- mu = TWO * prn() - ONE
- return
- end if
-
- ! determine number of incoming energies
- n = rxn % adist % n_energy
-
- ! find energy bin and calculate interpolation factor -- if the energy is
- ! outside the range of the tabulated energies, choose the first or last bins
- if (E < rxn % adist % energy(1)) then
- i = 1
- r = ZERO
- elseif (E > rxn % adist % energy(n)) then
- i = n - 1
- r = ONE
- else
- i = binary_search(rxn % adist % energy, n, E)
- r = (E - rxn % adist % energy(i)) / &
- (rxn % adist % energy(i+1) - rxn % adist % energy(i))
- end if
-
- ! Sample between the ith and (i+1)th bin
- if (r > prn()) i = i + 1
-
- ! check whether this is a 32-equiprobable bin or a tabular distribution
- lc = rxn % adist % location(i)
- type = rxn % adist % type(i)
- if (type == ANGLE_ISOTROPIC) then
- mu = TWO * prn() - ONE
- elseif (type == ANGLE_32_EQUI) then
- ! sample cosine bin
- xi = prn()
- k = 1 + int(32.0_8*xi)
-
- ! calculate cosine
- mu0 = rxn % adist % data(lc + k)
- mu1 = rxn % adist % data(lc + k+1)
- mu = mu0 + (32.0_8 * xi - k + ONE) * (mu1 - mu0)
-
- elseif (type == ANGLE_TABULAR) then
- interp = int(rxn % adist % data(lc + 1))
- NP = int(rxn % adist % data(lc + 2))
-
- ! determine outgoing cosine bin
- xi = prn()
- lc = lc + 2
- c_k = rxn % adist % data(lc + 2*NP + 1)
- do k = 1, NP - 1
- c_k1 = rxn % adist % data(lc + 2*NP + k+1)
- if (xi < c_k1) exit
- c_k = c_k1
- end do
-
- ! check to make sure k is <= NP - 1
- k = min(k, NP - 1)
-
- p0 = rxn % adist % data(lc + NP + k)
- mu0 = rxn % adist % data(lc + k)
- if (interp == HISTOGRAM) then
- ! Histogram interpolation
- if (p0 > ZERO) then
- mu = mu0 + (xi - c_k)/p0
- else
- mu = mu0
- end if
-
- elseif (interp == LINEAR_LINEAR) then
- ! Linear-linear interpolation
- p1 = rxn % adist % data(lc + NP + k+1)
- mu1 = rxn % adist % data(lc + k+1)
-
- frac = (p1 - p0)/(mu1 - mu0)
- if (frac == ZERO) then
- mu = mu0 + (xi - c_k)/p0
- else
- mu = mu0 + (sqrt(max(ZERO, p0*p0 + 2*frac*(xi - c_k))) - p0)/frac
- end if
- else
- ! call write_particle_restart(p)
- call fatal_error("Unknown interpolation type: " // trim(to_str(interp)))
- end if
-
- ! Because of floating-point roundoff, it may be possible for mu to be
- ! outside of the range [-1,1). In these cases, we just set mu to exactly
- ! -1 or 1
-
- if (abs(mu) > ONE) mu = sign(ONE,mu)
-
- else
- ! call write_particle_restart(p)
- call fatal_error("Unknown angular distribution type: " &
- &// trim(to_str(type)))
- end if
-
- end function sample_angle
-
-!===============================================================================
-! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary
-! neutron from a collision or for a prompt/delayed fission neutron
-!===============================================================================
-
- recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
- type(DistEnergy), intent(in) :: edist
- real(8), intent(in) :: E_in ! incoming energy of neutron
- real(8), intent(out) :: E_out ! outgoing energy
- real(8), intent(inout), optional :: mu_out ! outgoing cosine of angle
- real(8), intent(in), optional :: A ! mass number of nuclide
- real(8), intent(in), optional :: Q ! Q-value of reaction
-
- integer :: i ! index on incoming energy grid
- integer :: k ! sampled index on outgoing grid
- integer :: l ! sampled index on incoming grid
- integer :: n_sample ! number of rejections
- integer :: lc ! dummy index
- integer :: NR ! number of interpolation regions
- integer :: NE ! number of energies
- integer :: NET ! number of outgoing energies
- integer :: INTTp ! combination of INTT and ND
- integer :: INTT ! 1 = histogram, 2 = linear-linear
- integer :: JJ ! 1 = histogram, 2 = linear-linear
- 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
-
- real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
- real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
- real(8) :: c_k, c_k1 ! cumulative probability
-
- real(8) :: KM_A ! Kalbach-Mann parameter R
- real(8) :: KM_R ! Kalbach-Mann parameter R
- real(8) :: A_k, A_k1 ! Kalbach-Mann A on outgoing grid l
- real(8) :: R_k, R_k1 ! Kalbach-Mann R on outgoing grid l
-
- real(8) :: Watt_a, Watt_b ! Watt spectrum parameters
-
- real(8) :: mu_k ! angular cosine in bin k
- real(8) :: mu_k1 ! angular cosine in bin k+1
- real(8) :: p_k ! angular pdf in bin k
- real(8) :: p_k1 ! angular pdf in bin k+1
-
- real(8) :: r ! interpolation factor on incoming energy
- real(8) :: frac ! interpolation factor on outgoing energy
- real(8) :: U ! restriction energy
- real(8) :: T ! nuclear temperature
-
- real(8) :: Ap ! total mass ratio for n-body dist
- integer :: n_bodies ! number of bodies for n-body dist
- real(8) :: E_max ! parameter for n-body dist
- real(8) :: x, y, v ! intermediate variables for n-body dist
- real(8) :: r1, r2, r3, r4, r5, r6
- logical :: histogram_interp ! use histogram interpolation on incoming energy
-
- ! ==========================================================================
- ! SAMPLE ENERGY DISTRIBUTION IF THERE ARE MULTIPLE
-
- if (associated(edist % next)) then
- p_valid = interpolate_tab1(edist % p_valid, E_in)
-
- if (prn() > 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
-
- ! Determine which secondary energy distribution law to use
- select case (edist % law)
- case (1)
- ! =======================================================================
- ! TABULAR EQUIPROBABLE ENERGY BINS
-
- ! read number of interpolation regions, incoming energies, and outgoing
- ! energies
- NR = int(edist % data(1))
- NE = int(edist % data(2 + 2*NR))
- NET = int(edist % data(3 + 2*NR + NE))
- if (NR > 0) then
- ! call write_particle_restart(p)
- call fatal_error("Multiple interpolation regions not supported while &
- &attempting to sample equiprobable energy bins.")
- end if
-
- ! determine index on incoming energy grid and interpolation factor
- lc = 2 + 2*NR
- i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
- r = (E_in - edist%data(lc+i)) / &
- (edist%data(lc+i+1) - edist%data(lc+i))
-
- ! Sample outgoing energy bin
- r1 = prn()
- k = 1 + int(NET * r1)
-
- ! Determine E_1 and E_K
- lc = 3 + 3*NR + NE + (i-1)*NET
- E_i_1 = edist % data(lc + 1)
- E_i_K = edist % data(lc + NET)
-
- lc = 3 + 3*NR + NE + i*NET
- E_i1_1 = edist % data(lc + 1)
- E_i1_K = edist % data(lc + NET)
-
- E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
- E_K = E_i_K + r*(E_i1_K - E_i_K)
-
- ! Randomly select between the outgoing table for incoming energy E_i and
- ! E_(i+1)
- if (prn() < r) then
- l = i + 1
- else
- l = i
- end if
-
- ! Determine E_l_k and E_l_k+1
- lc = 3 + 2*NR + NE + (l-1)*NET
- E_l_k = edist % data(lc+k)
- E_l_k1 = edist % data(lc+k+1)
-
- ! Determine E' (denoted here as E_out)
- r2 = prn()
- E_out = E_l_k + r2*(E_l_k1 - E_l_k)
-
- ! Now interpolate between incident energy bins i and i + 1
- if (l == i) then
- E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
- else
- E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
- end if
-
- case (3)
- ! =======================================================================
- ! INELASTIC LEVEL SCATTERING
-
- E_out = edist%data(2) * (E_in - edist%data(1))
-
- case (4)
- ! =======================================================================
- ! CONTINUOUS TABULAR DISTRIBUTION
-
- ! read number of interpolation regions and incoming energies
- NR = int(edist % data(1))
- NE = int(edist % data(2 + 2*NR))
- if (NR == 1) then
- histogram_interp = (edist % data(3) == 1)
- else if (NR > 1) then
- ! call write_particle_restart(p)
- call fatal_error("Multiple interpolation regions not supported while &
- &attempting to sample continuous tabular distribution.")
- else
- histogram_interp = .false.
- end if
-
- ! find energy bin and calculate interpolation factor -- if the energy is
- ! outside the range of the tabulated energies, choose the first or last
- ! bins
- lc = 2 + 2*NR
- if (E_in < edist % data(lc+1)) then
- i = 1
- r = ZERO
- elseif (E_in > edist % data(lc+NE)) then
- i = NE - 1
- r = ONE
- else
- i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
- r = (E_in - edist%data(lc+i)) / &
- (edist%data(lc+i+1) - edist%data(lc+i))
- end if
-
- ! Sample between the ith and (i+1)th bin
- if (histogram_interp) then
- l = i
- else
- r2 = prn()
- if (r > r2) then
- l = i + 1
- else
- l = i
- end if
- end if
-
- ! interpolation for energy E1 and EK
- lc = int(edist%data(2 + 2*NR + NE + i))
- NP = int(edist%data(lc + 2))
- E_i_1 = edist%data(lc + 2 + 1)
- E_i_K = edist%data(lc + 2 + NP)
-
- lc = int(edist%data(2 + 2*NR + NE + i + 1))
- NP = int(edist%data(lc + 2))
- E_i1_1 = edist%data(lc + 2 + 1)
- E_i1_K = edist%data(lc + 2 + NP)
-
- E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
- E_K = E_i_K + r*(E_i1_K - E_i_K)
-
- ! determine location of outgoing energies, pdf, cdf for E(l)
- lc = int(edist % data(2 + 2*NR + NE + l))
-
- ! determine type of interpolation and number of discrete lines
- INTTp = int(edist % data(lc + 1))
- NP = int(edist % data(lc + 2))
- if (INTTp > 10) then
- INTT = mod(INTTp,10)
- ND = (INTTp - INTT)/10
- else
- INTT = INTTp
- ND = 0
- end if
-
- if (ND > 0) then
- ! discrete lines present
- ! call write_particle_restart(p)
- call fatal_error("Discrete lines in continuous tabular distributed not &
- &yet supported")
- end if
-
- ! determine outgoing energy bin
- r1 = prn()
- lc = lc + 2 ! start of EOUT
- c_k = edist % data(lc + 2*NP + 1)
- do k = 1, NP - 1
- c_k1 = edist % data(lc + 2*NP + k+1)
- if (r1 < c_k1) exit
- c_k = c_k1
- end do
-
- ! check to make sure k is <= NP - 1
- k = min(k, NP - 1)
-
- E_l_k = edist % data(lc+k)
- p_l_k = edist % data(lc+NP+k)
- if (INTT == HISTOGRAM) then
- ! Histogram interpolation
- if (p_l_k > ZERO) then
- E_out = E_l_k + (r1 - c_k)/p_l_k
- else
- E_out = E_l_k
- end if
-
- elseif (INTT == LINEAR_LINEAR) then
- ! Linear-linear interpolation
- E_l_k1 = edist % data(lc+k+1)
- p_l_k1 = edist % data(lc+NP+k+1)
-
- frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
- if (frac == ZERO) then
- E_out = E_l_k + (r1 - c_k)/p_l_k
- else
- E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
- 2*frac*(r1 - c_k))) - p_l_k)/frac
- end if
- else
- ! call write_particle_restart(p)
- call fatal_error("Unknown interpolation type: " // trim(to_str(INTT)))
- end if
-
- ! Now interpolate between incident energy bins i and i + 1
- if (.not. histogram_interp) then
- if (l == i) then
- E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
- else
- E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
- end if
- end if
-
- case (5)
- ! =======================================================================
- ! GENERAL EVAPORATION SPECTRUM
-
- case (7)
- ! =======================================================================
- ! MAXWELL FISSION SPECTRUM
-
- ! read number of interpolation regions and incoming energies
- NR = int(edist % data(1))
- NE = int(edist % data(2 + 2*NR))
-
- ! determine nuclear temperature from tabulated function
- T = interpolate_tab1(edist % data, E_in)
-
- ! determine restriction energy
- lc = 2 + 2*NR + 2*NE
- U = edist % data(lc + 1)
-
- n_sample = 0
- do
- ! sample maxwell fission spectrum
- E_out = maxwell_spectrum(T)
-
- ! accept energy based on restriction energy
- if (E_out <= E_in - U) exit
-
- ! check for large number of rejections
- n_sample = n_sample + 1
- if (n_sample == MAX_SAMPLE) then
- ! call write_particle_restart(p)
- call fatal_error("Too many rejections on Maxwell fission spectrum.")
- end if
- end do
-
- case (9)
- ! =======================================================================
- ! EVAPORATION SPECTRUM
-
- ! read number of interpolation regions and incoming energies
- NR = int(edist % data(1))
- NE = int(edist % data(2 + 2*NR))
-
- ! determine nuclear temperature from tabulated function
- T = interpolate_tab1(edist % data, E_in)
-
- ! determine restriction energy
- lc = 2 + 2*NR + 2*NE
- U = edist % data(lc + 1)
-
- y = (E_in - U)/T
- v = 1 - exp(-y)
-
- ! sample outgoing energy based on evaporation spectrum probability
- ! density function
- n_sample = 0
- do
- x = -log((1 - v*prn())*(1 - v*prn()))
- if (x <= y) exit
-
- ! check for large number of rejections
- n_sample = n_sample + 1
- if (n_sample == MAX_SAMPLE) then
- ! call write_particle_restart(p)
- call fatal_error("Too many rejections on evaporation spectrum.")
- end if
- end do
-
- E_out = x*T
-
- case (11)
- ! =======================================================================
- ! ENERGY-DEPENDENT WATT SPECTRUM
-
- ! read number of interpolation regions and incoming energies for
- ! parameter 'a'
- NR = int(edist % data(1))
- NE = int(edist % data(2 + 2*NR))
-
- ! determine Watt parameter 'a' from tabulated function
- Watt_a = interpolate_tab1(edist % data, E_in)
-
- ! determine Watt parameter 'b' from tabulated function
- lc = 2 + 2*(NR + NE)
- Watt_b = interpolate_tab1(edist % data, E_in, lc + 1)
-
- ! read number of interpolation regions and incoming energies for
- ! parameter 'a'
- NR = int(edist % data(lc + 1))
- NE = int(edist % data(lc + 2 + 2*NR))
-
- ! determine restriction energy
- lc = lc + 2 + 2*(NR + NE)
- U = edist % data(lc + 1)
-
- n_sample = 0
- do
- ! Sample energy-dependent Watt fission spectrum
- E_out = watt_spectrum(Watt_a, Watt_b)
-
- ! accept energy based on restriction energy
- if (E_out <= E_in - U) exit
-
- ! check for large number of rejections
- n_sample = n_sample + 1
- if (n_sample == MAX_SAMPLE) then
- ! call write_particle_restart(p)
- call fatal_error("Too many rejections on Watt spectrum.")
- end if
- end do
-
- case (44)
- ! =======================================================================
- ! KALBACH-MANN CORRELATED SCATTERING
-
- if (.not. present(mu_out)) then
- ! call write_particle_restart(p)
- call fatal_error("Law 44 called without giving mu_out as argument.")
- end if
-
- ! read number of interpolation regions and incoming energies
- NR = int(edist % data(1))
- NE = int(edist % data(2 + 2*NR))
- if (NR > 0) then
- ! call write_particle_restart(p)
- call fatal_error("Multiple interpolation regions not supported while &
- &attempting to sample Kalbach-Mann distribution.")
- end if
-
- ! find energy bin and calculate interpolation factor -- if the energy is
- ! outside the range of the tabulated energies, choose the first or last
- ! bins
- lc = 2 + 2*NR
- if (E_in < edist % data(lc+1)) then
- i = 1
- r = ZERO
- elseif (E_in > edist % data(lc+NE)) then
- i = NE - 1
- r = ONE
- else
- i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
- r = (E_in - edist%data(lc+i)) / &
- (edist%data(lc+i+1) - edist%data(lc+i))
- end if
-
- ! Sample between the ith and (i+1)th bin
- r2 = prn()
- if (r > r2) then
- l = i + 1
- else
- l = i
- end if
-
- ! determine endpoints on grid i
- lc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i
- NP = int(edist%data(lc + 2))
- E_i_1 = edist%data(lc + 2 + 1)
- E_i_K = edist%data(lc + 2 + NP)
-
- ! determine endpoints on grid i+1
- lc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1
- NP = int(edist%data(lc + 2))
- E_i1_1 = edist%data(lc + 2 + 1)
- E_i1_K = edist%data(lc + 2 + NP)
-
- E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
- E_K = E_i_K + r*(E_i1_K - E_i_K)
-
- ! determine location of outgoing energies, pdf, cdf for E(l)
- lc = int(edist % data(2 + 2*NR + NE + l))
-
- ! determine type of interpolation and number of discrete lines
- INTTp = int(edist % data(lc + 1))
- NP = int(edist % data(lc + 2))
- if (INTTp > 10) then
- INTT = mod(INTTp,10)
- ND = (INTTp - INTT)/10
- else
- INTT = INTTp
- ND = 0
- end if
-
- if (ND > 0) then
- ! discrete lines present
- ! call write_particle_restart(p)
- call fatal_error("Discrete lines in continuous tabular distributed not &
- &yet supported")
- end if
-
- ! determine outgoing energy bin
- r1 = prn()
- lc = lc + 2 ! start of EOUT
- c_k = edist % data(lc + 2*NP + 1)
- do k = 1, NP - 1
- c_k1 = edist % data(lc + 2*NP + k+1)
- if (r1 < c_k1) exit
- c_k = c_k1
- end do
-
- ! check to make sure k is <= NP - 1
- k = min(k, NP - 1)
-
- E_l_k = edist % data(lc+k)
- p_l_k = edist % data(lc+NP+k)
- if (INTT == HISTOGRAM) then
- ! Histogram interpolation
- if (p_l_k > ZERO) then
- E_out = E_l_k + (r1 - c_k)/p_l_k
- else
- E_out = E_l_k
- end if
-
- ! Determine Kalbach-Mann parameters
- KM_R = edist % data(lc + 3*NP + k)
- KM_A = edist % data(lc + 4*NP + k)
-
- elseif (INTT == LINEAR_LINEAR) then
- ! Linear-linear interpolation
- E_l_k1 = edist % data(lc+k+1)
- p_l_k1 = edist % data(lc+NP+k+1)
-
- ! Find E prime
- frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
- if (frac == ZERO) then
- E_out = E_l_k + (r1 - c_k)/p_l_k
- else
- E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
- 2*frac*(r1 - c_k))) - p_l_k)/frac
- end if
-
- ! Determine Kalbach-Mann parameters
- R_k = edist % data(lc + 3*NP + k)
- R_k1 = edist % data(lc + 3*NP + k+1)
- A_k = edist % data(lc + 4*NP + k)
- A_k1 = edist % data(lc + 4*NP + k+1)
-
- KM_R = R_k + (R_k1 - R_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k)
- KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k)
- else
- ! call write_particle_restart()
- call fatal_error("Unknown interpolation type: " // trim(to_str(INTT)))
- end if
-
- ! Now interpolate between incident energy bins i and i + 1
- if (l == i) then
- E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
- else
- E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
- end if
-
- ! Sampled correlated angle from Kalbach-Mann parameters
- r3 = prn()
- r4 = prn()
- if (r3 > KM_R) then
- T = (TWO*r4 - ONE) * sinh(KM_A)
- mu_out = log(T + sqrt(T*T + ONE))/KM_A
- else
- mu_out = log(r4*exp(KM_A) + (ONE - r4)*exp(-KM_A))/KM_A
- end if
-
- case (61)
- ! =======================================================================
- ! CORRELATED ENERGY AND ANGLE DISTRIBUTION
-
- if (.not. present(mu_out)) then
- ! call write_particle_restart()
- call fatal_error("Law 61 called without giving mu_out as argument.")
- end if
-
- ! read number of interpolation regions and incoming energies
- NR = int(edist % data(1))
- NE = int(edist % data(2 + 2*NR))
- if (NR > 0) then
- ! call write_particle_restart()
- call fatal_error("Multiple interpolation regions not supported while &
- &attempting to sample correlated energy-angle distribution.")
- end if
-
- ! find energy bin and calculate interpolation factor -- if the energy is
- ! outside the range of the tabulated energies, choose the first or last
- ! bins
- lc = 2 + 2*NR
- if (E_in < edist % data(lc+1)) then
- i = 1
- r = ZERO
- elseif (E_in > edist % data(lc+NE)) then
- i = NE - 1
- r = ONE
- else
- i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
- r = (E_in - edist%data(lc+i)) / &
- (edist%data(lc+i+1) - edist%data(lc+i))
- end if
-
- ! Sample between the ith and (i+1)th bin
- r2 = prn()
- if (r > r2) then
- l = i + 1
- else
- l = i
- end if
-
- ! determine endpoints on grid i
- lc = int(edist%data(2+2*NR+NE + i)) ! start of LDAT for i
- NP = int(edist%data(lc + 2))
- E_i_1 = edist%data(lc + 2 + 1)
- E_i_K = edist%data(lc + 2 + NP)
-
- ! determine endpoints on grid i+1
- lc = int(edist%data(2+2*NR+NE + i+1)) ! start of LDAT for i+1
- NP = int(edist%data(lc + 2))
- E_i1_1 = edist%data(lc + 2 + 1)
- E_i1_K = edist%data(lc + 2 + NP)
-
- E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
- E_K = E_i_K + r*(E_i1_K - E_i_K)
-
- ! determine location of outgoing energies, pdf, cdf for E(l)
- lc = int(edist % data(2 + 2*NR + NE + l))
-
- ! determine type of interpolation and number of discrete lines
- INTTp = int(edist % data(lc + 1))
- NP = int(edist % data(lc + 2))
- if (INTTp > 10) then
- INTT = mod(INTTp,10)
- ND = (INTTp - INTT)/10
- else
- INTT = INTTp
- ND = 0
- end if
-
- if (ND > 0) then
- ! discrete lines present
- ! call write_particle_restart()
- call fatal_error("Discrete lines in continuous tabular distributed not &
- &yet supported")
- end if
-
- ! determine outgoing energy bin
- r1 = prn()
- lc = lc + 2 ! start of EOUT
- c_k = edist % data(lc + 2*NP + 1)
- do k = 1, NP - 1
- c_k1 = edist % data(lc + 2*NP + k+1)
- if (r1 < c_k1) exit
- c_k = c_k1
- end do
-
- ! check to make sure k is <= NP - 1
- k = min(k, NP - 1)
-
- E_l_k = edist % data(lc+k)
- p_l_k = edist % data(lc+NP+k)
- if (INTT == HISTOGRAM) then
- ! Histogram interpolation
- if (p_l_k > ZERO) then
- E_out = E_l_k + (r1 - c_k)/p_l_k
- else
- E_out = E_l_k
- end if
-
- elseif (INTT == LINEAR_LINEAR) then
- ! Linear-linear interpolation
- E_l_k1 = edist % data(lc+k+1)
- p_l_k1 = edist % data(lc+NP+k+1)
-
- ! Find E prime
- frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
- if (frac == ZERO) then
- E_out = E_l_k + (r1 - c_k)/p_l_k
- else
- E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
- 2*frac*(r1 - c_k))) - p_l_k)/frac
- end if
- else
- ! call write_particle_restart()
- call fatal_error("Unknown interpolation type: " // trim(to_str(INTT)))
- end if
-
- ! Now interpolate between incident energy bins i and i + 1
- if (l == i) then
- E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
- else
- E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
- end if
-
- ! Find correlated angular distribution for closest outgoing energy bin
- if (r1 - c_k < c_k1 - r1) then
- lc = int(edist % data(lc + 3*NP + k))
- else
- lc = int(edist % data(lc + 3*NP + k + 1))
- end if
-
- ! Check if angular distribution is isotropic
- if (lc == 0) then
- mu_out = TWO * prn() - ONE
- return
- end if
-
- ! interpolation type and number of points in angular distribution
- JJ = int(edist % data(lc + 1))
- NP = int(edist % data(lc + 2))
-
- ! determine outgoing cosine bin
- r3 = prn()
- lc = lc + 2
- c_k = edist % data(lc + 2*NP + 1)
- do k = 1, NP - 1
- c_k1 = edist % data(lc + 2*NP + k+1)
- if (r3 < c_k1) exit
- c_k = c_k1
- end do
-
- ! check to make sure k is <= NP - 1
- k = min(k, NP - 1)
-
- p_k = edist % data(lc + NP + k)
- mu_k = edist % data(lc + k)
- if (JJ == HISTOGRAM) then
- ! Histogram interpolation
- if (p_k > ZERO) then
- mu_out = mu_k + (r3 - c_k)/p_k
- else
- mu_out = mu_k
- end if
-
- elseif (JJ == LINEAR_LINEAR) then
- ! Linear-linear interpolation
- p_k1 = edist % data(lc + NP + k+1)
- mu_k1 = edist % data(lc + k+1)
-
- frac = (p_k1 - p_k)/(mu_k1 - mu_k)
- if (frac == ZERO) then
- mu_out = mu_k + (r3 - c_k)/p_k
- else
- mu_out = mu_k + (sqrt(p_k*p_k + 2*frac*(r3 - c_k))-p_k)/frac
- end if
- else
- ! call write_particle_restart()
- call fatal_error("Unknown interpolation type: " // trim(to_str(JJ)))
- end if
-
- case (66)
- ! =======================================================================
- ! N-BODY PHASE SPACE DISTRIBUTION
-
- ! read number of bodies in phase space and total mass ratio
- n_bodies = int(edist % data(1))
- Ap = edist % data(2)
-
- ! determine E_max parameter
- E_max = (Ap - ONE)/Ap * (A/(A+ONE)*E_in + Q)
-
- ! x is essentially a Maxwellian distribution
- x = maxwell_spectrum(ONE)
-
- select case (n_bodies)
- case (3)
- y = maxwell_spectrum(ONE)
- case (4)
- r1 = prn()
- r2 = prn()
- r3 = prn()
- y = -log(r1*r2*r3)
- case (5)
- r1 = prn()
- r2 = prn()
- r3 = prn()
- r4 = prn()
- r5 = prn()
- r6 = prn()
- y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/TWO*r6)**2
- end select
-
- ! now determine v and E_out
- v = x/(x+y)
- E_out = E_max * v
-
- case (67)
- ! =======================================================================
- ! LABORATORY ENERGY-ANGLE LAW
-
- end select
-
- end subroutine sample_energy
-
end module physics
diff --git a/src/search.F90 b/src/search.F90
index 0c345471c2..f338105f4d 100644
--- a/src/search.F90
+++ b/src/search.F90
@@ -1,7 +1,6 @@
module search
use constants
- use error, only: fatal_error
implicit none
diff --git a/src/secondary_correlated.F90 b/src/secondary_correlated.F90
new file mode 100644
index 0000000000..c0289d55eb
--- /dev/null
+++ b/src/secondary_correlated.F90
@@ -0,0 +1,148 @@
+module secondary_correlated
+
+ use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
+ use distribution_univariate, only: DistributionContainer
+ use secondary_header, only: AngleEnergy
+ use random_lcg, only: prn
+ use search, only: binary_search
+
+!===============================================================================
+! CORRELATEDANGLEENERGY represents a correlated angle-energy distribution. This
+! corresponds to ACE law 61 and ENDF File 6, LAW=1, LANG/=2.
+!===============================================================================
+
+ type AngleEnergyTable
+ integer :: interpolation
+ integer :: n_discrete
+ real(8), allocatable :: e_out(:)
+ real(8), allocatable :: p(:)
+ real(8), allocatable :: c(:)
+ type(DistributionContainer), allocatable :: angle(:)
+ end type AngleEnergyTable
+
+ type, extends(AngleEnergy) :: CorrelatedAngleEnergy
+ integer :: n_region ! number of interpolation regions
+ integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
+ integer, allocatable :: interpolation(:) ! interpolation region codes
+ real(8), allocatable :: energy_in(:) ! incoming energies
+ type(AngleEnergyTable), allocatable :: table(:) ! outgoing E/mu distributions
+ contains
+ procedure :: sample => correlated_sample
+ end type CorrelatedAngleEnergy
+
+contains
+
+ subroutine correlated_sample(this, E_in, E_out, mu)
+ class(CorrelatedAngleEnergy), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8), intent(out) :: E_out ! sampled outgoing energy
+ real(8), intent(out) :: mu ! sapmled scattering cosine
+
+ integer :: i, k, l ! indices
+ integer :: n_energy_in ! number of incoming energies
+ integer :: n_energy_out ! number of outgoing energies
+ real(8) :: r ! interpolation factor on incoming energy
+ real(8) :: r1 ! random number on [0,1)
+ real(8) :: frac ! interpolation factor on outgoing energy
+ 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
+ real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
+ real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
+ real(8) :: c_k, c_k1 ! cumulative probability
+
+ ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ ! Before the secondary distribution refactor, an isotropic polar cosine was
+ ! always sampled but then overwritten with the polar cosine sampled from the
+ ! correlated distribution. To preserve the random number stream, we keep
+ ! this dummy sampling here but can remove it later (will change answers)
+ mu = TWO*prn() - ONE
+ ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ ! find energy bin and calculate interpolation factor -- if the energy is
+ ! outside the range of the tabulated energies, choose the first or last bins
+ n_energy_in = size(this%energy_in)
+ if (E_in < this%energy_in(1)) then
+ i = 1
+ r = ZERO
+ elseif (E_in > this%energy_in(n_energy_in)) then
+ i = n_energy_in - 1
+ r = ONE
+ else
+ i = binary_search(this%energy_in, n_energy_in, E_in)
+ r = (E_in - this%energy_in(i)) / &
+ (this%energy_in(i+1) - this%energy_in(i))
+ end if
+
+ ! Sample between the ith and (i+1)th bin
+ if (r > prn()) then
+ l = i + 1
+ else
+ l = i
+ end if
+
+ ! interpolation for energy E1 and EK
+ n_energy_out = size(this%table(i)%e_out)
+ E_i_1 = this%table(i)%e_out(1)
+ E_i_K = this%table(i)%e_out(n_energy_out)
+
+ n_energy_out = size(this%table(i+1)%e_out)
+ E_i1_1 = this%table(i+1)%e_out(1)
+ E_i1_K = this%table(i+1)%e_out(n_energy_out)
+
+ E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
+ E_K = E_i_K + r*(E_i1_K - E_i_K)
+
+ ! determine outgoing energy bin
+ n_energy_out = size(this%table(l)%e_out)
+ r1 = prn()
+ c_k = this%table(l)%c(1)
+ do k = 1, n_energy_out - 1
+ c_k1 = this%table(l)%c(k+1)
+ if (r1 < c_k1) exit
+ c_k = c_k1
+ end do
+
+ ! check to make sure k is <= NP - 1
+ k = min(k, n_energy_out - 1)
+
+ E_l_k = this%table(l)%e_out(k)
+ p_l_k = this%table(l)%p(k)
+ if (this%table(l)%interpolation == HISTOGRAM) then
+ ! Histogram interpolation
+ if (p_l_k > ZERO) then
+ E_out = E_l_k + (r1 - c_k)/p_l_k
+ else
+ E_out = E_l_k
+ end if
+
+ elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
+ ! Linear-linear interpolation
+ E_l_k1 = this%table(l)%e_out(k+1)
+ p_l_k1 = this%table(l)%p(k+1)
+
+ frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
+ if (frac == ZERO) then
+ E_out = E_l_k + (r1 - c_k)/p_l_k
+ else
+ E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
+ TWO*frac*(r1 - c_k))) - p_l_k)/frac
+ end if
+ end if
+
+ ! Now interpolate between incident energy bins i and i + 1
+ if (l == i) then
+ E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
+ else
+ E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
+ end if
+
+ ! Find correlated angular distribution for closest outgoing energy bin
+ if (r1 - c_k < c_k1 - r1) then
+ mu = this%table(l)%angle(k)%obj%sample()
+ else
+ mu = this%table(l)%angle(k + 1)%obj%sample()
+ end if
+ end subroutine correlated_sample
+
+end module secondary_correlated
diff --git a/src/secondary_header.F90 b/src/secondary_header.F90
new file mode 100644
index 0000000000..d858074129
--- /dev/null
+++ b/src/secondary_header.F90
@@ -0,0 +1,78 @@
+module secondary_header
+
+ use endf_header, only: Tab1
+ use interpolation, only: interpolate_tab1
+ use random_lcg, only: prn
+
+!===============================================================================
+! ANGLEENERGY (abstract) defines a correlated or uncorrelated angle-energy
+! distribution that is a function of incoming energy. Each derived type must
+! implement a sample() subroutine that returns an outgoing energy and scattering
+! cosine given an incoming energy.
+!===============================================================================
+
+ type, abstract :: AngleEnergy
+ contains
+ procedure(iSampleAngleEnergy), deferred :: sample
+ end type AngleEnergy
+
+ abstract interface
+ subroutine iSampleAngleEnergy(this, E_in, E_out, mu)
+ import AngleEnergy
+ class(AngleEnergy), intent(in) :: this
+ real(8), intent(in) :: E_in
+ real(8), intent(out) :: E_out
+ real(8), intent(out) :: mu
+ end subroutine iSampleAngleEnergy
+ end interface
+
+ type :: AngleEnergyContainer
+ class(AngleEnergy), allocatable :: obj
+ end type AngleEnergyContainer
+
+!===============================================================================
+! SECONDARYDISTRIBUTION stores multiple angle-energy distributions, each of
+! which has a given probability of occurring for a given incoming energy. In
+! general, most secondary distributions only have one angle-energy distribution,
+! but for some cases (e.g., (n,2n) in certain nuclides) multiple distinct
+! distributions exist.
+!===============================================================================
+
+ type :: SecondaryDistribution
+ type(Tab1), allocatable :: applicability(:)
+ type(AngleEnergyContainer), allocatable :: distribution(:)
+ contains
+ procedure :: sample => secondary_sample
+ end type SecondaryDistribution
+
+contains
+
+ subroutine secondary_sample(this, E_in, E_out, mu)
+ class(SecondaryDistribution), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8), intent(out) :: E_out ! sampled outgoing energy
+ real(8), intent(out) :: mu ! sampled scattering cosine
+
+ integer :: n ! number of angle-energy distributions
+ real(8) :: p_valid ! probability that given distribution is valid
+
+ n = size(this%applicability)
+ if (n > 1) then
+ do i = 1, n
+ ! Determine probability that i-th energy distribution is sampled
+ p_valid = interpolate_tab1(this%applicability(i), E_in)
+
+ ! If i-th distribution is sampled, sample energy from the distribution
+ if (prn() <= p_valid) then
+ call this%distribution(i)%obj%sample(E_in, E_out, mu)
+ exit
+ end if
+ end do
+ else
+ ! If only one distribution is present, go ahead and sample it
+ call this%distribution(1)%obj%sample(E_in, E_out, mu)
+ end if
+
+ end subroutine secondary_sample
+
+end module secondary_header
diff --git a/src/secondary_kalbach.F90 b/src/secondary_kalbach.F90
new file mode 100644
index 0000000000..5e6949206e
--- /dev/null
+++ b/src/secondary_kalbach.F90
@@ -0,0 +1,164 @@
+module secondary_kalbach
+
+ use constants, only: ZERO, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
+ use secondary_header, only: AngleEnergy
+ use random_lcg, only: prn
+ use search, only: binary_search
+
+!===============================================================================
+! KalbachMann represents a correlated angle-energy distribution with the angular
+! distribution represented using Kalbach-Mann systematics. This corresponds to
+! ACE law 44 and ENDF File 6, LAW=1, LANG=2.
+!===============================================================================
+
+ type KalbachMannTable
+ integer :: n_discrete
+ integer :: interpolation
+ real(8), allocatable :: e_out(:)
+ real(8), allocatable :: p(:)
+ real(8), allocatable :: c(:)
+ real(8), allocatable :: r(:)
+ real(8), allocatable :: a(:)
+ end type KalbachMannTable
+
+ type, extends(AngleEnergy) :: KalbachMann
+ integer :: n_region ! number of interpolation regions
+ integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
+ integer, allocatable :: interpolation(:) ! interpolation region codes
+ real(8), allocatable :: energy_in(:) ! incoming energies
+ type(KalbachMannTable), allocatable :: table(:) ! outgoing E/mu parameters
+ contains
+ procedure :: sample => kalbachmann_sample
+ end type KalbachMann
+
+contains
+
+ subroutine kalbachmann_sample(this, E_in, E_out, mu)
+ class(KalbachMann), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8), intent(out) :: E_out ! sampled outgoing energy
+ real(8), intent(out) :: mu ! sampled scattering cosine
+
+ integer :: i, k, l ! indices
+ integer :: n_energy_in ! number of incoming energies
+ integer :: n_energy_out ! number of outgoing energies
+ real(8) :: r ! interpolation factor on incoming energy
+ real(8) :: r1 ! random number on [0,1)
+ real(8) :: frac ! interpolation factor on outgoing energy
+ 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
+ real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
+ real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
+ real(8) :: c_k, c_k1 ! cumulative probability
+ real(8) :: km_r, km_a ! Kalbach-Mann parameters
+ real(8) :: T
+
+ ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ ! Before the secondary distribution refactor, an isotropic polar cosine was
+ ! always sampled but then overwritten with the polar cosine sampled from the
+ ! correlated distribution. To preserve the random number stream, we keep
+ ! this dummy sampling here but can remove it later (will change answers)
+ mu = TWO*prn() - ONE
+ ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ ! find energy bin and calculate interpolation factor -- if the energy is
+ ! outside the range of the tabulated energies, choose the first or last bins
+ n_energy_in = size(this%energy_in)
+ if (E_in < this%energy_in(1)) then
+ i = 1
+ r = ZERO
+ elseif (E_in > this%energy_in(n_energy_in)) then
+ i = n_energy_in - 1
+ r = ONE
+ else
+ i = binary_search(this%energy_in, n_energy_in, E_in)
+ r = (E_in - this%energy_in(i)) / &
+ (this%energy_in(i+1) - this%energy_in(i))
+ end if
+
+ ! Sample between the ith and (i+1)th bin
+ if (r > prn()) then
+ l = i + 1
+ else
+ l = i
+ end if
+
+ ! interpolation for energy E1 and EK
+ n_energy_out = size(this%table(i)%e_out)
+ E_i_1 = this%table(i)%e_out(1)
+ E_i_K = this%table(i)%e_out(n_energy_out)
+
+ n_energy_out = size(this%table(i+1)%e_out)
+ E_i1_1 = this%table(i+1)%e_out(1)
+ E_i1_K = this%table(i+1)%e_out(n_energy_out)
+
+ E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
+ E_K = E_i_K + r*(E_i1_K - E_i_K)
+
+ ! determine outgoing energy bin
+ n_energy_out = size(this%table(l)%e_out)
+ r1 = prn()
+ c_k = this%table(l)%c(1)
+ do k = 1, n_energy_out - 1
+ c_k1 = this%table(l)%c(k+1)
+ if (r1 < c_k1) exit
+ c_k = c_k1
+ end do
+
+ ! check to make sure k is <= NP - 1
+ k = min(k, n_energy_out - 1)
+
+ E_l_k = this%table(l)%e_out(k)
+ p_l_k = this%table(l)%p(k)
+ if (this%table(l)%interpolation == HISTOGRAM) then
+ ! Histogram interpolation
+ if (p_l_k > ZERO) then
+ E_out = E_l_k + (r1 - c_k)/p_l_k
+ else
+ E_out = E_l_k
+ end if
+
+ ! Determine Kalbach-Mann parameters
+ km_r = this%table(l)%r(k)
+ km_a = this%table(l)%a(k)
+
+ elseif (this%table(l)%interpolation == LINEAR_LINEAR) then
+ ! Linear-linear interpolation
+ E_l_k1 = this%table(l)%e_out(k+1)
+ p_l_k1 = this%table(l)%p(k+1)
+
+ frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
+ if (frac == ZERO) then
+ E_out = E_l_k + (r1 - c_k)/p_l_k
+ else
+ E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
+ TWO*frac*(r1 - c_k))) - p_l_k)/frac
+ end if
+
+ ! Determine Kalbach-Mann parameters
+ km_r = this%table(l)%r(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
+ (this%table(l)%r(k+1) - this%table(l)%r(k))
+ km_a = this%table(l)%a(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * &
+ (this%table(l)%a(k+1) - this%table(l)%a(k))
+ end if
+
+ ! Now interpolate between incident energy bins i and i + 1
+ if (l == i) then
+ E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
+ else
+ E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
+ end if
+
+ ! Sampled correlated angle from Kalbach-Mann parameters
+ if (prn() > km_r) then
+ T = (TWO*prn() - ONE) * sinh(km_a)
+ mu = log(T + sqrt(T*T + ONE))/km_a
+ else
+ r1 = prn()
+ mu = log(r1*exp(km_a) + (ONE - r1)*exp(-km_a))/km_a
+ end if
+
+ end subroutine kalbachmann_sample
+
+end module secondary_kalbach
diff --git a/src/secondary_uncorrelated.F90 b/src/secondary_uncorrelated.F90
new file mode 100644
index 0000000000..22a56aa127
--- /dev/null
+++ b/src/secondary_uncorrelated.F90
@@ -0,0 +1,48 @@
+module secondary_uncorrelated
+
+ use angle_distribution, only: AngleDistribution
+ use constants, only: ONE, TWO
+ use energy_distribution, only: EnergyDistribution
+ use secondary_header, only: AngleEnergy
+ use random_lcg, only: prn
+
+!===============================================================================
+! UNCORRELATEDANGLEENERGY represents an uncorrelated angle-energy
+! distribution. This corresponds to when an energy distribution is given in ENDF
+! File 5/6 and an angular distribution is given in ENDF File 4.
+!===============================================================================
+
+ type, extends(AngleEnergy) :: UncorrelatedAngleEnergy
+ logical :: fission = .false.
+ type(AngleDistribution) :: angle
+ class(EnergyDistribution), allocatable :: energy
+ contains
+ procedure :: sample => uncorrelated_sample
+ end type UncorrelatedAngleEnergy
+
+contains
+
+ subroutine uncorrelated_sample(this, E_in, E_out, mu)
+ class(UncorrelatedAngleEnergy), intent(in) :: this
+ real(8), intent(in) :: E_in ! incoming energy
+ real(8), intent(out) :: E_out ! sampled outgoing energy
+ real(8), intent(out) :: mu ! sampled scattering cosine
+
+ ! Sample cosine of scattering angle
+ if (this%fission) then
+ ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ ! For fission, the angle is not used, so just assign a dummy value
+ mu = ONE
+ ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ elseif (allocated(this%angle%energy)) then
+ mu = this%angle%sample(E_in)
+ else
+ ! no angle distribution given => assume isotropic for all energies
+ mu = TWO*prn() - ONE
+ end if
+
+ ! Sample outgoing energy
+ E_out = this%energy%sample(E_in)
+ end subroutine uncorrelated_sample
+
+end module secondary_uncorrelated
diff --git a/tests/test_energy_laws/geometry.xml b/tests/test_energy_laws/geometry.xml
new file mode 100644
index 0000000000..c42f45597d
--- /dev/null
+++ b/tests/test_energy_laws/geometry.xml
@@ -0,0 +1,5 @@
+
+
+
+
+
diff --git a/tests/test_energy_laws/materials.xml b/tests/test_energy_laws/materials.xml
new file mode 100644
index 0000000000..97f9b84be1
--- /dev/null
+++ b/tests/test_energy_laws/materials.xml
@@ -0,0 +1,11 @@
+
+
+ 71c
+
+
+
+
+
+
+
+
diff --git a/tests/test_energy_laws/results_true.dat b/tests/test_energy_laws/results_true.dat
new file mode 100644
index 0000000000..48eb6bc81b
--- /dev/null
+++ b/tests/test_energy_laws/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+2.130076E+00 1.938907E-03
diff --git a/tests/test_energy_laws/settings.xml b/tests/test_energy_laws/settings.xml
new file mode 100644
index 0000000000..1c3f444f0f
--- /dev/null
+++ b/tests/test_energy_laws/settings.xml
@@ -0,0 +1,11 @@
+
+
+
+ 10
+ 5
+ 1000
+
+
+
+
+
diff --git a/tests/test_energy_laws/test_energy_laws.py b/tests/test_energy_laws/test_energy_laws.py
new file mode 100644
index 0000000000..0e593876fc
--- /dev/null
+++ b/tests/test_energy_laws/test_energy_laws.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+
+"""The purpose of this test is to provide coverage of energy distributions that
+are not covered in other tests. It has a single material with the following
+nuclides:
+
+U-233: Only nuclide that has a Watt fission spectrum
+
+H-2: Only nuclide that has an N-body phase space distribution, in this case for
+(n,2n)
+
+Na-23: Has an evaporation spectrum and also has reactions that have multiple
+angle-energy distributions, so it provides coverage for both of those
+situations.
+
+Ta-181: One of a few nuclides that has reactions with Kalbach-Mann distributions
+that use linear-linear interpolation.
+
+"""
+
+import glob
+import os
+import sys
+sys.path.insert(0, os.pardir)
+from testing_harness import TestHarness
+
+
+if __name__ == '__main__':
+ harness = TestHarness('statepoint.10.*')
+ harness.main()