{
 "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",
    "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": [
    "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 = 5*u.second\n",
    "params = dict(\n",
    "    temperature=T_K*u.K\n",
    ")"
   ]
  },
  {
   "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, 3, 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 = NOBr0_M*np.exp(-kref*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')] = NOBr0_M - NOBr_ref\n",
    "    ref[:, odesys.names.index('NO')] = NOBr0_M - NOBr_ref\n",
    "    assert np.allclose(cmp, ref)\n",
    "    axes[0].plot(t_sec, cmp)\n",
    "    axes[1].plot(t_sec, cmp - ref)\n",
    "    res.plot_invariant_violations(ax=axes[2])"
   ]
  },
  {
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}