In [1]:

    

%matplotlib inline
import pylab as pl


    

In [2]:

    

import os
if not os.path.exists('oh2co-h2.dat'):
    import urllib
    urllib.urlretrieve('http://home.strw.leidenuniv.nl/~moldata/datafiles/oh2co-h2.dat')


    

In [3]:

    

import pyradex
import numpy as np
tau = np.empty([20,4])


    

In [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






    
Out[4]:





88

In [5]:

    

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]


    

In [6]:

    

pl.loglog(densities,tau,'o')


    

    
Out[6]:





[<matplotlib.lines.Line2D at 0x16459a390>,
 <matplotlib.lines.Line2D at 0x164879450>,
 <matplotlib.lines.Line2D at 0x164879650>,
 <matplotlib.lines.Line2D at 0x1648797d0>]






    

In [7]:

    

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


    

    
Out[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])]

In [8]:

    

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)))


    

    
Out[8]:





[<matplotlib.lines.Line2D at 0x164e7ef50>]






    

In [9]:

    

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])


    

    
Out[9]:





[<matplotlib.lines.Line2D at 0x165037e10>]






    

In [22]:

    

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


    

In [23]:

    

tau1 = buildgrid()


    

In [24]:

    

pl.loglog(np.logspace(0,8,20), tau1)


    

    
Out[24]:





[<matplotlib.lines.Line2D at 0x1656d9e10>]






    

In [ ]: