Allow TTB implementation to handle any photon cutoff energy

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

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

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

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