From d4237668d20b88e9cc213cfa14844cd9f1f7246f Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 31 May 2018 20:52:54 -0500 Subject: [PATCH] Allow TTB implementation to handle any photon cutoff energy --- src/input_xml.F90 | 7 ++++++ src/material_header.F90 | 2 +- src/photon_header.F90 | 51 +++++++++++++++++++++++++++++++++++++---- src/photon_physics.F90 | 34 +++++++++++++-------------- 4 files changed, 72 insertions(+), 22 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index fa077549c..4d29cdc81 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4388,6 +4388,13 @@ contains if (allocated(ttb_k_grid)) deallocate(ttb_k_grid) end do + ! Determine if minimum/maximum energy for bremsstrahlung is greater/less + ! than the current minimum/maximum + if (size(ttb_e_grid) >= 1) then + energy_min(PHOTON) = max(energy_min(PHOTON), ttb_e_grid(1)) + energy_max(PHOTON) = min(energy_max(PHOTON), ttb_e_grid(size(ttb_e_grid))) + end if + ! Take logarithm of energies since they are log-log interpolated ttb_e_grid = log(ttb_e_grid) end if diff --git a/src/material_header.F90 b/src/material_header.F90 index 460103b65..6cdff7397 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -872,7 +872,7 @@ contains do j = 2, n_e ! Set last element of PDF to small non-zero value to enable log-log ! interpolation - this % pdf(j,j) = 1.0e-6_8 * this % pdf(j-1,j) + this % pdf(j,j) = exp(-500.0_8) ! Loop over photon energies c = ZERO diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 35b506a7f..bb48c7f40 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -131,7 +131,11 @@ contains character(3), allocatable :: designators(:) real(8) :: a real(8) :: c + real(8) :: f + real(8) :: y + real(8), allocatable :: electron_energy(:) real(8), allocatable :: matrix(:,:) + real(8), allocatable :: dcs(:,:) ! Get name of nuclide from group name_len = len(this % name) @@ -340,10 +344,8 @@ contains call close_dataset(dset_id) ! Get energy grids used for bremsstrahlung DCS and for stopping powers - if (.not. allocated(ttb_e_grid)) then - allocate(ttb_e_grid(n_e)) - call read_dataset(ttb_e_grid, rgroup, 'electron_energy') - end if + allocate(electron_energy(n_e)) + call read_dataset(electron_energy, rgroup, 'electron_energy') if (.not. allocated(ttb_k_grid)) then allocate(ttb_k_grid(n_k)) call read_dataset(ttb_k_grid, rgroup, 'photon_energy') @@ -360,6 +362,47 @@ contains call read_attribute(this % I, rgroup, 'I') call close_group(rgroup) end if + + ! Truncate the bremsstrahlung data at the cutoff energy + if (energy_cutoff(PHOTON) > electron_energy(1)) then + i_grid = binary_search(electron_energy, n_e, energy_cutoff(PHOTON)) + + ! calculate interpolation factor + f = (log(energy_cutoff(PHOTON)) - log(electron_energy(i_grid))) / & + (log(electron_energy(i_grid+1)) - log(electron_energy(i_grid))) + + ! Interpolate collision stopping power at the cutoff energy and + ! truncate + y = exp(log(this % stopping_power_collision(i_grid)) + & + f*(log(this % stopping_power_collision(i_grid+1)) - & + log(this % stopping_power_collision(i_grid)))) + this % stopping_power_collision = & + [y, this % stopping_power_collision(i_grid+1:n_e)] + + ! Interpolate radiative stopping power at the cutoff energy and + ! truncate + y = exp(log(this % stopping_power_radiative(i_grid)) + & + f*(log(this % stopping_power_radiative(i_grid+1)) - & + log(this % stopping_power_radiative(i_grid)))) + this % stopping_power_radiative = & + [y, this % stopping_power_radiative(i_grid+1:n_e)] + + ! Interpolate bremsstrahlung DCS at the cutoff energy and truncate + allocate(dcs(n_k, n_e-i_grid+1)) + do i = 1, n_k + y = exp(log(this % dcs(i,i_grid)) + & + f*(log(this % dcs(i,i_grid+1)) - log(this % dcs(i,i_grid)))) + dcs(i,:) = [y, this % dcs(i,i_grid+1:n_e)] + end do + call move_alloc(dcs, this % dcs) + + electron_energy = [energy_cutoff(PHOTON), electron_energy(i_grid+1:n_e)] + end if + + ! Set incident particle energy grid + if (.not. allocated(ttb_e_grid)) then + call move_alloc(electron_energy, ttb_e_grid) + end if end if ! Take logarithm of energies and cross sections since they are log-log diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index 88c68dd18..1a2e12379 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -536,13 +536,13 @@ contains !p % E = 100.0e6_8 - !if (p % E < energy_cutoff(PHOTON)) return - if (p % E < energy_cutoff(PHOTON)) then - open(unit=13, file="energies.txt", action="write", position="append") - write(13,*) p % E, 0 - close(13) - return - end if + if (p % E < energy_cutoff(PHOTON)) return + !if (p % E < energy_cutoff(PHOTON)) then + ! open(unit=13, file="energies.txt", action="write", position="append") + ! write(13,*) p % E, 0 + ! close(13) + ! return + !end if ! Get bremsstrahlung data for this material mat => ttb(p % material) @@ -573,13 +573,13 @@ contains n = int(y + prn()) E_lost = ZERO - !if (n == 0) return - if (n == 0) then - open(unit=13, file="energies.txt", action="write", position="append") - write(13,*) p % E, n - close(13) - return - end if + if (n == 0) return + !if (n == 0) then + ! open(unit=13, file="energies.txt", action="write", position="append") + ! write(13,*) p % E, n + ! close(13) + ! return + !end if ! Sample index of the tabulated PDF in the energy grid, j or j+1 if (prn() <= f .or. j == 1) then @@ -628,9 +628,9 @@ contains end if end do - open(unit=13, file="energies.txt", action="write", position="append") - write(13,*) p % E, n, photon_energies(:n) - close(13) + !open(unit=13, file="energies.txt", action="write", position="append") + !write(13,*) p % E, n, photon_energies(:n) + !close(13) end subroutine thick_target_bremsstrahlung