Martin-Synge.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Martin-Synge algorithm\n",
    "\n",
    "\n",
    "Marek Szymanski"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Version 1.1:\n",
    "- units are min/cm/M\n",
    "- parameters from email\n",
    "- tinj parameter : injection time\n",
    "- saving to file: result-Martin-Synge.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate\n",
    "import numpy as np\n",
    "from sparsegrad.base import *\n",
    "from sparsegrad.forward import nvalue, seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "params_type = namedtuple('params_type', ['cinj', 'tinj', 'F', 'a', 'b', 'L', 'u'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cinj_prescribed(t, params):\n",
    "    return where(t<=params.tinj, params.cinj, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_langmuir(c, params):\n",
    "    return params.a*c/(1.+params.b*c)\n",
    "def dfdc_langmuir(c, params):\n",
    "    return params.a/(params.b*c + 1)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ode(t, x, params):\n",
    "    c = stack(cinj_prescribed(t, params), x)\n",
    "    N = len(nvalue(x)) + 1\n",
    "    dz = params.L / (N-1)\n",
    "    deriv_z = (c[1:]-c[:-1])/dz\n",
    "    # dq/dt = dq/dc*dc/dt assuming const params\n",
    "    return -params.u*deriv_z / (1.+params.F*dfdc_langmuir(x, params))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_table = params_type(L=15, \n",
    "                           u=20.62, \n",
    "                           F=0.42,\n",
    "                           cinj=0.1,\n",
    "                           a=15, \n",
    "                           b=50,\n",
    "                           tinj=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "x0 = np.zeros(N)\n",
    "def jac(t, x):\n",
    "    return ode(t, seed(x), params_table).dvalue.tocsr()\n",
    "def func(t, x):\n",
    "    return ode(t, x,params_table)\n",
    "result = scipy.integrate.solve_ivp(func, [0, 10], x0, method='Radau', jac=jac, rtol=1e-6, atol=1e-8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1512059d30>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xt81PWd7/HXJ/eES4AkCCYgIBENVkUj1Xqp1lqxF9CtnuLZtm7Xx3FPj+62x/bRg2cfx+56jo+tuz29bKue46l27bYrUtpu0+qKtpS1WuXmHTAYQSACEggQIJlMZvI5f8wvdhxyGchMfpPJ+/l4+GDm9/v+5vf5tZm88/1+fxdzd0RERArCLkBERHKDAkFERAAFgoiIBBQIIiICKBBERCSgQBAREUCBICIiAQWCiIgACgQREQkUhV3AiaiurvZZs2aFXYaIyKixcePG/e5ek07bURUIs2bNYsOGDWGXISIyapjZjnTbashIREQABYKIiAQUCCIiAigQREQkoEAQERFAgSAiIgEFgoiIAAoEGUEvbDvAM1vbwi5DRAagQJAR8Ye39vO5h9by+YfX8c8vpH2djIiMIAWCjIgH1rzF1AllXDynim891UxXNB52SSKSQoEgWdcdi7P+7XaubjiFOz52Bgc7e/jpxl1hlyUiKRQIknUv7zxEpKeXD51eReNpkzm3rpJH1ykQRHKNAkGybv3b7QB8cE4VZsanzj2VLXs62HHgWMiViUgyBYJk3fb9nUybWEZleTEA18yfBsCTr+8NsywRSaFAkKzbdbCTmVMq3ns/Y0oFH6itZNUmBYJILlEgSNbtau+kbkr5+5ZdOa+Gl3cd4nBXT0hViUiqtALBzBaZWbOZtZjZsn7Wl5rZY8H6tWY2K2ndncHyZjO7Jlg2z8xeTvqvw8y+nKmDktzRHYuztyPCjMkV71t+2Rk19Do8/9b+kCoTkVRDBoKZFQL3AdcCDcBNZtaQ0uwW4KC7zwW+DdwbbNsALAXmA4uA+82s0N2b3f08dz8PuADoBH6RoWOSHPLOwS7ced+QEcB5MyYxvrSIZ95UIIjkinR6CAuBFnff5u5RYDmwJKXNEuCR4PVK4Cozs2D5cnfvdvftQEvwecmuAt5yd12+mod2HewCEvMGyYoLC7j49Cqe2dqGu4dRmoikSCcQaoHkk8Zbg2X9tnH3GHAYqEpz26XAo+mXLKPJ7kOJQKidXH7cusvqq2k92MXO9s6RLktE+pFOIFg/y1L/pBuozaDbmlkJsBj46YA7N7vVzDaY2Ya2Nt0YbbQ5cLQbgOrxJcet+9DpVQCs3dY+ojWJSP/SCYRWYEbS+zpg90BtzKwIqATa09j2WuBFd393oJ27+4Pu3ujujTU1NWmUK7lk/9EoE0qLKC0qPG7d6TXjqR5fwgvbDoRQmYikSicQ1gP1ZjY7+It+KdCU0qYJuDl4fQOw2hMDw03A0uAspNlAPbAuabub0HBRXms/FqWqn94BgJnxwTlVvLDtgOYRRHLAkIEQzAncDqwCtgAr3H2Tmd1tZouDZg8BVWbWAtwBLAu23QSsADYDTwK3uXscwMwqgKuBn2f2kCSXHDjWzZRx/QcCwEVzqth9OMKu9q4RrEpE+lOUTiN3fwJ4ImXZXUmvI8CNA2x7D3BPP8s7SUw8Sx47cDR63BlGyS6eMwVIPDxnZtXA7UQk+3SlsmTVgWPRfieU+2geQSR3KBAka3p7nfZj0UGHjDSPIJI7FAiSNUejMeK9zuSKgQMBNI8gkisUCJI1HcGN6yaWFQ/aLnkeQUTCo0CQrOnoigEwsXzwcxc0jyCSGxQIkjUdkfR6CJpHEMkNCgTJmveGjMoHDwTQPIJILlAgSNZ0RBJDRpVpBILmEUTCp0CQrEl3Uhk0jyCSCxQIkjV9cwjjy4a+IF7zCCLhUyBI1hzu6mFCaRGFBf3dBf14ffMIej6CSDgUCJI1HV2xtCaU+2geQSRcCgTJmqPdPYwvTev+iUDyPIIemCMSBgWCZE1nNE5F6fEPxhmI5hFEwqVAkKzpjMYZV5J+DwES8wh7NI8gEgoFgmTNse4YFSXp9xBA8wgiYVIgSNZ0RuOMO4E5BNA8gkiYFAiSNZ3RGOUn2EPQPIJIeNIKBDNbZGbNZtZiZsv6WV9qZo8F69ea2aykdXcGy5vN7Jqk5ZPMbKWZvWFmW8zs4kwckOSOY91xxp1gIIDmEUTCMmQgmFkhcB9wLdAA3GRmDSnNbgEOuvtc4NvAvcG2DcBSYD6wCLg/+DyA7wJPuvuZwLnAluEfjuSK3l6nqydOxQlOKoPmEUTCkk4PYSHQ4u7b3D0KLAeWpLRZAjwSvF4JXGVmFixf7u7d7r4daAEWmtlE4HLgIQB3j7r7oeEfjuSKrp44AONO4LTTPppHEAlHOoFQC+xKet8aLOu3jbvHgMNA1SDbzgHagB+a2Utm9gMzG3dSRyA56Vg0cafTk+khaB5BJBzpBEJ/N6JJ/ZYO1Gag5UXA+cAD7r4AOAYcNzcBYGa3mtkGM9vQ1taWRrmSCzq7Ez2EEz3ttI/mEURGXjqB0ArMSHpfB+weqI2ZFQGVQPsg27YCre6+Nli+kkRAHMfdH3T3RndvrKmpSaNcyQXD6SGA5hFEwpBOIKwH6s1stpmVkJgkbkpp0wTcHLy+AVjtib5+E7A0OAtpNlAPrHP3vcAuM5sXbHMVsHmYxyI5pCt68nMIoHkEkTAM+eebu8fM7HZgFVAIPOzum8zsbmCDuzeRmBz+ZzNrIdEzWBpsu8nMVpD4ZR8DbnP3ePDRfwn8JAiZbcAXMnxsEqJj0b4ho5PrIaTOIyTOURCRbErr2+ruTwBPpCy7K+l1BLhxgG3vAe7pZ/nLQOOJFCujR2d335DRyfUQIDGP8Pire9hxoJNZ1TrnQCTbdKWyZEVfD+FEb26X7OI5VQA899b+jNQkIoNTIEhWdPZNKp/kHALA6TXjmF5ZxnMtCgSRkaBAkKzozEAPwcy4ZG41z7UcIN6r6xFEsk2BIFnR2R3DDMqKh/cjdll9NYe7enj9ncMZqkxEBqJAkKw4Fo1TUVw47LODLplbDcCzGjYSyToFgmRFZzRGxQk+C6E/1eNLOWv6RJ59U4Egkm0KBMmKk731dX8unVvFxh0H37vYTUSyQ4EgWdEZPblbX/fn0voaovFe1r2tq5ZFskmBIFnRGT3x5ykPZOGsKZQUFvDsm7q5oUg2KRAkK45F4xmZQwAoLynkgtMm83vNI4hklQJBsqKzO0ZFcWZ6CACX1lfzxt4jtB3pzthnisj7KRAkKxKPz8xcIFxWnzj99A+6jYVI1igQJCsiPb2UZrCHMP/USiZVFPPMVgWCSLYoECQrIj1xyjMYCIUFxmX1Nfz71jZ6dRsLkaxQIEhWRHriw75tRaorzqhh/9FuNu/pyOjnikiCAkEyrifeS6zXM9pDAPjwvMQjVH/3xr6Mfq6IJCgQJOMiPYkrissyHAjV40s5p66SNVt1PYJINigQJOMiPb3A8O902p8r5k3lpZ0HOdQZzfhni4x1aX1jzWyRmTWbWYuZLetnfamZPRasX2tms5LW3Rksbzaza5KWv21mr5nZy2a2IRMHI7khWz0EgCvm1dDr8IwuUhPJuCEDwcwKgfuAa4EG4CYza0hpdgtw0N3nAt8G7g22bQCWAvOBRcD9wef1udLdz3N3PVs5j2QzEM6tm8TkimLWaB5BJOPS6SEsBFrcfZu7R4HlwJKUNkuAR4LXK4GrLHEj/CXAcnfvdvftQEvweZLHuoJAyPSkMiROP738DJ1+KpIN6QRCLbAr6X1rsKzfNu4eAw4DVUNs68BTZrbRzG4daOdmdquZbTCzDW1tmkwcDf44h5D5QAC4ct5UDhyL8pqeoiaSUekEQn+PvEr902ygNoNte4m7n09iKOo2M7u8v527+4Pu3ujujTU1NWmUK2HrGzIqL8nOOQuXn1GDGaxp1h8IIpmUzje2FZiR9L4O2D1QGzMrAiqB9sG2dfe+f/cBv0BDSXmjb8iotCg7PYQp40o4t24Sv2vWPIJIJqUTCOuBejObbWYlJCaJm1LaNAE3B69vAFa7uwfLlwZnIc0G6oF1ZjbOzCYAmNk44GPA68M/HMkFf+whZCcQIHG20Suth2g/ptNPRTJlyEAI5gRuB1YBW4AV7r7JzO42s8VBs4eAKjNrAe4AlgXbbgJWAJuBJ4Hb3D0OnAI8a2avAOuAx939ycwemoQlm2cZ9bly3lTc4RldpCaSMWk9wcTdnwCeSFl2V9LrCHDjANveA9yTsmwbcO6JFiujQ9+kcjbOMurzgdpKqsaVsKZ5H9ctSD3HQUROhq5Uloz7Yw8hez9eBQXGh8+oYc3WNmLx3qztR2QsUSBIxvVNKpdlaVK5z9UNp3Cos4f1bx/M6n5ExgoFgmRcpKeXkqICCgr6O+s4cy4/o4aSogKe3vxuVvcjMlYoECTjIj1xyoqy/6M1rrSIS+dW8/SWvSROahOR4VAgSMZFeuJZPeU02dUNp7CrvYs39h4Zkf2J5DMFgmRcV088q6ecJrvqrKmYoWEjkQxQIEjGZfp5yoOZOqGMBTMmKRBEMkCBIBkX6emldIQCAeDqhmm89s5hdh/qGrF9iuQjBYJkXFdPnPIsXoOQ6uqGUwD4zRb1EkSGQ4EgGdc9gnMIAHOnjmdOzTgNG4kMkwJBMq6rJ571i9JSXd1wCs+/dYDDXT0jul+RfKJAkIyL9PSO2GmnfT7WMI1Yr7NGt8QWOWkKBMm4SE88q/cx6s+CGZOoHl+qYSORYVAgSMZ19cSz9nCcgRQUGB89ayprmtvojsVHdN8i+UKBIBnXHcKQEcA186dxtDvGcy37R3zfIvlAgSAZFe91ovHeEZ9UBrhkbjUTy4r49at7RnzfIvlAgSAZNRLPQhhISVEB18yfxtOb39WwkchJSOtba2aLzKzZzFrMbFk/60vN7LFg/Vozm5W07s5gebOZXZOyXaGZvWRmvx7ugUhuGInnKQ/m4+dM50gkxrNvathI5EQNGQhmVgjcB1wLNAA3mVlDSrNbgIPuPhf4NnBvsG0DsBSYDywC7g8+r8+XSDynWfJEJJZ4elkYQ0YAl5xeTWV5MY9r2EjkhKXTQ1gItLj7NnePAsuBJSltlgCPBK9XAleZmQXLl7t7t7tvB1qCz8PM6oBPAD8Y/mFIrujrIZSGMGQEiWGjjzWcomEjkZOQzre2FtiV9L41WNZvG3ePAYeBqiG2/Q7wNUAPxM0jXdG+OYRweggAnzhnOke6Y/x+q4aNRE5EOoHQ33MQUx9PNVCbfpeb2SeBfe6+ccidm91qZhvMbENbW9vQ1Uqo+v4qH6nbX/fnkrnBsNFrGjYSORHpBEIrMCPpfR2we6A2ZlYEVALtg2x7CbDYzN4mMQT1ETP7cX87d/cH3b3R3RtramrSKFfCFOkJ5hBCDITiwgKumZ8YNuobwhKRoaUTCOuBejObbWYlJCaJm1LaNAE3B69vAFZ74iG3TcDS4Cyk2UA9sM7d73T3OnefFXzeanf/bAaOR0IW5mmnyT5xzqkc7Y7xe51tJJK2Ib+1wZzA7cAqEmcErXD3TWZ2t5ktDpo9BFSZWQtwB7As2HYTsALYDDwJ3Obu+pMtj3X1hD+HAPCh06uYVFHM46+mdmZFZCBF6TRy9yeAJ1KW3ZX0OgLcOMC29wD3DPLZa4A16dQhua9vyCjMOQQIho0apvH4a3uCm+2FW4/IaKArlSWjwj7tNNknzpnO0e4Ya5p1MoJIOsL/1kpeieTIkBEkho2qx5fwy5ffCbsUkVFBgSAZ9V4ghHSlcrKiwgI+de6p/HbLPg536klqIkNRIEhGRXp6KSwwigv7uwRl5F2/oJZovJcnXtc1CSJDUSBIRkV64pQVFZC4c0n4PlBbyek14/jFSxo2EhmKAkEyKhLLrTN6zIzrF9Sybns7u9o7wy5HJKcpECSjuqK9ORUIAEvOS9w+q+kVXZMgMhgFgmRUooeQWz9WM6ZUsHDWFH7+YiuJC+hFpD+59c2VUa87Ry8Cu/78Wt5qO8br73SEXYpIzlIgSEZFenJvyAjg42dPp6SwQJPLIoNQIEhGdfXk3pARQGVFMR85cypNr+wmFtcjOET6k3vfXBnVIj3x0O9jNJDrz69l/9Funm3RHVBF+qNAkIyK9MQpzdFAuGJeDZMqilm5sTXsUkRykgJBMirS05sTt63oT2lRIdedV8tTm97l4LFo2OWI5BwFgmRUJEfnEPp85sIZROO9/KtueCdynNz95sqolMtzCABnTZ/IuXWVLF+3S9ckiKRQIEhGRWK5edppss9cOJPmd4/wSuvhsEsRySkKBMmYnngv8V7P6SEjgE+dO53y4kIeW78z7FJEckpa31wzW2RmzWbWYmbL+llfamaPBevXmtmspHV3BsubzeyaYFmZma0zs1fMbJOZ/W2mDkjCkyvPUx7KhLJiPnHOdJpe3s2x7ljY5YjkjCEDwcwKgfuAa4EG4CYza0hpdgtw0N3nAt8G7g22bQCWAvOBRcD9wed1Ax9x93OB84BFZnZRZg5JwpJLT0sbytILZ3AsGufx1/ScBJE+6fQQFgIt7r7N3aPAcmBJSpslwCPB65XAVZa4If4SYLm7d7v7dqAFWOgJR4P2xcF/muEb5bp7ElcAj4ZAuOC0yZxeM47H1u8KuxSRnJFOINQCyd+a1mBZv23cPQYcBqoG29bMCs3sZWAf8LS7r+1v52Z2q5ltMLMNbW16WHou+2MPIbfnECDxnITPXDiDjTsO0rLvSNjliOSEdL65/T36KvWv+YHaDLitu8fd/TygDlhoZmf3t3N3f9DdG929saamJo1yJSxdOfQ85XT8yfl1FBWYegkigXQCoRWYkfS+Dkh90sh7bcysCKgE2tPZ1t0PAWtIzDHIKBYJhozKS0ZHIFSPL+XqhlNYubH1vd6NyFiWTiCsB+rNbLaZlZCYJG5KadME3By8vgFY7YmrfpqApcFZSLOBemCdmdWY2SQAMysHPgq8MfzDkTCNpiGjPp+7+DQOdvbwKz1NTWToQAjmBG4HVgFbgBXuvsnM7jazxUGzh4AqM2sB7gCWBdtuAlYAm4EngdvcPQ5MB35nZq+SCJyn3f3XmT00GWl9gVA6SoaMAC6eU0X91PH86PkdunJZxryidBq5+xPAEynL7kp6HQFuHGDbe4B7Upa9Ciw40WIlt42W6xCSmRmfv/g0/scvN/HyrkMsmDk57JJEQjN6+vaS87pH2RxCn+vPr2N8aRE/en5H2KWIhEqBIBkTifWdZTS6fqzGlxZxwwV1PP7qHtqOdIddjkhoRtc3V3LaaLpSOdVnLzqNaLxX9zeSMU2BIBnTFR09Vyqnmjt1PJfVV/PjF3bqmcsyZikQJGMisTjFhUZhQX/XI+a+z188i70dEZ7e/G7YpYiEQoEgGdMVjVNRktaJaznpI2dOpXZSOY88/3bYpYiEQoEgGdMZjVExys4wSlZYYHzu4tN4YVs7m3br4Tky9igQJGOOReOj7pTTVDddOJOKkkJ+8PvtYZciMuIUCJIxiSGj0R0IlRXFLL1wJr96ZTe7D3WFXY7IiFIgSMYkhoxG7xxCnz+/dBYOPPysegkytigQJGM686CHAFA3uYJPnjOdR9ft5HBXT9jliIwYBYJkTL4EAsCtl8/hWDTOv6zVhWoydigQJGM6u/NjyAhg/qmVXDq3mh8+t53umJ6VIGODAkEyprMnf3oIkOgl7DvSzS9f1rMSZGxQIEjGdObBaafJLquv5qzpE/l/z2zTsxJkTFAgSEbE4r1EY72My5MhI0g8K+HWy2fz5r6j/HbLvrDLEck6BYJkRGdwp9N8GjIC+OQ5pzJjSjn/uPpN9RIk76UVCGa2yMyazazFzJb1s77UzB4L1q81s1lJ6+4Mljeb2TXBshlm9jsz22Jmm8zsS5k6IAlHVzQRCPk0ZARQXFjA7VfO5dXWw6xpbgu7HJGsGjIQzKwQuA+4FmgAbjKzhpRmtwAH3X0u8G3g3mDbBmApMB9YBNwffF4M+Iq7nwVcBNzWz2fKKHKsOwaQV0NGff7k/DrqJpfznd+qlyD5LZ0ewkKgxd23uXsUWA4sSWmzBHgkeL0SuMrMLFi+3N273X070AIsdPc97v4igLsfAbYAtcM/HAlLZ572ECDRS7jtyrm8susQ/75VvQTJX+kEQi2wK+l9K8f/8n6vjbvHgMNAVTrbBsNLC4C16ZctuaYvEPJtDqHPp8+vo3ZSOd9VL0HyWDqB0N/TTlK/EQO1GXRbMxsP/Az4srt39Ltzs1vNbIOZbWhr019nuaozmhgyypcL01KVFCV6CS/tPMQzb+4PuxyRrEgnEFqBGUnv64DUK3Xea2NmRUAl0D7YtmZWTCIMfuLuPx9o5+7+oLs3untjTU1NGuVKGLryvIcAcMMFQS/hN1vVS5C8lE4grAfqzWy2mZWQmCRuSmnTBNwcvL4BWO2Jb0wTsDQ4C2k2UA+sC+YXHgK2uPu3MnEgEq5jYyAQSooK+OIVp/PizkP8Xr0EyUNDBkIwJ3A7sIrE5O8Kd99kZneb2eKg2UNAlZm1AHcAy4JtNwErgM3Ak8Bt7h4HLgE+B3zEzF4O/vt4ho9NRlBXng8Z9bmxMdFL+OZTzfT2qpcg+SWtb6+7PwE8kbLsrqTXEeDGAba9B7gnZdmz9D+/IKPUWOghAJQWFXLH1WfwlZ++wq9f28Pic08NuySRjNGVypIR7512WpzfgQBw3YJazpo+kX9Y9YbuhCp5RYEgGdEVjVFWXEBBQf53/AoLjDuvPZNd7V38+AU9L0HyhwJBMuJod4wJZcVhlzFiLj+jhsvqq/ne6jf1VDXJGwoEyYiOrhgTyvJ7QjnVf1t0Joc6e3hgzVthlyKSEQoEyYiOSM+Y6iEAnF1byfULann4ue28c6gr7HJEhk2BIBnREYkxcYz1EAC+8rEzwOFbT20NuxSRYVMgSEYcifQwcYz1EADqJlfwZ5fM4ucvtfL6O4fDLkdkWBQIkhFHIjEmlo+9HgLAbVfMZUpFCXf98nVdrCajmgJBMqKja+zNIfSprCjmzo+fxYs7D/HTjbuG3kAkRykQZNiisV66Y71MKB2bPQSAT59fy8JZU/i7f3uD9mPRsMsROSkKBBm2I5HEefgTy8dmDwHAzPif153NkUiMe//tjbDLETkpCgQZto5I4sZ2Y+06hFTzpk3glktn89iGXWzc0R52OSInTIEgw/ZeD2GMziEk+9JV9UybWMZf/+J1YvHesMsROSEKBBm2ji71EPqMKy3i659q4I29R3jk+R1hlyNyQhQIMmyaQ3i/RWdP48Nn1PCtp5rZ1d4ZdjkiaVMgyLAd0RzC+5gZ/+u6szEzvvrTV3RtgowaCgQZto6ghzBWr0Poz4wpFdz1qQbWbm/n4ee2h12OSFoUCDJsHZEYZozp6xD6c+MFdXz0rFP4+1XNbH33SNjliAwprUAws0Vm1mxmLWa2rJ/1pWb2WLB+rZnNSlp3Z7C82cyuSVr+sJntM7PXM3EgEp4jkR7GlxSNiYfjnAgz4xuf/gATSov48vKXicZ01pHktiEDwcwKgfuAa4EG4CYza0hpdgtw0N3nAt8G7g22bQCWAvOBRcD9wecB/FOwTEa5Q509VFZouKg/1eNL+bs/+QCb93Twj799M+xyRAaVTg9hIdDi7tvcPQosB5aktFkCPBK8XglcZWYWLF/u7t3uvh1oCT4Pd38G0NU7eeDAsShV40rCLiNnfWz+NG64oI7717SwccfBsMsRGVA6gVALJN+xqzVY1m8bd48Bh4GqNLcdlJndamYbzGxDW1vbiWwqI+TgsSiTFQiD+vqnGpheWc6Xlr/EQd3rSHJUOoHQ38Bw6nl0A7VJZ9tBufuD7t7o7o01NTUnsqmMkPZjUaYoEAY1oayY7//HBezr6Oavlr9EXKeiSg5KJxBagRlJ7+uA3QO1MbMioJLEcFA628ood+BYN1MqFAhDWTBzMncvmc/v39zPN59qDrsckeOkEwjrgXozm21mJSQmiZtS2jQBNwevbwBWu7sHy5cGZyHNBuqBdZkpXXJBVzROpKeXKeMVCOlYunAmNy2cyQNr3uLfXtsTdjki7zNkIARzArcDq4AtwAp332Rmd5vZ4qDZQ0CVmbUAdwDLgm03ASuAzcCTwG3uHgcws0eB54F5ZtZqZrdk9tBkJOw/2g2gSeUT8DeLG1gwcxJf/ekrvKnrEySHWOIP+dGhsbHRN2zYEHYZkmTjjoN8+oE/8MMvXMiV86aGXc6osfdwhE9+71kmlBXxy9sv0Z1iJWvMbKO7N6bTVlcqy7Ds64gAMHVCaciVjC7TKsu4/0/PZ1d7J3/xo41EeuJhlySiQJDh2XckMWR0ysSykCsZfRbOnsI3bzyXF7Yf4PZ/eYkePT9BQqZAkGHZdyRCUYHpLKOTdN2CWu5ePJ/fbHmXr618VXdGlVDpbmQyLHsPd1MzoVT3MRqGz108i45IjH9Y1cyEsiL+dvF8Ehf6i4wsBYIMS+vBTuoml4ddxqj3X644ncNdPTz4zDYqy4v5ysfmhV2SjEEKBBmW1oNdLJw9JewyRj0z485rz6Sjq4fvrW6hO9bLskVnquclI0qBICctFu9lb0eE2knqIWSCmXHP9R+gpKiAB5/Zxp7DEb554zmUFhUOvbFIBigQ5KTtPhQh3usaMsqgwgLjbxfPZ3plOfc++QZtRyL83881UqnnVcsI0FlGctLeajsKwOlTx4dcSX4xM754xel85zPnsXHHQf7D/3me3Ye6wi5LxgAFgpy0ln2JQJhbo0DIhusW1PJPX1jI7kNdLLnvOZ7Zqtu/S3YpEOSkvbnvCFXjSvQshCy6ZG41K7/4ISaVF/P5h9dx968266pmyRoFgpy019/poOHUiWGXkffmTZvAr/7yUv7sQ7N4+LntXHffc7yxtyPssiQPKRDkpER64mx99wjn1FWGXcqYUFZcyN8sns8Pv3Ah+49GWfz95/j+6jfpjMbCLk3yiAJBTspLOw8R63XOmzE57FLGlCvnTWXVly/jqjOn8s2ntnL536/hkT+8TTSm+yDJ8CkQ5KT84a39FBj8urhTAAAHvklEQVR8cI4uShtpVeNLeeCzF/CzL17M6TXj+HrTJj7yv9ewcmOrHs0pw6JAkJOyatNeGk+bovv4h+iC06aw/NaL+NGfL2RyRQlf/ekrXPOdZ/jJ2h0c6oyGXZ6MQgoEOWGvth5i67tH+dS508MuZcwzMy4/o4am2y/hgT89n0Iz/voXr3PhPb/hP/1oA4+/ukdnJUna0rpS2cwWAd8FCoEfuPs3UtaXAj8CLgAOAJ9x97eDdXcCtwBx4K/cfVU6nym563urWxhfWsR1C2rDLkUCZsa1H5jOorOnsWl3B//60js0vbKbpze/y4TSIq6efwoXzanigtMmM6d6nO6mKv0aMhDMrBC4D7gaaAXWm1mTu29OanYLcNDd55rZUuBe4DNm1gAsBeYDpwK/MbMzgm2G+kzJQT9+YQdPb36Xry2axwQNF+UcM+Ps2krOrq3kzo+fxfNvHeBfX36Hpze/y89ffAeASRXFnD9zMhecNpmGUycyc0oFdZPLdc8kSauHsBBocfdtAGa2HFgCJP/yXgL8TfB6JfB9S/wJsgRY7u7dwHYzawk+jzQ+U3LEkUgPzXuP8Oi6XfzsxVY+fEYNt142J+yyZAiFBcal9dVcWl9Nb6+zbf9RNu44yIs7DrFx50FWv7HvvbZmMG1iGTMmVzBjSgXVE0qYVF7CpIpiJpUXU1lRzMSyYspLCiktKqC0qJCy4sS/xYWmHkeeSCcQaoFdSe9bgQ8O1MbdY2Z2GKgKlr+Qsm3fOMNQn5kxn/ze74n0JE7Lc3//WRjHnZPhA78dals/blsfeN0QJ4MMtq/B9tP/+sH2O/i23bFejnYnznUvLSrgLy6fw1evmUdRoaafRpOCAmPu1AnMnTqBz1w4E4BDnVFa9h1lZ3snO9s72dXexa72Tv7w1n4OHIumfSprgUFxYQGFBUaBGWZgwT4LzCiwRM+lwAjeJ9okr+vbJpfkUshNqShhxX++OOv7SScQ+vtfJfXXykBtBlre32+Tfn9FmtmtwK0AM2fOHLjKQcytGU9PPOnjU6pKLTL1B8Het+7kt01dYClrh/7sE9h2kB2fyH6KCwuYXlnGzCnjuHhOFZUVGibKF5MqSmicNYXGWcefOuzuRHp6OdQV5VBnDwc7o3R0xeiOxenu6U38G+sl0pP4NxrrpdedXoded9wTn9H3vve99+9vk/w+p+RYORPKRubG1OnspRWYkfS+Dtg9QJtWMysCKoH2IbYd6jMBcPcHgQcBGhsbT+r/pu8sXXAym4mMWWZGeUkh5SXlTK/U7c3HinT6/euBejObbWYlJCaJm1LaNAE3B69vAFZ7YsyjCVhqZqVmNhuoB9al+ZkiIjKChuwhBHMCtwOrSJwi+rC7bzKzu4EN7t4EPAT8czBp3E7iFzxBuxUkJotjwG3uHgfo7zMzf3giIpIuS528zGWNjY2+YcOGsMsQERk1zGyjuzem01anioiICKBAEBGRgAJBREQABYKIiAQUCCIiAoyys4zMrA3YcZKbVwP7M1jOaKBjzn9j7XhBx3yiTnP3mnQajqpAGA4z25DuqVf5Qsec/8ba8YKOOZs0ZCQiIoACQUREAmMpEB4Mu4AQ6Jjz31g7XtAxZ82YmUMQEZHBjaUegoiIDCLvA8HMFplZs5m1mNmysOvJNjObYWa/M7MtZrbJzL4Udk0jxcwKzewlM/t12LWMBDObZGYrzeyN4P/v7D9SK2Rm9l+Dn+vXzexRMysLu6ZMM7OHzWyfmb2etGyKmT1tZm8G/07Oxr7zOhDMrBC4D7gWaABuMrOGcKvKuhjwFXc/C7gIuG0MHHOfLwFbwi5iBH0XeNLdzwTOJc+P3cxqgb8CGt39bBK3zl8ablVZ8U/AopRly4Dfuns98NvgfcbldSAAC4EWd9/m7lFgObAk5Jqyyt33uPuLwesjJH5J1A6+1ehnZnXAJ4AfhF3LSDCzicDlJJ5FgrtH3f1QuFWNiCKgPHgyYwUDPGlxNHP3Z0g8VybZEuCR4PUjwHXZ2He+B0ItsCvpfStj4JdjHzObBSwA1oZbyYj4DvA1IL0nw49+c4A24IfBMNkPzGxc2EVlk7u/A3wT2AnsAQ67+1PhVjViTnH3PZD4ow+Ymo2d5HsgHPeoeXLu8dnZYWbjgZ8BX3b3jrDrySYz+ySwz903hl3LCCoCzgcecPcFwDGyNIyQK4Jx8yXAbOBUYJyZfTbcqvJLvgdCKzAj6X0dedjFTGVmxSTC4Cfu/vOw6xkBlwCLzextEsOCHzGzH4dbUta1Aq3u3tf7W0kiIPLZR4Ht7t7m7j3Az4EPhVzTSHnXzKYDBP/uy8ZO8j0Q1gP1ZjbbzEpITEA1hVxTVpmZkRhX3uLu3wq7npHg7ne6e527zyLx//Fqd8/rvxzdfS+wy8zmBYuuIvHs8ny2E7jIzCqCn/OryPOJ9CRNwM3B65uBX2ZjJ0XZ+NBc4e4xM7sdWEXijISH3X1TyGVl2yXA54DXzOzlYNl/d/cnQqxJsuMvgZ8Ef+xsA74Qcj1Z5e5rzWwl8CKJs+leIg+vWjazR4ErgGozawW+DnwDWGFmt5AIxhuzsm9dqSwiIpD/Q0YiIpImBYKIiAAKBBERCSgQREQEUCCIiEhAgSAiIoACQUREAgoEEREB4P8DY5+O7nihzj4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1511ff6588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "time, c_out = result.t, result.y[-1,:]\n",
    "plt.plot(time, c_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt('result-Martin-Synge.txt', np.vstack((time, c_out)).transpose(), header='time [min] c_out [M]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5c2a006065a44b6b80bd51b2ec047607",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='i', max=515), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import ipywidgets\n",
    "\n",
    "def plot(i=ipywidgets.IntSlider(min=0,max=len(result.t)-1)):\n",
    "    print('time = ', result.t[i])\n",
    "    plt.plot(result.y[:N-2,i])\n",
    "    plt.show()\n",
    "\n",
    "ipywidgets.interact(plot)"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}