diff --git a/colradpy/ionization_balance_class.py b/colradpy/ionization_balance_class.py index d5afa8b..55f3d17 100644 --- a/colradpy/ionization_balance_class.py +++ b/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 diff --git a/colradpy/saha_boltzmann_lte_pops.py b/colradpy/saha_boltzmann_lte_pops.py index a51f97d..b3f2473 100644 --- a/colradpy/saha_boltzmann_lte_pops.py +++ b/colradpy/saha_boltzmann_lte_pops.py @@ -26,7 +26,7 @@ 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 @@ -34,11 +34,11 @@ def boltzmann_population(w,energy,temp,norm_tot=False): 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)