diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 4d5a798a81..0d2bea01c9 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -95,6 +95,7 @@ input_xml.o: xml-fortran/templates/settings_t.o input_xml.o: xml-fortran/templates/tallies_t.o interpolation.o: constants.o +interpolation.o: endf_header.o interpolation.o: error.o interpolation.o: global.o interpolation.o: search.o diff --git a/src/OBJECTS b/src/OBJECTS index db515d0b08..a80d124184 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -6,6 +6,7 @@ datatypes.o \ datatypes_header.o \ doppler.o \ endf.o \ +endf_header.o \ energy_grid.o \ error.o \ fileio.o \ diff --git a/src/endf_header.f90 b/src/endf_header.f90 new file mode 100644 index 0000000000..bee61aa592 --- /dev/null +++ b/src/endf_header.f90 @@ -0,0 +1,18 @@ +module endf_header + + implicit none + +!=============================================================================== +! TAB1 represents a one-dimensional interpolable function +!=============================================================================== + + type Tab1 + integer :: n_regions ! # of interpolation regions + integer, allocatable :: nbt(:) ! values separating interpolation regions + integer, allocatable :: int(:) ! interpolation scheme + integer :: n_pairs ! # of pairs of (x,y) values + real(8), allocatable :: x(:) ! values of abscissa + real(8), allocatable :: y(:) ! values of ordinate + end type Tab1 + +end module endf_header diff --git a/src/interpolation.f90 b/src/interpolation.f90 index 231007e738..e9d1f5337a 100644 --- a/src/interpolation.f90 +++ b/src/interpolation.f90 @@ -1,24 +1,29 @@ module interpolation use constants - use error, only: fatal_error - use global, only: message - use search, only: binary_search + use endf_header, only: Tab1 + use error, only: fatal_error + use global, only: message + use search, only: binary_search implicit none + interface interpolate_tab1 + module procedure interpolate_tab1_array, interpolate_tab1_object + end interface interpolate_tab1 + contains !=============================================================================== -! INTERPOLATE_TAB1 interpolates a function between two points based on +! INTERPOLATE_TAB1_ARRAY interpolates a function between two points based on ! particular interpolation scheme. The data needs to be organized as a ENDF TAB1 ! type function containing the interpolation regions, break points, and ! tabulated x's and y's. !=============================================================================== - function interpolate_tab1(data, x, loc_start) result(y) + function interpolate_tab1_array(data, x, loc_start) result(y) - real(8), intent(in) :: data(:) ! array of data whose + real(8), intent(in) :: data(:) ! array of data real(8), intent(in) :: x ! x value to find y at integer, intent(in), optional :: loc_start ! starting location in data real(8) :: y ! y(x) @@ -108,6 +113,84 @@ contains y = exp((1-r)*log(y0) + r*log(y1)) end select - end function interpolate_tab1 + end function interpolate_tab1_array + +!=============================================================================== +! INTERPOLATE_TAB1_OBJECT interpolates a function between two points based on +! particular interpolation scheme. The data needs to be organized as a ENDF TAB1 +! type function containing the interpolation regions, break points, and +! tabulated x's and y's. +!=============================================================================== + + function interpolate_tab1_object(obj, x) result(y) + + type(Tab1), intent(in) :: obj ! ENDF Tab1 interpolable function + real(8), intent(in) :: x ! x value to find y at + real(8) :: y ! y(x) + + integer :: i ! bin in which to interpolate + integer :: n_regions ! number of interpolation regions + integer :: n_pairs ! number of tabulated values + integer :: interp ! ENDF interpolation scheme + real(8) :: r ! interpolation factor + real(8) :: x0, x1 ! bounding x values + real(8) :: y0, y1 ! bounding y values + + ! determine number of interpolation regions and pairs + n_regions = obj % n_regions + n_pairs = obj % n_pairs + + ! find which bin the abscissa is in -- if the abscissa is outside the + ! tabulated range, the first or last point is chosen, i.e. no interpolation + ! is done outside the energy range + if (x < obj % x(1)) then + y = obj % y(1) + return + elseif (x > obj % x(n_pairs)) then + y = obj % y(n_pairs) + return + else + i = binary_search(obj % x, n_pairs, x) + end if + + ! determine interpolation scheme + if (n_regions == 0) then + interp = LINEAR_LINEAR + elseif (n_regions == 1) then + interp = obj % int(1) + elseif (n_regions > 1) then + message = "Multiple interpolation regions not yet supported." + call fatal_error() + end if + + ! handle special case of histogram interpolation + if (interp == HISTOGRAM) then + y = obj % y(i) + return + end if + + ! determine bounding values + x0 = obj % x(i) + x1 = obj % x(i + 1) + y0 = obj % y(i) + y1 = obj % y(i + 1) + + ! determine interpolation factor and interpolated value + select case (interp) + case (LINEAR_LINEAR) + r = (x - x0)/(x1 - x0) + y = (1 - r)*y0 + r*y1 + case (LINEAR_LOG) + r = (log(x) - log(x0))/(log(x1) - log(x0)) + y = (1 - r)*y0 + r*y1 + case (LOG_LINEAR) + r = (x - x0)/(x1 - x0) + y = exp((1-r)*log(y0) + r*log(y1)) + case (LOG_LOG) + r = (log(x) - log(x0))/(log(x1) - log(x0)) + y = exp((1-r)*log(y0) + r*log(y1)) + end select + + end function interpolate_tab1_object end module interpolation