{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "from chempy import ReactionSystem\n",
    "from chempy.kinetics.ode import get_odesys\n",
    "sp.init_printing()\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ReactionSystem?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rsys = ReactionSystem.from_string(\"\"\"\n",
    "-> H; 'p'\n",
    "2 H -> H2; 'k2'\n",
    "\"\"\", checks=('substance_keys', 'duplicate', 'duplicate_names'))\n",
    "rsys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H, k, p, c1, t, H0  = sp.symbols('H k p c1 t H0', real=True, positive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dHdt = p - 2*k*H**2\n",
    "dHdt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analytic = sp.sqrt(p/k/2)*sp.tanh((c1 + t)*sp.sqrt(2*k*p))\n",
    "analytic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(analytic.diff(t) - dHdt.subs(H, analytic)).simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol1, sol2 = sp.solve(analytic.subs(t, 0) - H0, c1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res1 = analytic.subs(c1, sol1).simplify()\n",
    "res1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res2 = analytic.subs(c1, sol2).simplify()\n",
    "res2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_python_code(expr, be='np'):\n",
    "    cses, new_expr = sp.cse(expr)\n",
    "    snippet = '\\n'.join(['%s = %s' % cse for cse in cses] + ['return %s' % new_expr[0]]).replace(\n",
    "        'sqrt', '%s.sqrt' % be).replace(\n",
    "        'log', '%s.log' % be).replace(\n",
    "        'tanh', '%s.tanh' % be).replace(\n",
    "        'exp', '%s.exp' % be)\n",
    "    return snippet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(get_python_code(res2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(get_python_code(res2.rewrite('exp'), ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1 = sp.lambdify([t, p, k, H0], res1)\n",
    "f2 = sp.lambdify([t, p, k, H0], res2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "odesys, extra = get_odesys(rsys, include_params=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c0 = {'H': 42e-6, 'H2': 17e-6}\n",
    "params = {'p': 314/3600*.998*2*1.036426856e-7, 'k2': 53/60}\n",
    "result = odesys.integrate(7*3600, c0, params,\n",
    "                          integrator='cvode')\n",
    "result.plot()\n",
    "plt.plot(result.xout, f2(result.xout, params['p'], params['k2'], c0['H']), ls=':')\n",
    "plt.gca().set_ylim([0, 300e-6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result.params"
   ]
  },
  {
   "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
}