fix the bug of interp_rates_adf11 in ionization_balance_class.py

colradpy/ionization_balance_class.py

@@ -4,18 +4,27 @@
 from scipy.interpolate import RectBivariateSpline
 from colradpy.read_adf11 import *
 
-def interp_rates_adf11(logged_temp,logged_dens,temp,dens,logged_gcr): # y are logged_temp and logged_dens args alongside temp and dens? seems redundant?
+def interp_rates_adf11(adf11_temp,adf11_dens,temp,dens,logged_gcr): # y are logged_temp and logged_dens args alongside temp and dens? seems redundant?
     # Dane did some optimization, Curt seemed to cobble this together (Dane optimized just this block of code as a self-contained unit)
     # parallelization should be implemented if more performance is needed
+
+    '''
+    param adf11_temp: temperature grid from adf11 file
+    param adf11_dens: density grid from adf11 file
+    param temp: temperature grid to interpolate onto
+    param dens: density grid to interpolate onto
+    param logged_gcr: the logged GCR rates from adf11 file
+    return: interpolated GCR rates on the temp and dens grid
+    '''
     gcr_arr = np.zeros( (np.shape(logged_gcr)[0],np.shape(logged_gcr)[1],len(temp),len(dens)) )
 
     logged_temp, logged_dens = np.log10(temp), np.log10(dens) # pre-compute to save time
-    
+
     for i in range(0,np.shape(logged_gcr)[0]):
         for j in range(0,np.shape(logged_gcr)[1]):
-            interp_gcr = RectBivariateSpline(logged_temp,
-                logged_dens,
-                logged_gcr[i,j,:,:],
+            interp_gcr = RectBivariateSpline(adf11_temp,
+                adf11_dens,
+                logged_gcr[i,j,:,:]
             )
             gcr_arr[i,j] = interp_gcr(logged_temp, logged_dens) # array size matching works when Dane tests it
     return 10**gcr_arr

colradpy/saha_boltzmann_lte_pops.py

@@ -26,19 +26,19 @@ def boltzmann_population(w,energy,temp,norm_tot=False):
     """
 
     stat_weight = (w[1:]*2 + 1) / (w[0]*2 + 1)
-    boltz_pop = stat_weight[:,None]*np.exp(-(engergy[1:,None] - energy[0,None])/(temp*8065.313546745508))
+    boltz_pop = stat_weight[:,None]*np.exp(-(energy[1:,None] - energy[0,None])/(temp*8065.313546745508))
 
 
     return boltz_pop
 
 
 def saha_population(ci, w_z, w_z1, energy_z, energy_z1, ftemp, norm_tot=False):
 
-    z = boltzmann_population(w_z,energy_z,temp,norm_tot = norm_tot)
-    z1 = boltzmann_population(w_z1,energy_z1,temp,norm_tot = norm_tot)
+    z = boltzmann_population(w_z,energy_z,ftemp,norm_tot = norm_tot)
+    z1 = boltzmann_population(w_z1,energy_z1,ftemp,norm_tot = norm_tot)
 
 
-    saha = 2*z1/z/ne *( 2*np.pi*temp*8065.313546745508/h**2)**(3/2)
+    saha = 2*z1/z/ne *( 2*np.pi*ftemp*8065.313546745508/h**2)**(3/2)