improved fang 2010 integration limits

Author: julesvanirsel

src/ionization/fang.f90

@@ -68,17 +68,29 @@ elemental real(wp) function fang2010_spectrum(Q0_keV, E0_keV, Tn, massden_gcm3,
 qtot = 0
 do k=1,lbins
   !! energy bin limits
-  bin_lb = -1.0_wp ! lower limit from Fang 2008, 2010
   select case (diff_num_flux)
   case (0) ! Maxwellian
-    bin_ub = log10(20*E0_keV) ! gives 1.12e-7 x max(phi)
+    bin_lb = max(-1.0_wp, log10(E0_keV / 300.0_wp))
+    bin_ub = min(3.0_wp, log10(8.0_wp * E0_keV))
+  case (1) ! kappa distribution
+    bin_lb = max(-1.0_wp, log10(E0_keV / 300.0_wp))
+    bin_ub = min(3.0_wp, log10(8.0_wp * E0_keV * (kappa + 8.0_wp) / kappa))
+  case (2) ! multiple Maxwellians
+    bin_lb = max(-1.0_wp, log10(E0_keV / 300.0_wp))
+    bin_ub = min(3.0_wp, log10(8.0_wp * E0_keV * bimax_frac))
   case (3) ! accelerated Maxwellian
-    bin_lb = log10(E0_keV) ! minimum energy of spectrum
-    bin_ub = log10(10*max(E0_keV, E0_char_keV))
+    bin_lb = max(-1.0_wp, log10(E0_keV)) ! minimum energy of spectrum
+    if (E0_keV < E0_char_keV) then
+      bin_ub = min(3.0_wp, log10(8.0_wp * E0_char_keV))
+    else
+      bin_ub = min(3.0_wp, log10(6.0_wp * E0_char_keV + E0_keV))
+    endif
   case (4) ! Evans, D. S. (1974) 2 keV acc., 0.8 keV temp
-    bin_ub = log10(10*2.0_wp)
+    bin_lb = -1.0_wp
+    bin_ub = log10(40.0_wp)
   case default
-    bin_ub = 3 ! upper limit from Fang 2008, 2010
+    bin_lb = -1.0_wp ! lower limit from Fang 2008, 2010
+    bin_ub = 3.0_wp ! upper limit from Fang 2008, 2010
   end select
 
   !! log bins, midpoint rule
@@ -105,27 +117,25 @@ elemental real(wp) function fang2010_spectrum(Q0_keV, E0_keV, Tn, massden_gcm3,
   !! user choice of differential energy flux
   select case (diff_num_flux)
   case (0) ! Maxwellian
-    phi_keV = (Q0_keV/(2*E0_keV**3)) * Ebin_keV * exp(-1*Ebin_keV/E0_keV) ! 1/keV/s/cm^2
+    phi_keV = (Q0_keV/(2*E0_keV**2)) * (Ebin_keV/E0_keV) * exp(-1*Ebin_keV/E0_keV) ! 1/keV/s/cm^2
   case (1) ! kappa distribution
     if (kappa <= 2) then
       error stop 'ERROR:ionize_fang:fang2010_spectrum: for finite <E>, kappa must be greater than 2'
     end if
-    phi_keV = (Q0_keV/(2*E0_keV**3)) * (gamma(kappa+1)/(gamma(kappa-1/2)*(kappa**(3/2)))) &
-      * Ebin_keV * (1+Ebin_keV/(kappa*E0_keV))**(-1-1*kappa)
+    phi_keV = (Q0_keV/(2*E0_keV**2)) * (Ebin_keV/E0_keV) * ((kappa-1)*(kappa-2)/kappa**2) * (1+Ebin_keV/E0_keV/kappa)**(-1-1*kappa)
   case (2) ! multiple Maxwellians
-    if (bimax_frac <= 0) then
-      error stop 'ERROR:ionize_fang:fang2010_spectrum: bimax_frac must be greater than 0'
+    if (bimax_frac <= 1) then
+      error stop 'ERROR:ionize_fang:fang2010_spectrum: bimax_frac must be greater than 1'
     end if
-    phi_keV = (Q0_keV/(2*E0_keV**3)) * (Ebin_keV/bimax_frac**2) * exp(-1*Ebin_keV/(bimax_frac*E0_keV)) ! secondary Maxwellian
-    phi_keV = phi_keV + (Q0_keV/(2*E0_keV**3)) * Ebin_keV * exp(-1*Ebin_keV/E0_keV) ! add primary Maxwellian
-    phi_keV = phi_keV/2 ! normalize
+    phi_keV = (Q0_keV/(2*E0_keV**2)) * (Ebin_keV/E0_keV) * (1/(1+bimax_frac**2)) &
+      * (exp(-1*Ebin_keV/E0_keV) + (1/bimax_frac)*exp(-1*Ebin_keV/E0_keV/bimax_frac))
   case (3) ! accelerated Maxwellian
     if (E0_char_keV < 0) then
       error stop 'ERROR:ionize_fang:fang2010_spectrum: W0_char must be greater than 0'
     end if
-    if (Ebin_keV<E0_keV) then ! no electrons with energies less than U0
+    if (Ebin_keV<E0_keV) then ! no electrons with energies less than potential drop
       phi_keV = 0
-    else ! E0_kev = acc. potential, E0_char_keV = thermal energy
+    else ! E0_keV = acc. potential, E0_char_keV = thermal energy
       phi_keV = (Q0_keV/(E0_char_keV**2+(E0_char_keV+E0_keV)**2)) * (Ebin_keV/E0_char_keV) * exp(-1*(Ebin_keV-E0_keV)/E0_char_keV)
     end if
   case (4) ! Data used from Evans, D. S. (1974)