An incomplete example showing attempts to fit the RADEX parameter space with low-order polynomials
%matplotlib inline
import pylab as pl
import os
if not os.path.exists('oh2co-h2.dat'):
import urllib
urllib.urlretrieve('http://home.strw.leidenuniv.nl/~moldata/datafiles/oh2co-h2.dat')
import pyradex
import numpy as np
tau = np.empty([20,4])
abundance = 10**-8.5
R = pyradex.Radex(species='oh2co-h2', abundance=abundance, column=1e14, temperature=20)
R.run_radex()
/Users/adam/repos/astroquery/astroquery/lamda/__init__.py:21: UserWarning: Experimental: LAMDA has not yet been refactored to have its API match the rest of astroquery. warnings.warn("Experimental: LAMDA has not yet been refactored to have its API match the rest of astroquery.") WARNING: Assuming the density is n(H_2). [pyradex.core] WARNING:astropy:Assuming the density is n(H_2). /Users/adam/repos/pyradex/pyradex/core.py:483: UserWarning: Using a default ortho-to-para ratio (which will only affect species for which independent ortho & para collision rates are given) warnings.warn("Using a default ortho-to-para ratio (which " WARNING: Assuming the density is n(H_2). [pyradex.core] WARNING:astropy:Assuming the density is n(H_2). /Users/adam/repos/pyradex/pyradex/core.py:800: RuntimeWarning: invalid value encountered in divide frac_level_diff = level_diff/self.level_population
88
fortho = 1e-3
densities = np.logspace(0,8,20)
for jj,dd in enumerate(densities):
R.density = {'oH2':dd*fortho,'pH2':dd*(1-fortho)}
R.abundance = abundance # reset column to the appropriate value
R.run_radex(reuse_last=False, reload_molfile=True)
tau[jj,:] = R.tau[:4]
pl.loglog(densities,tau,'o')
[<matplotlib.lines.Line2D at 0x16459a390>, <matplotlib.lines.Line2D at 0x164879450>, <matplotlib.lines.Line2D at 0x164879650>, <matplotlib.lines.Line2D at 0x1648797d0>]
fits_7 = [np.polyfit(np.log10(densities),np.log10(t),7) for t in tau.T]
fits_11 = [np.polyfit(np.log10(densities),np.log10(t),11) for t in tau.T]
fits_7
[array([ 4.86015567e-04, -1.45822371e-02, 1.69408821e-01, -9.47638215e-01, 2.57322263e+00, -3.05652747e+00, 2.30492326e+00, -3.97239753e+00]), array([ 2.39146864e-04, -7.13724406e-03, 8.25362832e-02, -4.60721266e-01, 1.25940091e+00, -1.55001397e+00, 1.66597043e+00, -2.90593168e+00]), array([ 5.87632284e-05, -2.77338874e-03, 4.41794386e-02, -3.16452339e-01, 1.05884924e+00, -1.50480111e+00, 1.74907731e+00, -4.85548112e+00]), array([ 1.47711558e-04, -4.58305035e-03, 5.49716696e-02, -3.17766208e-01, 9.00971669e-01, -1.16643424e+00, 1.52717654e+00, -2.93472028e+00])]
pl.loglog(densities,tau[:,0],'o')
pl.loglog(densities,10**np.polyval(fits_7[0],np.log10(densities)))
pl.loglog(densities,10**np.polyval(fits_11[0],np.log10(densities)))
[<matplotlib.lines.Line2D at 0x164e7ef50>]
pl.semilogx(densities,abs(tau[:,0]-10**np.polyval(fits_7[0],np.log10(densities)))/tau[:,0])
pl.semilogx(densities,abs(tau[:,0]-10**np.polyval(fits_11[0],np.log10(densities)))/tau[:,0])
[<matplotlib.lines.Line2D at 0x165037e10>]
def buildgrid(densities=np.logspace(0,8,20), abundance=10**-8.5, fortho=1e-3, **kwargs):
tau = np.zeros_like(densities)
R = pyradex.Radex(species='co', abundance=abundance, collider_densities={'oH2':10,'pH2':10},
temperature=20)
for jj,dd in enumerate(densities):
R.density = {'oH2':dd*fortho,'pH2':dd*(1-fortho)}
R.abundance = abundance # reset column to the appropriate value
R.run_radex(**kwargs)
tau[jj] = R.tau[0]
return tau
tau1 = buildgrid()
pl.loglog(np.logspace(0,8,20), tau1)
[<matplotlib.lines.Line2D at 0x1656d9e10>]