Fix yiled calc if max energy is small

numaryu · Sep 27, 2024 · 2d507aa · 2d507aa
1 parent 1b21586
commit 2d507aa
Showing 1 changed file with 50 additions and 38 deletions.
diff --git 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