diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index d4167bcb60..49605d4851 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -32,6 +32,7 @@ fileio.o: string.o fission.o: constants.o fission.o: cross_section_header.o fission.o: error.o +fission.o: interpolation.o fission.o: search.o geometry.o: constants.o diff --git a/src/fission.f90 b/src/fission.f90 index 9bc4807e63..cd979562ad 100644 --- a/src/fission.f90 +++ b/src/fission.f90 @@ -3,6 +3,7 @@ module fission use constants use cross_section_header, only: Nuclide use error, only: fatal_error + use interpolation, only: interpolate_tab1 use search, only: binary_search implicit none @@ -21,13 +22,8 @@ contains real(8) :: nu ! number of total neutrons emitted per fission integer :: i ! loop index - integer :: j ! index on nu energy grid / precursor group - integer :: loc ! index before start of energies/nu values integer :: NC ! number of polynomial coefficients - integer :: NR ! number of interpolation regions - integer :: NE ! number of energies tabulated real(8) :: c ! polynomial coefficient - real(8) :: f ! interpolation factor character(MAX_LINE_LEN) :: msg ! error message if (nuc % nu_t_type == NU_NONE) then @@ -44,38 +40,8 @@ contains nu = nu + c * E**i end do elseif (nuc % nu_t_type == NU_TABULAR) then - ! determine number of interpolation regions -- as far as I can tell, no - ! nu data has multiple interpolation regions. Furthermore, it seems all - ! are lin-lin. - NR = int(nuc % nu_t_data(1)) - if (NR /= 0) then - msg = "Multiple interpolation regions not supported while & - &attempting to determine total nu." - call fatal_error(msg) - end if - - ! determine number of energies - loc = 2 + 2*NR - NE = int(nuc % nu_t_data(loc)) - - ! do binary search over tabuled energies to determine appropriate index - ! and interpolation factor - if (E < nuc % nu_t_data(loc+1)) then - j = 1 - f = ZERO - elseif (E > nuc % nu_t_data(loc+NE)) then - j = NE - 1 - f = ONE - else - j = binary_search(nuc % nu_t_data(loc+1), NE, E) - f = (E - nuc % nu_t_data(loc+j)) / & - & (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j)) - end if - - ! determine nu total - loc = loc + NE - nu = nuc % nu_t_data(loc+j) + f * & - & (nuc % nu_t_data(loc+j+1) - nuc % nu_t_data(loc+j)) + ! use ENDF interpolation laws to determine nu + nu = interpolate_tab1(nuc % nu_t_data, E) end if end function nu_total @@ -92,13 +58,8 @@ contains real(8) :: nu ! number of prompt neutrons emitted per fission integer :: i ! loop index - integer :: j ! index on nu energy grid / precursor group - integer :: loc ! index before start of energies/nu values integer :: NC ! number of polynomial coefficients - integer :: NR ! number of interpolation regions - integer :: NE ! number of energies tabulated real(8) :: c ! polynomial coefficient - real(8) :: f ! interpolation factor character(MAX_LINE_LEN) :: msg ! error message if (nuc % nu_p_type == NU_NONE) then @@ -118,36 +79,8 @@ contains nu = nu + c * E**i end do elseif (nuc % nu_p_type == NU_TABULAR) then - ! determine number of interpolation regions - NR = int(nuc % nu_p_data(1)) - if (NR /= 0) then - msg = "Multiple interpolation regions not supported while & - &attempting to determine prompt nu." - call fatal_error(msg) - end if - - ! determine number of energies - loc = 2 + 2*NR - NE = int(nuc % nu_p_data(loc)) - - ! do binary search over tabuled energies to determine appropriate index - ! and interpolation factor - if (E < nuc % nu_p_data(loc+1)) then - j = 1 - f = ZERO - elseif (E > nuc % nu_p_data(loc+NE)) then - j = NE - 1 - f = ONE - else - j = binary_search(nuc % nu_p_data(loc+1), NE, E) - f = (E - nuc % nu_p_data(loc+j)) / & - & (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j)) - end if - - ! determine nu total - loc = loc + NE - nu = nuc % nu_p_data(loc+j) + f * & - & (nuc % nu_p_data(loc+j+1) - nuc % nu_p_data(loc+j)) + ! use ENDF interpolation laws to determine nu + nu = interpolate_tab1(nuc % nu_p_data, E) end if end function nu_prompt @@ -163,46 +96,11 @@ contains real(8), intent(in) :: E ! energy of incoming neutron real(8) :: nu ! number of delayed neutrons emitted per fission - integer :: j ! index on nu energy grid / precursor group - integer :: loc ! index before start of energies/nu values - integer :: NR ! number of interpolation regions - integer :: NE ! number of energies tabulated - real(8) :: f ! interpolation factor - character(MAX_LINE_LEN) :: msg ! error message - if (nuc % nu_d_type == NU_NONE) then nu = ZERO elseif (nuc % nu_d_type == NU_TABULAR) then - ! determine number of interpolation regions - NR = int(nuc % nu_d_data(1)) - if (NR /= 0) then - msg = "Multiple interpolation regions not supported while & - &attempting to determine delayed nu." - call fatal_error(msg) - end if - - ! determine number of energies - loc = 2 + 2*NR - NE = int(nuc % nu_d_data(loc)) - - ! do binary search over tabuled energies to determine appropriate index - ! and interpolation factor - if (E < nuc % nu_d_data(loc+1)) then - j = 1 - f = ZERO - elseif (E > nuc % nu_d_data(loc+NE)) then - j = NE - 1 - f = ONE - else - j = binary_search(nuc % nu_d_data(loc+1), NE, E) - f = (E - nuc % nu_d_data(loc+j)) / & - & (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j)) - end if - - ! determine nu total - loc = loc + NE - nu = nuc % nu_d_data(loc+j) + f * & - & (nuc % nu_d_data(loc+j+1) - nuc % nu_d_data(loc+j)) + ! use ENDF interpolation laws to determine nu + nu = interpolate_tab1(nuc % nu_d_data, E) end if end function nu_delayed