From 56b5e733ac029258df1d1eefc94950498df930f6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 20 Jul 2018 07:36:26 -0500 Subject: [PATCH] Remove unused angle/energy distribution Fortran files --- CMakeLists.txt | 6 - src/angle_distribution.F90 | 130 -------- src/energy_distribution.F90 | 583 --------------------------------- src/nuclide_header.F90 | 2 - src/physics.F90 | 1 - src/product_header.F90 | 153 --------- src/secondary_kalbach.F90 | 278 ---------------- src/secondary_nbody.F90 | 82 ----- src/secondary_uncorrelated.F90 | 98 ------ 9 files changed, 1333 deletions(-) delete mode 100644 src/angle_distribution.F90 delete mode 100644 src/energy_distribution.F90 delete mode 100644 src/product_header.F90 delete mode 100644 src/secondary_kalbach.F90 delete mode 100644 src/secondary_nbody.F90 delete mode 100644 src/secondary_uncorrelated.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index e54c6c3678..a5a760a898 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -288,7 +288,6 @@ set_target_properties(faddeeva PROPERTIES add_library(libopenmc SHARED src/algorithm.F90 - src/angle_distribution.F90 src/angleenergy_header.F90 src/bank_header.F90 src/api.F90 @@ -306,7 +305,6 @@ add_library(libopenmc SHARED src/eigenvalue.F90 src/endf.F90 src/endf_header.F90 - src/energy_distribution.F90 src/error.F90 src/geometry.F90 src/geometry_header.F90 @@ -334,7 +332,6 @@ add_library(libopenmc SHARED src/physics_mg.F90 src/plot.F90 src/plot_header.F90 - src/product_header.F90 src/progress_header.F90 src/pugixml/pugixml_f.F90 src/random_lcg.F90 @@ -342,9 +339,6 @@ add_library(libopenmc SHARED src/relaxng src/sab_header.F90 src/secondary_correlated.F90 - src/secondary_kalbach.F90 - src/secondary_nbody.F90 - src/secondary_uncorrelated.F90 src/set_header.F90 src/settings.F90 src/simulation_header.F90 diff --git a/src/angle_distribution.F90 b/src/angle_distribution.F90 deleted file mode 100644 index a1f20e7799..0000000000 --- a/src/angle_distribution.F90 +++ /dev/null @@ -1,130 +0,0 @@ -module angle_distribution - - use algorithm, only: binary_search - use constants, only: ZERO, ONE, HISTOGRAM, LINEAR_LINEAR - use distribution_univariate, only: DistributionContainer, Tabular - use hdf5_interface, only: read_attribute, get_shape, read_dataset, & - open_dataset, close_dataset, HID_T, HSIZE_T - use random_lcg, only: prn - - 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 - procedure :: from_hdf5 => angle_from_hdf5 - 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 - - subroutine angle_from_hdf5(this, group_id) - class(AngleDistribution), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer :: i, j - integer :: n - integer :: n_energy - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(1), dims2(2) - integer, allocatable :: offsets(:) - integer, allocatable :: interp(:) - real(8), allocatable :: temp(:,:) - - ! Get incoming energies - dset_id = open_dataset(group_id, 'energy') - call get_shape(dset_id, dims) - n_energy = int(dims(1), 4) - allocate(this % energy(n_energy)) - allocate(this % distribution(n_energy)) - call read_dataset(this % energy, dset_id) - call close_dataset(dset_id) - - ! Get outgoing energy distribution data - dset_id = open_dataset(group_id, 'mu') - call read_attribute(offsets, dset_id, 'offsets') - call read_attribute(interp, dset_id, 'interpolation') - call get_shape(dset_id, dims2) - allocate(temp(dims2(1), dims2(2))) - call read_dataset(temp, dset_id) - call close_dataset(dset_id) - - do i = 1, n_energy - ! Determine number of outgoing energies - j = offsets(i) - if (i < n_energy) then - n = offsets(i+1) - j - else - n = size(temp, 1) - j - end if - - ! Create and initialize tabular distribution - allocate(Tabular :: this % distribution(i) % obj) - select type (mudist => this % distribution(i) % obj) - type is (Tabular) - mudist % interpolation = interp(i) - allocate(mudist % x(n), mudist % p(n), mudist % c(n)) - mudist % x(:) = temp(j+1:j+n, 1) - mudist % p(:) = temp(j+1:j+n, 2) - - ! To get answers that match ACE data, for now we still use the tabulated - ! CDF values that were passed through to the HDF5 library. At a later - ! time, we can remove the CDF values from the HDF5 library and - ! reconstruct them using the PDF - if (.true.) then - mudist % c(:) = temp(j+1:j+n, 3) - else - call mudist % initialize(temp(j+1:j+n, 1), temp(j+1:j+n, 2), interp(i)) - end if - end select - - j = j + n - end do - end subroutine angle_from_hdf5 - -end module angle_distribution diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 deleted file mode 100644 index 2be945d80e..0000000000 --- a/src/energy_distribution.F90 +++ /dev/null @@ -1,583 +0,0 @@ -module energy_distribution - - use algorithm, only: binary_search - use constants, only: ZERO, ONE, HALF, TWO, PI, HISTOGRAM, LINEAR_LINEAR - use endf_header, only: Tabulated1D - use hdf5_interface - use math, only: maxwell_spectrum, watt_spectrum - use random_lcg, only: prn - -!=============================================================================== -! 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(energy_distribution_sample_), deferred :: sample - procedure(energy_distribution_from_hdf5_), deferred :: from_hdf5 - end type EnergyDistribution - - abstract interface - function energy_distribution_sample_(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 energy_distribution_sample_ - - subroutine energy_distribution_from_hdf5_(this, group_id) - import EnergyDistribution - import HID_T - class(EnergyDistribution), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - end subroutine energy_distribution_from_hdf5_ - 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 - procedure :: from_hdf5 => equiprobable_from_hdf5 - end type TabularEquiprobable - -!=============================================================================== -! DISCRETEPHOTON gives the energy distribution for a discrete photon (usually -! used for photon production from an incident-neutron reaction) -!=============================================================================== - - type, extends(EnergyDistribution) :: DiscretePhoton - integer :: primary_flag - real(8) :: energy - real(8) :: A - contains - procedure :: sample => discrete_photon_sample - procedure :: from_hdf5 => discrete_photon_from_hdf5 - end type DiscretePhoton - -!=============================================================================== -! 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 - procedure :: from_hdf5 => level_inelastic_from_hdf5 - 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(:) - type(CTTable), allocatable :: distribution(:) - contains - procedure :: sample => continuous_sample - procedure :: from_hdf5 => continuous_from_hdf5 - 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(Tabulated1D) :: theta ! incoming-energy-dependent parameter - real(8) :: u ! restriction energy - contains - procedure :: sample => maxwellenergy_sample - procedure :: from_hdf5 => maxwellenergy_from_hdf5 - end type MaxwellEnergy - -!=============================================================================== -! EVAPORATION represents an evaporation spectrum corresponding to ACE law 9 and -! ENDF File 5, LF=9. -!=============================================================================== - - type, extends(EnergyDistribution) :: Evaporation - type(Tabulated1D) :: theta - real(8) :: u - contains - procedure :: sample => evaporation_sample - procedure :: from_hdf5 => evaporation_from_hdf5 - 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(Tabulated1D) :: a - type(Tabulated1D) :: b - real(8) :: u - contains - procedure :: sample => watt_sample - procedure :: from_hdf5 => watt_from_hdf5 - end type WattEnergy - -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 - - subroutine equiprobable_from_hdf5(this, group_id) - class(TabularEquiprobable), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - end subroutine equiprobable_from_hdf5 - - function discrete_photon_sample(this, E_in) result(E_out) - class(DiscretePhoton), intent(in) :: this - real(8), intent(in) :: E_in - real(8) :: E_out - - if (this % primary_flag == 2) then - E_out = this % energy + this % A/(this % A + 1)*E_in - else - E_out = this % energy - end if - end function discrete_photon_sample - - subroutine discrete_photon_from_hdf5(this, group_id) - class(DiscretePhoton), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - call read_attribute(this % primary_flag, group_id, 'primary_flag') - call read_attribute(this % energy, group_id, 'energy') - call read_attribute(this % A, group_id, 'atomic_weight_ratio') - end subroutine discrete_photon_from_hdf5 - - 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 - - subroutine level_inelastic_from_hdf5(this, group_id) - class(LevelInelastic), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - call read_attribute(this%threshold, group_id, 'threshold') - call read_attribute(this%mass_ratio, group_id, 'mass_ratio') - end subroutine level_inelastic_from_hdf5 - - 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) - if (E_in < this%energy(1)) then - i = 1 - r = ZERO - elseif (E_in > this%energy(n_energy_in)) then - i = n_energy_in - 1 - r = ONE - else - i = binary_search(this%energy, n_energy_in, E_in) - r = (E_in - this%energy(i)) / & - (this%energy(i+1) - this%energy(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%distribution(i)%e_out) - E_i_1 = this%distribution(i)%e_out(1) - E_i_K = this%distribution(i)%e_out(n_energy_out) - - n_energy_out = size(this%distribution(i+1)%e_out) - E_i1_1 = this%distribution(i+1)%e_out(1) - E_i1_K = this%distribution(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%distribution(l)%e_out) - r1 = prn() - c_k = this%distribution(l)%c(1) - do k = 1, n_energy_out - 1 - c_k1 = this%distribution(l)%c(k+1) - if (r1 < c_k1) exit - c_k = c_k1 - end do - - ! Check to make sure 1 <= k <= NP - 1 - k = max(1, min(k, n_energy_out - 1)) - - E_l_k = this%distribution(l)%e_out(k) - p_l_k = this%distribution(l)%p(k) - if (this%distribution(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%distribution(l)%interpolation == LINEAR_LINEAR) then - ! Linear-linear interpolation - E_l_k1 = this%distribution(l)%e_out(k+1) - p_l_k1 = this%distribution(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 .and. n_energy_out > 1) 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 - - subroutine continuous_from_hdf5(this, group_id) - class(ContinuousTabular), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer :: i, j, k - integer :: n - integer :: n_energy - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(1), dims2(2) - integer, allocatable :: temp(:,:) - integer, allocatable :: offsets(:) - integer, allocatable :: interp(:) - integer, allocatable :: n_discrete(:) - real(8), allocatable :: eout(:,:) - - ! Open incoming energy dataset - dset_id = open_dataset(group_id, 'energy') - - ! Get interpolation parameters - call read_attribute(temp, dset_id, 'interpolation') - allocate(this%breakpoints(size(temp, 1))) - allocate(this%interpolation(size(temp, 1))) - this%breakpoints(:) = temp(:, 1) - this%interpolation(:) = temp(:, 2) - this%n_region = size(this%breakpoints) - - ! Get incoming energies - call get_shape(dset_id, dims) - n_energy = int(dims(1), 4) - allocate(this%energy(n_energy)) - allocate(this%distribution(n_energy)) - call read_dataset(this%energy, dset_id) - call close_dataset(dset_id) - - ! Get outgoing energy distribution data - dset_id = open_dataset(group_id, 'distribution') - call read_attribute(offsets, dset_id, 'offsets') - call read_attribute(interp, dset_id, 'interpolation') - call read_attribute(n_discrete, dset_id, 'n_discrete_lines') - call get_shape(dset_id, dims2) - allocate(eout(dims2(1), dims2(2))) - call read_dataset(eout, dset_id) - call close_dataset(dset_id) - - do i = 1, n_energy - ! Determine number of outgoing energies - j = offsets(i) - if (i < n_energy) then - n = offsets(i+1) - j - else - n = size(eout, 1) - j - end if - - associate (d => this % distribution(i)) - ! Assign interpolation scheme and number of discrete lines - d % interpolation = interp(i) - d % n_discrete = n_discrete(i) - - ! Allocate arrays for energies and PDF/CDF - allocate(d % e_out(n)) - allocate(d % p(n)) - allocate(d % c(n)) - - ! Copy data - d % e_out(:) = eout(j+1:j+n, 1) - d % p(:) = eout(j+1:j+n, 2) - - ! To get answers that match ACE data, for now we still use the tabulated - ! CDF values that were passed through to the HDF5 library. At a later - ! time, we can remove the CDF values from the HDF5 library and - ! reconstruct them using the PDF - if (.true.) then - d % c(:) = eout(j+1:j+n, 3) - else - ! Calculate cumulative distribution function -- discrete portion - do k = 1, n_discrete(i) - if (k == 1) then - d % c(k) = d % p(k) - else - d % c(k) = d % c(k-1) + d % p(k) - end if - end do - - ! Continuous portion - do k = d % n_discrete + 1, n - if (k == d % n_discrete + 1) then - d % c(k) = sum(d % p(1:d % n_discrete)) - else - if (d % interpolation == HISTOGRAM) then - d % c(k) = d % c(k-1) + d % p(k-1) * & - (d % e_out(k) - d % e_out(k-1)) - elseif (d % interpolation == LINEAR_LINEAR) then - d % c(k) = d % c(k-1) + HALF*(d % p(k-1) + d % p(k)) * & - (d % e_out(k) - d % e_out(k-1)) - end if - end if - end do - - ! Normalize density and distribution functions - d % p(:) = d % p(:)/d % c(n) - d % c(:) = d % c(:)/d % c(n) - end if - end associate - end do - end subroutine continuous_from_hdf5 - - 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 = this % theta % evaluate(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 - - subroutine maxwellenergy_from_hdf5(this, group_id) - class(MaxwellEnergy), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer(HID_T) :: dset_id - - call read_attribute(this%u, group_id, 'u') - dset_id = open_dataset(group_id, 'theta') - call this%theta%from_hdf5(dset_id) - call close_dataset(dset_id) - end subroutine maxwellenergy_from_hdf5 - - 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 = this % theta % evaluate(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 - - subroutine evaporation_from_hdf5(this, group_id) - class(Evaporation), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer(HID_T) :: dset_id - - call read_attribute(this%u, group_id, 'u') - dset_id = open_dataset(group_id, 'theta') - call this%theta%from_hdf5(dset_id) - call close_dataset(dset_id) - end subroutine evaporation_from_hdf5 - - 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 = this % a % evaluate(E_in) - - ! Determine Watt parameter 'b' from tabulated function - b = this % b % evaluate(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 - - subroutine watt_from_hdf5(this, group_id) - class(WattEnergy), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer(HID_T) :: dset_id - - call read_attribute(this%u, group_id, 'u') - - dset_id = open_dataset(group_id, 'a') - call this%a%from_hdf5(dset_id) - call close_dataset(dset_id) - - dset_id = open_dataset(group_id, 'b') - call this%b%from_hdf5(dset_id) - call close_dataset(dset_id) - end subroutine watt_from_hdf5 - -end module energy_distribution diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 319bf93ab0..8cd225bc98 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -17,11 +17,9 @@ module nuclide_header FIT_T, FIT_A, FIT_F, MultipoleArray use message_passing use multipole_header, only: MultipoleArray - use product_header, only: AngleEnergyContainer use random_lcg, only: prn, future_prn, prn_set_stream use reaction_header, only: Reaction use sab_header, only: SAlphaBeta, sab_tables - use secondary_uncorrelated, only: UncorrelatedAngleEnergy use settings use stl_vector, only: VectorInt, VectorReal use string diff --git a/src/physics.F90 b/src/physics.F90 index dab6bb7f9e..d5f6302a82 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -18,7 +18,6 @@ module physics use random_lcg, only: prn, advance_prn_seed, prn_set_stream use reaction_header, only: Reaction use sab_header, only: sab_tables - use secondary_uncorrelated, only: UncorrelatedAngleEnergy use settings use simulation_header use string, only: to_str diff --git a/src/product_header.F90 b/src/product_header.F90 deleted file mode 100644 index a5e23f4b6c..0000000000 --- a/src/product_header.F90 +++ /dev/null @@ -1,153 +0,0 @@ -module product_header - - use angleenergy_header, only: AngleEnergyContainer - use constants, only: ZERO, MAX_WORD_LEN, EMISSION_PROMPT, EMISSION_DELAYED, & - EMISSION_TOTAL, NEUTRON, PHOTON - use endf_header, only: Tabulated1D, Function1D, Polynomial - use hdf5_interface, only: read_attribute, open_group, close_group, & - open_dataset, close_dataset, read_dataset, HID_T - use random_lcg, only: prn - use secondary_correlated, only: CorrelatedAngleEnergy - use secondary_kalbach, only: KalbachMann - use secondary_nbody, only: NBodyPhaseSpace - use secondary_uncorrelated, only: UncorrelatedAngleEnergy - use string, only: to_str - -!=============================================================================== -! REACTIONPRODUCT stores a data for a reaction product including its yield and -! angle-energy distributions, each of which has a given probability of occurring -! for a given incoming energy. In general, most products only have one -! angle-energy distribution, but for some cases (e.g., (n,2n) in certain -! nuclides) multiple distinct distributions exist. -!=============================================================================== - - type :: ReactionProduct - integer :: particle - integer :: emission_mode ! prompt, delayed, or total emission - real(8) :: decay_rate ! Decay rate for delayed neutron precursors - class(Function1D), pointer :: yield => null() ! Energy-dependent neutron yield - type(Tabulated1D), allocatable :: applicability(:) - type(AngleEnergyContainer), allocatable :: distribution(:) - contains - procedure :: sample => reactionproduct_sample - procedure :: from_hdf5 => reactionproduct_from_hdf5 - end type ReactionProduct - -contains - - subroutine reactionproduct_sample(this, E_in, E_out, mu) - class(ReactionProduct), 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 ! loop counter - integer :: n ! number of angle-energy distributions - real(8) :: prob ! cumulative probability - real(8) :: c ! sampled cumulative probability - - n = size(this%applicability) - if (n > 1) then - prob = ZERO - c = prn() - do i = 1, n - ! Determine probability that i-th energy distribution is sampled - prob = prob + this % applicability(i) % evaluate(E_in) - - ! If i-th distribution is sampled, sample energy from the distribution - if (c <= prob) 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 reactionproduct_sample - - subroutine reactionproduct_from_hdf5(this, group_id) - class(ReactionProduct), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer :: i - integer :: n - integer(HID_T) :: dgroup - integer(HID_T) :: app - integer(HID_T) :: yield - character(MAX_WORD_LEN) :: temp - - ! Read particle type - call read_attribute(temp, group_id, 'particle') - select case (temp) - case ('neutron') - this % particle = NEUTRON - case ('photon') - this % particle = PHOTON - end select - - ! Read emission mode and decay rate - call read_attribute(temp, group_id, 'emission_mode') - select case (temp) - case ('prompt') - this % emission_mode = EMISSION_PROMPT - case ('delayed') - this % emission_mode = EMISSION_DELAYED - case ('total') - this % emission_mode = EMISSION_TOTAL - end select - - ! Read decay rate for delayed emission - if (this % emission_mode == EMISSION_DELAYED) then - call read_attribute(this % decay_rate, group_id, 'decay_rate') - end if - - ! Read secondary particle yield - yield = open_dataset(group_id, 'yield') - call read_attribute(temp, yield, 'type') - select case (temp) - case ('Tabulated1D') - allocate(Tabulated1D :: this % yield) - case ('Polynomial') - allocate(Polynomial :: this % yield) - end select - call this % yield % from_hdf5(yield) - call close_dataset(yield) - - call read_attribute(n, group_id, 'n_distribution') - allocate(this%applicability(n)) - allocate(this%distribution(n)) - - do i = 1, n - dgroup = open_group(group_id, trim('distribution_' // to_str(i - 1))) - - ! Read applicability - if (n > 1) then - app = open_dataset(dgroup, 'applicability') - call this%applicability(i)%from_hdf5(app) - call close_dataset(app) - end if - - ! Read type of distribution and allocate accordingly - call read_attribute(temp, dgroup, 'type') - select case (temp) - case ('uncorrelated') - allocate(UncorrelatedAngleEnergy :: this%distribution(i)%obj) - case ('correlated') - allocate(CorrelatedAngleEnergy :: this%distribution(i)%obj) - case ('nbody') - allocate(NBodyPhaseSpace :: this%distribution(i)%obj) - case ('kalbach-mann') - allocate(KalbachMann :: this%distribution(i)%obj) - end select - - ! Read distribution data - call this%distribution(i)%obj%from_hdf5(dgroup) - - call close_group(dgroup) - end do - - end subroutine reactionproduct_from_hdf5 - -end module product_header diff --git a/src/secondary_kalbach.F90 b/src/secondary_kalbach.F90 deleted file mode 100644 index 3315174e1d..0000000000 --- a/src/secondary_kalbach.F90 +++ /dev/null @@ -1,278 +0,0 @@ -module secondary_kalbach - - use algorithm, only: binary_search - use angleenergy_header, only: AngleEnergy - use constants, only: ZERO, HALF, ONE, TWO, HISTOGRAM, LINEAR_LINEAR - use hdf5_interface - use random_lcg, only: prn - -!=============================================================================== -! 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(:) ! incoming energies - type(KalbachMannTable), allocatable :: distribution(:) ! outgoing E/mu parameters - contains - procedure :: sample => kalbachmann_sample - procedure :: from_hdf5 => kalbachmann_from_hdf5 - 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) - if (E_in < this%energy(1)) then - i = 1 - r = ZERO - elseif (E_in > this%energy(n_energy_in)) then - i = n_energy_in - 1 - r = ONE - else - i = binary_search(this%energy, n_energy_in, E_in) - r = (E_in - this%energy(i)) / & - (this%energy(i+1) - this%energy(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%distribution(i)%e_out) - E_i_1 = this%distribution(i)%e_out(1) - E_i_K = this%distribution(i)%e_out(n_energy_out) - - n_energy_out = size(this%distribution(i+1)%e_out) - E_i1_1 = this%distribution(i+1)%e_out(1) - E_i1_K = this%distribution(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%distribution(l)%e_out) - r1 = prn() - c_k = this%distribution(l)%c(1) - do k = 1, n_energy_out - 1 - c_k1 = this%distribution(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%distribution(l)%e_out(k) - p_l_k = this%distribution(l)%p(k) - if (this%distribution(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%distribution(l)%r(k) - km_a = this%distribution(l)%a(k) - - elseif (this%distribution(l)%interpolation == LINEAR_LINEAR) then - ! Linear-linear interpolation - E_l_k1 = this%distribution(l)%e_out(k+1) - p_l_k1 = this%distribution(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%distribution(l)%r(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * & - (this%distribution(l)%r(k+1) - this%distribution(l)%r(k)) - km_a = this%distribution(l)%a(k) + (E_out - E_l_k)/(E_l_k1 - E_l_k) * & - (this%distribution(l)%a(k+1) - this%distribution(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 - - subroutine kalbachmann_from_hdf5(this, group_id) - class(KalbachMann), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer :: i, j, k - integer :: n - integer :: n_energy - integer(HID_T) :: dset_id - integer(HSIZE_T) :: dims(1), dims2(2) - integer, allocatable :: temp(:,:) - integer, allocatable :: offsets(:) - integer, allocatable :: interp(:) - integer, allocatable :: n_discrete(:) - real(8), allocatable :: eout(:,:) - - ! Open incoming energy dataset - dset_id = open_dataset(group_id, 'energy') - - ! Get interpolation parameters - call read_attribute(temp, dset_id, 'interpolation') - allocate(this%breakpoints(size(temp, 1))) - allocate(this%interpolation(size(temp, 1))) - this%breakpoints(:) = temp(:, 1) - this%interpolation(:) = temp(:, 2) - this%n_region = size(this%breakpoints) - - ! Get incoming energies - call get_shape(dset_id, dims) - n_energy = int(dims(1), 4) - allocate(this%energy(n_energy)) - allocate(this%distribution(n_energy)) - call read_dataset(this%energy, dset_id) - call close_dataset(dset_id) - - ! Get outgoing energy distribution data - dset_id = open_dataset(group_id, 'distribution') - call read_attribute(offsets, dset_id, 'offsets') - call read_attribute(interp, dset_id, 'interpolation') - call read_attribute(n_discrete, dset_id, 'n_discrete_lines') - call get_shape(dset_id, dims2) - allocate(eout(dims2(1), dims2(2))) - call read_dataset(eout, dset_id) - call close_dataset(dset_id) - - do i = 1, n_energy - ! Determine number of outgoing energies - j = offsets(i) - if (i < n_energy) then - n = offsets(i+1) - j - else - n = size(eout, 1) - j - end if - - associate (d => this%distribution(i)) - ! Assign interpolation scheme and number of discrete lines - d % interpolation = interp(i) - d % n_discrete = n_discrete(i) - - ! Allocate arrays for energies and PDF/CDF - allocate(d % e_out(n)) - allocate(d % p(n)) - allocate(d % c(n)) - allocate(d % r(n)) - allocate(d % a(n)) - - ! Copy data - d % e_out(:) = eout(j+1:j+n, 1) - d % p(:) = eout(j+1:j+n, 2) - d % c(:) = eout(j+1:j+n, 3) - d % r(:) = eout(j+1:j+n, 4) - d % a(:) = eout(j+1:j+n, 5) - - - ! To get answers that match ACE data, for now we still use the tabulated - ! CDF values that were passed through to the HDF5 library. At a later - ! time, we can remove the CDF values from the HDF5 library and - ! reconstruct them using the PDF - if (.false.) then - ! Calculate cumulative distribution function -- discrete portion - do k = 1, d % n_discrete - if (k == 1) then - d % c(k) = d % p(k) - else - d % c(k) = d % c(k-1) + d % p(k) - end if - end do - - ! Continuous portion - do k = d % n_discrete + 1, n - if (k == d % n_discrete + 1) then - d % c(k) = sum(d % p(1:d % n_discrete)) - else - if (d % interpolation == HISTOGRAM) then - d % c(k) = d % c(k-1) + d % p(k-1) * & - (d % e_out(k) - d % e_out(k-1)) - elseif (d % interpolation == LINEAR_LINEAR) then - d % c(k) = d % c(k-1) + HALF*(d % p(k-1) + d % p(k)) * & - (d % e_out(k) - d % e_out(k-1)) - end if - end if - end do - - ! Normalize density and distribution functions - d % p(:) = d % p(:)/d % c(n) - d % c(:) = d % c(:)/d % c(n) - end if - end associate - - j = j + n - end do - end subroutine kalbachmann_from_hdf5 - -end module secondary_kalbach diff --git a/src/secondary_nbody.F90 b/src/secondary_nbody.F90 deleted file mode 100644 index aee4b1ff78..0000000000 --- a/src/secondary_nbody.F90 +++ /dev/null @@ -1,82 +0,0 @@ -module secondary_nbody - - use angleenergy_header, only: AngleEnergy - use constants, only: ONE, TWO, PI - use hdf5_interface, only: read_attribute, HID_T - use math, only: maxwell_spectrum - use random_lcg, only: prn - -!=============================================================================== -! 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(AngleEnergy) :: NBodyPhaseSpace - integer :: n_bodies - real(8) :: mass_ratio - real(8) :: A - real(8) :: Q - contains - procedure :: sample => nbody_sample - procedure :: from_hdf5 => nbody_from_hdf5 - end type NBodyPhaseSpace - -contains - - subroutine nbody_sample(this, E_in, E_out, mu) - class(NBodyPhaseSpace), 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 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 - - ! By definition, the distribution of the angle is isotropic for an N-body - ! phase space distribution - mu = TWO*prn() - ONE - - ! 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 subroutine nbody_sample - - subroutine nbody_from_hdf5(this, group_id) - class(NBodyPhaseSpace), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - call read_attribute(this%mass_ratio, group_id, 'total_mass') - call read_attribute(this%n_bodies, group_id, 'n_particles') - call read_attribute(this%A, group_id, 'atomic_weight_ratio') - call read_attribute(this%Q, group_id, 'q_value') - end subroutine nbody_from_hdf5 - -end module secondary_nbody diff --git a/src/secondary_uncorrelated.F90 b/src/secondary_uncorrelated.F90 deleted file mode 100644 index e559c1e5ae..0000000000 --- a/src/secondary_uncorrelated.F90 +++ /dev/null @@ -1,98 +0,0 @@ -module secondary_uncorrelated - - use angle_distribution, only: AngleDistribution - use angleenergy_header, only: AngleEnergy - use constants, only: ONE, TWO, MAX_WORD_LEN - use energy_distribution, only: EnergyDistribution, LevelInelastic, & - ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, DiscretePhoton - use error, only: warning - use hdf5_interface, only: read_attribute, open_group, close_group, & - object_exists, HID_T - 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 - procedure :: from_hdf5 => uncorrelated_from_hdf5 - 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 - - subroutine uncorrelated_from_hdf5(this, group_id) - class(UncorrelatedAngleEnergy), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer(HID_T) :: energy_group - integer(HID_T) :: angle_group - character(MAX_WORD_LEN) :: type - - ! Check if angle group is present & read - if (object_exists(group_id, 'angle')) then - angle_group = open_group(group_id, 'angle') - call this%angle%from_hdf5(angle_group) - call close_group(angle_group) - end if - - ! Check if energy group is present & read - if (object_exists(group_id, 'energy')) then - energy_group = open_group(group_id, 'energy') - call read_attribute(type, energy_group, 'type') - select case (type) - case ('discrete_photon') - allocate(DiscretePhoton :: this%energy) - case ('level') - allocate(LevelInelastic :: this%energy) - case ('continuous') - allocate(ContinuousTabular :: this%energy) - case ('maxwell') - allocate(MaxwellEnergy :: this%energy) - case ('evaporation') - allocate(Evaporation :: this%energy) - case ('watt') - allocate(WattEnergy :: this%energy) - case default - call warning("Energy distribution type '" // trim(type) & - // "' not implemented.") - end select - - if (allocated(this % energy)) then - call this%energy%from_hdf5(energy_group) - end if - - call close_group(energy_group) - end if - end subroutine uncorrelated_from_hdf5 - -end module secondary_uncorrelated