In [ ]:

    

import sympy as sm
import matplotlib.pyplot as plt
import numpy as np
from chempy import ReactionSystem
from chempy.units import to_unitless, SI_base_registry as si, default_units as u, default_constants as const
from chempy.kinetics.ode import get_odesys
from chempy.kinetics.rates import RampedTemp
sm.init_printing()
%matplotlib inline


    

In [ ]:

    

symbs = t, k, m, A, B, C1 = sm.symbols('t k m A B C1')


    

In [ ]:

    

y = -sm.E**(B/(m + k*t))*k/(
     A*B*m - A*m**2 + A*B*k*t - 2*A*k*m*t - A*k**2*t**2 +
     sm.E**(B/(m + k*t))*k*C1 +
     A*B**2*sm.E**(B/(m + k*t))*sm.Ei(-(B/(m + k*t)))
)


    

In [ ]:

    

(y.diff(t)/y).simplify().expand().simplify().factor().powsimp(force=True)


    

In [ ]:

    

_C1, = sm.solve(y.subs(t, 0) - 1, C1)


    

In [ ]:

    

yunit0 = y.subs(C1, _C1).simplify()
yunit0


    

In [ ]:

    

print(sm.python(yunit0))


    

In [ ]:

    

from scipy.special import expi
f = sm.lambdify(symbs[:-1], yunit0, modules=['numpy', {'Ei': expi}])


    

In [ ]:

    

R = 8.314472
T_K = 290
dTdt_Ks = 3
kB = 1.3806504e-23
h = 6.62606896e-34
dH = 80e3
dS = 10
rsys1 = ReactionSystem.from_string("""
2 NO2 -> N2O4; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)
""".format(dH=dH, dS=dS))


    

In [ ]:

    

_A = kB/h*np.exp(dS/R)
_B = dH/R


    

In [ ]:

    

f(np.array([0, 1, 5, 20]), dTdt_Ks, T_K, _A, _B)


    

In [ ]:

    

NO2_M = 1.0
init_cond = dict(
    NO2=NO2_M*u.M,
    N2O4=0*u.M
)
t = 20*u.second


    

In [ ]:

    

def integrate_and_plot(rsys):
    odes, extra = get_odesys(rsys, unit_registry=si, constants=const, substitutions={
        'temperature': RampedTemp([T_K*u.K, dTdt_Ks*u.K/u.s])})
    fig, all_axes = plt.subplots(2, 3, figsize=(14, 6))
    for axes, odesys in zip(all_axes, [odes, odes.as_autonomous()]):
        res = odesys.integrate(t, init_cond, integrator='cvode')
        t_sec = to_unitless(res.xout, u.second)
        NO2_ref = f(t_sec, dTdt_Ks, T_K, _A, _B)
        cmp = to_unitless(res.yout, u.M)
        ref = np.empty_like(cmp)
        ref[:, odesys.names.index('NO2')] = NO2_ref
        ref[:, odesys.names.index('N2O4')] = (NO2_M - NO2_ref)/2
        axes[0].plot(t_sec, cmp)
        axes[1].plot(t_sec, cmp - ref)
        res.plot_invariant_violations(ax=axes[2])
        assert np.allclose(cmp, ref)
        print({k: v for k, v in res.info.items() if not k.startswith('internal')})


    

In [ ]:

    

integrate_and_plot(rsys1)


    

In [ ]:

    

rsys2 = ReactionSystem.from_string("""
2 NO2 -> N2O4; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol]))
""".format(dH=dH, dS=dS))


    

In [ ]:

    

integrate_and_plot(rsys2)


    

In [ ]: