Astrochem issue #53

input.ini

[files]
source = source.mdl
chem = osu2009.chm
# Physical paramaters
[phys]
chi = 1.0
cosmic = 1e-15
# Solver parameters
[solver]
ti = 1e-6
tf = 1e7
# Initial abundances
[abundances]
H2      = 0.5
He      = 0.14
N       = 2.14e-5
O       = 1.76e-4
C(+)    = 7.30e-5
S(+)    = 8.00e-8
Si(+)   = 8.00e-9
Fe(+)   = 3.00e-9
Na(+)   = 2.00e-9
Mg(+)   = 7.00e-9
P(+)    = 2.00e-10
Cl(+)   = 1.00e-9
F       = 6.68e-9
e(-)    = 7.31012e-5
# Output
[output]
abundances = CO,HCO(+)

show-bug.py

from astrochem import tools
from astrochem.wrapper import *
import numpy as np
import matplotlib.pyplot as plt

# Read the abundances computed by the command line version

time_co, abun_co_cl = tools.readabun("astrochem_output.h5", 'CO')
time_hcop, abun_hcop_cl = tools.readabun("astrochem_output.h5", 'HCO(+)')

# Compute the abundances with the Python version

p = phys()
p.chi = 1.0
p.cosmic = 1e-15

initial_abundances = {
    "H2": 0.5,
    "He": 0.14,
    "N": 2.14e-5,
    "O": 1.76e-4,
    "C(+)": 7.30e-5,
    "S(+)": 8.00e-8,
    "Si(+)": 8.00e-9,
    "Fe(+)": 3.00e-9,
    "Na(+)": 2.00e-9,
    "Mg(+)": 7.00e-9,
    "P(+)": 2.00e-10,
    "Cl(+)": 1.00e-9,
    "F": 6.68e-9,
    "e(-)": 7.31012e-5
}

c = cell(av=20, nh=1e4, tgas=10, tdust=10)
s = solver(c, "osu2009.chm", p, 1e-15, 1e-6, initial_abundances, 1e4, 0.)

abun_co = []
abun_hcop = []
for ti in time_co: # use the same time values than for the CL version
    abundances = s.solve(ti*365*24*3600, c) # time should be given in years !!!!
    abun_co.append(abundances['CO']*.9)
    abun_hcop.append(abundances['HCO(+)'])

# Plot the results

plt.xlabel("Time (years)")
plt.ylabel("Abundances")
plt.title("CO abundance")
plt.yscale('log')
plt.xscale('log')
plt.xlim([1e-6, 1e7])
plt.ylim([1e-15, 1e-4])
plt.plot(time_co, abun_co_cl, label="Command line")
plt.plot(time_co, abun_co, label="Python interface")
plt.legend()
fig = plt.gcf()
plt.show()

source.mdl

# Source model file example
# shell number, Av [mag], nH [cm^-3], Tgas [K], Tdust [K]
#
0	20.0	1e+04	10.0	10.0