{
 "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": [
    "t, t0, A, B, C1 = sm.symbols('t t0 A B C1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = C1/sm.E**((A*(((t + t0)*(-B + t + t0))/sm.E**(B/(t + t0)) - B**2*sm.Ei(-(B/(t + t0)))))/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y"
   ]
  },
  {
   "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": [
    "y.subs(t, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yunit0 = y.subs(C1, C1/y.subs(t, 0)).simplify()\n",
    "yunit0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import expi\n",
    "f = sm.lambdify([t, t0, A, B], yunit0, modules=['numpy', {'Ei': expi}])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 8.314472\n",
    "T_K = 290\n",
    "kB = 1.3806504e-23\n",
    "h = 6.62606896e-34\n",
    "dH = 80e3\n",
    "dS = 10\n",
    "rsys1 = ReactionSystem.from_string(\"\"\"\n",
    "NOBr -> NO + Br; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)\n",
    "\"\"\".format(dH=dH, dS=dS))\n",
    "kref = 20836643994.118652*T_K*np.exp(-(dH - T_K*dS)/(R*T_K))\n",
    "kref"
   ]
  },
  {
   "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]), 290, _A, _B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NOBr0_M = 0.7\n",
    "init_cond = dict(\n",
    "    NOBr=NOBr0_M*u.M,\n",
    "    NO=0*u.M,\n",
    "    Br=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, 1*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",
    "        NOBr_ref = NOBr0_M*f(t_sec, T_K, _A, _B)\n",
    "        cmp = to_unitless(res.yout, u.M)\n",
    "        ref = np.empty_like(cmp)\n",
    "        ref[:, odesys.names.index('NOBr')] = NOBr_ref\n",
    "        ref[:, odesys.names.index('Br')] = NOBr0_M - NOBr_ref\n",
    "        ref[:, odesys.names.index('NO')] = NOBr0_M - NOBr_ref\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",
    "NOBr -> NO + Br; 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)"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}