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

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

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

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

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

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)