{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 8.314472\n",
    "T_K = 300\n",
    "dH=80e3\n",
    "dS=10\n",
    "rsys1b = ReactionSystem.from_string(\"\"\"\n",
    "NO + Br -> NOBr; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)\n",
    "\"\"\".format(dH=dH, dS=dS))\n",
    "c0 = 1 # mol/dm3 === 1000 mol/m3\n",
    "kbref = 20836643994.118652*T_K*np.exp(-(dH - T_K*dS)/(R*T_K))/c0\n",
    "kbref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NO0_M = 1.5\n",
    "Br0_M = 0.7\n",
    "init_cond = dict(\n",
    "    NOBr=0*u.M,\n",
    "    NO=NO0_M*u.M,\n",
    "    Br=Br0_M*u.M\n",
    ")\n",
    "t = 5*u.second\n",
    "params = dict(\n",
    "    temperature=T_K*u.K\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytic_b(t):\n",
    "    U, V = NO0_M, Br0_M\n",
    "    d = U - V\n",
    "    return (U*(1 - np.exp(-kbref*t*d)))/(U/V - np.exp(-kbref*t*d))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_and_plot(rsys):\n",
    "    odesys, extra = get_odesys(rsys, unit_registry=si, constants=const)\n",
    "    fig, axes = plt.subplots(1, 4, figsize=(14, 6))\n",
    "    res = odesys.integrate(t, init_cond, params, integrator='cvode')\n",
    "    t_sec = to_unitless(res.xout, u.second)\n",
    "    NOBr_ref = analytic_b(t_sec)\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')] = Br0_M - NOBr_ref\n",
    "    ref[:, odesys.names.index('NO')] = NO0_M - NOBr_ref\n",
    "    axes[0].plot(t_sec, cmp)\n",
    "    axes[1].plot(t_sec, ref)\n",
    "    axes[2].plot(t_sec, cmp - ref)\n",
    "    res.plot_invariant_violations(ax=axes[3])\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(rsys1b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rsys2b = ReactionSystem.from_string(\"\"\"\n",
    "NO + Br -> NOBr; 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(rsys2b)"
   ]
  },
  {
   "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
}