import numpy as np

def k_CO_OH_to_CO2_H_JPL_2015(T,dens):
    """ Computes rate for CO + OH -> CO2 + H from JPL 2015 """
    rate1 = 1.15e-13*T/300.
    ak0   = 5.9e-33*(T/300.)**(-1.0)
    ak1   = 1.1e-12*(T/300.)**(1.3)
    rate2 = (ak0 * dens) / (1. + ak0 * dens / ak1)
    xpo2  = 1./(1. + np.log10((ak0 * dens) / ak1)**2)
    return rate1 + rate2 * 0.6**xpo2

def k_CO_OH_to_CO2_H_Joshi_2006(T,dens):
    """ Computes rate for CO + OH -> CO2 + H from Joshi and Wang 2006 """
    k1a0 = 1.34*2.5*dens                                  \
         * 1/(1/(3.62e-26*T**(-2.739)*np.exp(-20./T))  \
         + 1/(6.48e-33*T**(0.14)*np.exp(-57./T)))
    k1b0 = 1.17e-19*T**(2.053)*np.exp(139./T)          \
         + 9.56e-12*T**(-0.664)*np.exp(-167./T)
    k1ainf = 1.52e-17*T**(1.858)*np.exp(28.8/T)        \
           + 4.78e-8*T**(-1.851)*np.exp(-318./T)
    x = k1a0/(k1ainf - k1b0)
    y = k1b0/(k1ainf - k1b0)
    fc = 0.628*np.exp(-1223./T) + (1. - 0.628)*np.exp(-39./T)  \
       + np.exp(-T/255.)
    fx = fc**(1./(1. + (np.log(x))**2))
    k1a = k1a0*((1. + y)/(1. + x))*fx
    k1b = k1b0*(1./(1.+x))*fx
    return  k1a + k1b

k_lib = {'co + oh -> co2 + h':{'JPL2015':k_CO_OH_to_CO2_H_JPL_2015,
                               'Joshi2006':k_CO_OH_to_CO2_H_Joshi_2006}}