{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sm\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from chempy import ReactionSystem\n",
    "from chempy.units import to_unitless, SI_base_registry as si, default_units as u, default_constants as const\n",
    "from chempy.kinetics.ode import get_odesys\n",
    "from chempy.kinetics.rates import RampedTemp\n",
    "sm.init_printing()\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "symbs = t, k, m, A, B, C1 = sm.symbols('t k m A B C1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = -sm.E**(B/(m + k*t))*k/(\n",
    "     A*B*m - A*m**2 + A*B*k*t - 2*A*k*m*t - A*k**2*t**2 +\n",
    "     sm.E**(B/(m + k*t))*k*C1 +\n",
    "     A*B**2*sm.E**(B/(m + k*t))*sm.Ei(-(B/(m + k*t)))\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(y.diff(t)/y).simplify().expand().simplify().factor().powsimp(force=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_C1, = sm.solve(y.subs(t, 0) - 1, C1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yunit0 = y.subs(C1, _C1).simplify()\n",
    "yunit0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sm.python(yunit0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import expi\n",
    "f = sm.lambdify(symbs[:-1], yunit0, modules=['numpy', {'Ei': expi}])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 8.314472\n",
    "T_K = 290\n",
    "dTdt_Ks = 3\n",
    "kB = 1.3806504e-23\n",
    "h = 6.62606896e-34\n",
    "dH = 80e3\n",
    "dS = 10\n",
    "rsys1 = ReactionSystem.from_string(\"\"\"\n",
    "2 NO2 -> N2O4; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)\n",
    "\"\"\".format(dH=dH, dS=dS))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_A = kB/h*np.exp(dS/R)\n",
    "_B = dH/R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f(np.array([0, 1, 5, 20]), dTdt_Ks, T_K, _A, _B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NO2_M = 1.0\n",
    "init_cond = dict(\n",
    "    NO2=NO2_M*u.M,\n",
    "    N2O4=0*u.M\n",
    ")\n",
    "t = 20*u.second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_and_plot(rsys):\n",
    "    odes, extra = get_odesys(rsys, unit_registry=si, constants=const, substitutions={\n",
    "        'temperature': RampedTemp([T_K*u.K, dTdt_Ks*u.K/u.s])})\n",
    "    fig, all_axes = plt.subplots(2, 3, figsize=(14, 6))\n",
    "    for axes, odesys in zip(all_axes, [odes, odes.as_autonomous()]):\n",
    "        res = odesys.integrate(t, init_cond, integrator='cvode')\n",
    "        t_sec = to_unitless(res.xout, u.second)\n",
    "        NO2_ref = f(t_sec, dTdt_Ks, T_K, _A, _B)\n",
    "        cmp = to_unitless(res.yout, u.M)\n",
    "        ref = np.empty_like(cmp)\n",
    "        ref[:, odesys.names.index('NO2')] = NO2_ref\n",
    "        ref[:, odesys.names.index('N2O4')] = (NO2_M - NO2_ref)/2\n",
    "        axes[0].plot(t_sec, cmp)\n",
    "        axes[1].plot(t_sec, cmp - ref)\n",
    "        res.plot_invariant_violations(ax=axes[2])\n",
    "        assert np.allclose(cmp, ref)\n",
    "        print({k: v for k, v in res.info.items() if not k.startswith('internal')})    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate_and_plot(rsys1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rsys2 = ReactionSystem.from_string(\"\"\"\n",
    "2 NO2 -> N2O4; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol]))\n",
    "\"\"\".format(dH=dH, dS=dS))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate_and_plot(rsys2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}