From 2d507aab0408006ed15e006d4f3f77602699109d Mon Sep 17 00:00:00 2001
From: Ryusuke Numata <ryusuke.numata@gmail.com>
Date: Fri, 27 Sep 2024 14:45:59 +0900
Subject: [PATCH] Fix yiled calc if max energy is small

---
 src/class_orbital.f90 | 88 ++++++++++++++++++++++++-------------------
 1 file changed, 50 insertions(+), 38 deletions(-)

diff --git a/src/class_orbital.f90 b/src/class_orbital.f90
index 24d5165..8054118 100644
--- a/src/class_orbital.f90
+++ b/src/class_orbital.f90
@@ -429,46 +429,58 @@ subroutine calculate_yield(self)
 
     do io = 1, self%number
 
-       ne1 = egrid%grid_number(self%energy_ionize(io))
-       ne2 = egrid%grid_number(self%energy_singlet(io))
-       ne3 = egrid%grid_number(self%energy_triplet(io))
-
-       ! energy >=  energy_ionize
-       do ie = 1, ne1
-          self%total_cross_section_ionize_orbital(io, ie) = &
-               integrate_sigma(self, "direct", io, &
-               egrid%val(ie), self%energy_ionize(io), egrid%val(ie))
-
-          self%total_cross_section_singlet_orbital(io, ie) = &
-               integrate_sigma(self, "direct", io, &
-               egrid%val(ie), self%energy_singlet(io), self%energy_ionize(io)) &
-               + 0.5*integrate_sigma(self, "exchange", io, &
-               egrid%val(ie), self%energy_singlet(io), self%energy_ionize(io))
-
-          self%total_cross_section_triplet_orbital(io, ie) = &
-               0.5*integrate_sigma(self, "exchange", io, &
-               egrid%val(ie), self%energy_triplet(io), self%energy_ionize(io) )
-       end do
+       if (self%energy_ionize(io) < egrid%val_max) then
+          ne1 = egrid%grid_number(self%energy_ionize(io))
 
-       ! energy_ionize > energy >=  energy_singlet
-       do ie = ne1+1, ne2
-          self%total_cross_section_singlet_orbital(io, ie) = &
-               integrate_sigma(self, "direct", io, &
-               egrid%val(ie), self%energy_singlet(io), egrid%val(ie)) &
-               + 0.5*integrate_sigma(self, "exchange", io, &
-               egrid%val(ie), self%energy_singlet(io), egrid%val(ie))
-
-          self%total_cross_section_triplet_orbital(io, ie) = &
-               0.5*integrate_sigma(self, "exchange", io, &
-               egrid%val(ie), self%energy_triplet(io), egrid%val(ie))
-       end do
+          ! energy >=  energy_ionize
+          do ie = 1, ne1
+             self%total_cross_section_ionize_orbital(io, ie) = &
+                  integrate_sigma(self, "direct", io, &
+                  egrid%val(ie), self%energy_ionize(io), egrid%val(ie))
+
+             self%total_cross_section_singlet_orbital(io, ie) = &
+                  integrate_sigma(self, "direct", io, &
+                  egrid%val(ie), self%energy_singlet(io), self%energy_ionize(io)) &
+                  + 0.5*integrate_sigma(self, "exchange", io, &
+                  egrid%val(ie), self%energy_singlet(io), self%energy_ionize(io))
+
+             self%total_cross_section_triplet_orbital(io, ie) = &
+                  0.5*integrate_sigma(self, "exchange", io, &
+                  egrid%val(ie), self%energy_triplet(io), self%energy_ionize(io) )
+          end do
+       else
+          ne1 = 0
+       end if
 
-       ! singlet_energy > energy >=  energy_triplet
-       do ie = ne2+1, ne3
-          self%total_cross_section_triplet_orbital(io, ie) = &
-               0.5*integrate_sigma(self, "exchange", io, &
-               egrid%val(ie), self%energy_triplet(io), egrid%val(ie))
-       end do
+       if (self%energy_singlet(io) < egrid%val_max) then
+          ne2 = egrid%grid_number(self%energy_singlet(io))
+
+          ! energy_ionize > energy >=  energy_singlet
+          do ie = ne1+1, ne2
+             self%total_cross_section_singlet_orbital(io, ie) = &
+                  integrate_sigma(self, "direct", io, &
+                  egrid%val(ie), self%energy_singlet(io), egrid%val(ie)) &
+                  + 0.5*integrate_sigma(self, "exchange", io, &
+                  egrid%val(ie), self%energy_singlet(io), egrid%val(ie))
+
+             self%total_cross_section_triplet_orbital(io, ie) = &
+                  0.5*integrate_sigma(self, "exchange", io, &
+                  egrid%val(ie), self%energy_triplet(io), egrid%val(ie))
+          end do
+       else
+          ne2 = 0
+       end if
+
+       if (self%energy_triplet(io) < egrid%val_max) then
+          ne3 = egrid%grid_number(self%energy_triplet(io))
+
+          ! singlet_energy > energy >=  energy_triplet
+          do ie = ne2+1, ne3
+             self%total_cross_section_triplet_orbital(io, ie) = &
+                  0.5*integrate_sigma(self, "exchange", io, &
+                  egrid%val(ie), self%energy_triplet(io), egrid%val(ie))
+          end do
+       end if
 
     end do