dive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dive\n",
    "\n",
    "wrapped in a python class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    " - Dive profile is usually shown as a series of depth and time, in `MM:SS` format, points.\n",
    "   - Need to convert the latter into decimal minutes,\n",
    "   - and convert it to time and current depth.\n",
    "\n",
    " - constructor: initialize the model to ZH-L16C w/ 5-minute compartment,\n",
    "   - diving air,\n",
    "   - at sea level,\n",
    "   - using $ RQ = 0.9 $\n",
    "   - and fixed gradient factor 0.85\n",
    " \n",
    " - `segment()` takes the new depth and time and updates gas loadings.\n",
    " \n",
    " This also keeps track of the ceiling together with $ P_t $\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pprint\n",
    "import re\n",
    "\n",
    "import diyzhl\n",
    "\n",
    "#\n",
    "#\n",
    "class dive( object ) :\n",
    "    \n",
    "    # ctor\n",
    "    #\n",
    "    def __init__( self, verbose = False ) :\n",
    "\n",
    "        self._verbose = bool( verbose )\n",
    "        self._timepat = re.compile( r\"(?:(\\d{1,2}):)?(\\d{1,2}):(\\d{1,2})\" )\n",
    "        \n",
    "# air, sea level, USN RQ. S is surface pressure (const), P is current pressure (var)\n",
    "        self._T = 0\n",
    "        self._S = 1.0\n",
    "        self._P = self._S\n",
    "        self._Q = 0.79\n",
    "        self._RQ = 0.9\n",
    "        self._GFHi = 0.85\n",
    "        self._TCs = []\n",
    "\n",
    "# starting Pt (same for all TCs)\n",
    "        sp = diyzhl.palv( Pamb = self._P, Q = self._Q, RQ = self._RQ )\n",
    "\n",
    "# use ZH-L16Cb (skip over 4-minute TC)\n",
    "        for tc in diyzhl.ZHL16N.keys() :\n",
    "            if tc == 1 : continue\n",
    "            self._TCs.append( { \n",
    "                \"t\" : diyzhl.ZHL16N[tc][\"t\"], \n",
    "                \"a\" : diyzhl.ZHL16N[tc][\"a\"][\"C\"],\n",
    "                \"b\" : diyzhl.ZHL16N[tc][\"b\"],\n",
    "                \"P\" : sp\n",
    "            } )\n",
    "\n",
    "# init. ceiling\n",
    "        for i in range( len( self._TCs ) ) :\n",
    "            self._TCs[i][\"C\"] = diyzhl.buhlmann( Pn = self._TCs[i][\"P\"],\n",
    "                                         an = self._TCs[i][\"a\"],\n",
    "                                         bn = self._TCs[i][\"b\"],\n",
    "                                         gf = self._GFHi )\n",
    "            \n",
    "        if self._verbose :\n",
    "            pprint.pprint( self._TCs )\n",
    "            \n",
    "    # helpers for plotting\n",
    "    # (could actually do this in a more \"pythonic\" way but this is more obvious)\n",
    "    #\n",
    "    @property\n",
    "    def compartments( self ) :\n",
    "        rc = []\n",
    "        for i in range( len( self._TCs ) ) :\n",
    "            rc.append( self._TCs[i][\"t\"] )\n",
    "        return rc\n",
    "    \n",
    "    @property\n",
    "    def loadings( self ) :\n",
    "        rc = []\n",
    "        for i in range( len( self._TCs ) ) :\n",
    "            rc.append( self._TCs[i][\"P\"] )\n",
    "        return rc\n",
    "\n",
    "    # helper function: takes human-readable time string like \"1:30\" and returns minutes: 1.5\n",
    "    #\n",
    "    def _time( self, t = \"0:0\" ) :\n",
    "        if t is None : return 0\n",
    "        m = self._timepat.search( str( t ).strip() )\n",
    "        if not m : raise Exception( \"Invalid time string %s\" % (t,) )\n",
    "        rc = 0.0\n",
    "        if m.group( 1 ) is not None :\n",
    "            rc = float( m.group( 1 ) ) * 60.0\n",
    "        rc += float( m.group( 2 ) )\n",
    "        rc += float( m.group( 3 ) ) / 60.0\n",
    "        return round( rc, 1 )\n",
    "\n",
    "    # newdepth is new depth in 0.1 bar\n",
    "    # timestr is time as [hours:]minutes:seconds string. *it is the total elapsed* time\n",
    "    #\n",
    "    def segment( self, newdepth = 0.0, newtimestr = \"1:0\" ) :\n",
    "        assert float( newdepth ) >= 0.0\n",
    "        if float( newdepth ) == 0.0 :\n",
    "            newP = self._P\n",
    "        else :\n",
    "            newP = round( self._S + float( newdepth ) / 10, 1 )\n",
    "        t = self._time( newtimestr ) - self._T\n",
    "    \n",
    "        for i in range( len( self._TCs ) ) :\n",
    "            p = diyzhl.schreiner( Pi = self._TCs[i][\"P\"], \n",
    "                          Palv = diyzhl.palv( Pamb = self._P, Q = self._Q, RQ = self._RQ ), \n",
    "                          t = t, \n",
    "                          R = diyzhl.arr( d0 = self._P, dt = newP, t = t, Q = self._Q ),\n",
    "                          k = diyzhl.kay( Th = self._TCs[i][\"t\"] ) )\n",
    "            self._TCs[i][\"P\"] = p\n",
    "            self._TCs[i][\"C\"] = diyzhl.buhlmann( Pn = self._TCs[i][\"P\"],\n",
    "                                         an = self._TCs[i][\"a\"],\n",
    "                                         bn = self._TCs[i][\"b\"],\n",
    "                                         gf = self._GFHi )\n",
    "\n",
    "        self._P = newP\n",
    "        self._T += t\n",
    "    \n",
    "        if self._verbose :\n",
    "            sys.stdout.write( \"* At time %f, P %f:\\n\" % (self._T, self._P,) )\n",
    "            pprint.pprint( self._TCs )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animate\n",
    "\n",
    "The profile is from a real dive, simplified. \n",
    "\n",
    "(TODO: change bar colour when TC has a ceiling)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAD6tJREFUeJzt3X+MZXV9xvH3I7v+Qiu2jJUui6sVpUpUcIpY04aKJoiEbVKaYFqKVrOJEcXG/gBNMOWPBttGW4uRbIUCloAGUbcWq1Rt1aSiw5Zfy2rdqpURKoMoSLXabT/9454109m73DMzd+be/e77lUzmnnO+c+6zszvP/c6555xNVSFJasujJh1AkjR+lrskNchyl6QGWe6S1CDLXZIaZLlLUoNGlnuSxyb5YpLbkuxK8kdDxjwmyQeS7Elyc5ItaxFWktRPn5n7j4CXVtXzgRcApyU5ecmY1wLfrapnAu8C3jHemJKk5RhZ7jXwcLe4sftYeuXTVuCq7vH1wKlJMraUkqRl2dBnUJLDgFuAZwLvqaqblwzZBNwNUFV7kzwI/Axw/5L9bAO2ARx++OEvPO6441aXXpIOMbfccsv9VTUzalyvcq+q/wFekOQI4MNJjq+qOxcNGTZL3+++BlW1HdgOMDs7W3Nzc32eXpLUSfLvfcYt62yZqvoe8I/AaUs2zQObuyfeADwJeGA5+5YkjU+fs2Vmuhk7SR4HvAz48pJhO4Bzu8dnAZ8u70gmSRPT57DMUcBV3XH3RwEfrKqPJbkYmKuqHcDlwPuT7GEwYz97zRJLkkYaWe5VdTtwwpD1Fy16/F/Ab4w3miRppbxCVZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lq0MhyT7I5yWeS7E6yK8n5Q8ackuTBJLd2HxetTVxJUh8beozZC7ylqnYmeSJwS5KbququJeM+V1VnjD+iJGm5Rs7cq+reqtrZPf4+sBvYtNbBJEkrt6xj7km2ACcANw/Z/OIktyX5eJLnjiGbJGmF+hyWASDJE4APAW+uqoeWbN4JPK2qHk5yOvAR4Ngh+9gGbAM45phjVhxakvTIes3ck2xkUOzXVNUNS7dX1UNV9XD3+EZgY5Ijh4zbXlWzVTU7MzOzyuiSpAPpc7ZMgMuB3VX1zgOMeWo3jiQndfv9zjiDSpL663NY5iXAOcAdSW7t1r0VOAagqi4DzgJen2Qv8EPg7KqqNcgrSephZLlX1eeBjBhzKXDpuEJJklbHK1QlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUoJHlnmRzks8k2Z1kV5Lzh4xJkncn2ZPk9iQnrk1cSVIfG3qM2Qu8pap2JnkicEuSm6rqrkVjXgEc2328CHhv91mSNAEjZ+5VdW9V7ewefx/YDWxaMmwrcHUNfAE4IslRY08rSeplWcfck2wBTgBuXrJpE3D3ouV59n8BIMm2JHNJ5hYWFpaXVJLUW+9yT/IE4EPAm6vqoaWbh3xJ7beiantVzVbV7MzMzPKSSpJ661XuSTYyKPZrquqGIUPmgc2Llo8G7ll9PEnSSvQ5WybA5cDuqnrnAYbtAH67O2vmZODBqrp3jDklScvQ52yZlwDnAHckubVb91bgGICqugy4ETgd2AP8AHjN+KNKkvoaWe5V9XmGH1NfPKaAN4wrlCRpdbxCVZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoM2TDpAa7Zc8Hcr/tpvXPLKMSaRdCgbWe5JrgDOAO6rquOHbD8F+Cjw9W7VDVV18ThDHqpW80IBvlhIh7I+M/crgUuBqx9hzOeq6oyxJJIkrdrIY+5V9VnggXXIIkkak3G9ofriJLcl+XiS5x5oUJJtSeaSzC0sLIzpqSVJS42j3HcCT6uq5wN/CXzkQAOrantVzVbV7MzMzBieWpI0zKrLvaoeqqqHu8c3AhuTHLnqZJKkFVv1qZBJngp8u6oqyUkMXjC+s+pk6+hQOH3RM2+kQ0ufUyGvBU4BjkwyD7wd2AhQVZcBZwGvT7IX+CFwdlXVmiWWJI00styr6lUjtl/K4FRJSdKU8PYDktQgy12SGmS5S1KDLHdJapDlLkkN8pa/WpFD4doA6WDmzF2SGmS5S1KDLHdJatBBeczd+6RI0iNz5i5JDbLcJalBlrskNchyl6QGHZRvqKotXhAljZ8zd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGuR57mqKN5WTBpy5S1KDLHdJapDlLkkNstwlqUEjyz3JFUnuS3LnAbYnybuT7Elye5ITxx9TkrQcfc6WuRK4FLj6ANtfARzbfbwIeG/3WTroecdKHaxGztyr6rPAA48wZCtwdQ18ATgiyVHjCihJWr5xHHPfBNy9aHm+W7efJNuSzCWZW1hYGMNTS5KGGUe5Z8i6GjawqrZX1WxVzc7MzIzhqSVJw4yj3OeBzYuWjwbuGcN+JUkrNI7bD+wAzktyHYM3Uh+sqnvHsF+pKb45q/U0styTXAucAhyZZB54O7ARoKouA24ETgf2AD8AXrNWYSVJ/Yws96p61YjtBbxhbIkkSavmXSGlg5CHeDSKtx+QpAZZ7pLUIA/LSIc4/4OTNjlzl6QGOXOXNFa+2TsdLHdJU8sXipXzsIwkNciZu6RDwqH2xrHlLknLdDC8UHhYRpIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhrUq9yTnJbkK0n2JLlgyPZXJ1lIcmv38brxR5Uk9TXyv9lLchjwHuDlwDzwpSQ7ququJUM/UFXnrUFGSdIy9Zm5nwTsqaqvVdWPgeuArWsbS5K0Gn3KfRNw96Ll+W7dUr+e5PYk1yfZPJZ0kqQV6VPuGbKuliz/LbClqp4H/ANw1dAdJduSzCWZW1hYWF5SSVJvfcp9Hlg8Ez8auGfxgKr6TlX9qFv8K+CFw3ZUVduraraqZmdmZlaSV5LUQ59y/xJwbJKnJ3k0cDawY/GAJEctWjwT2D2+iJKk5Rp5tkxV7U1yHvAJ4DDgiqraleRiYK6qdgBvSnImsBd4AHj1GmaWJI0wstwBqupG4MYl6y5a9PhC4MLxRpMkrZRXqEpSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBvco9yWlJvpJkT5ILhmx/TJIPdNtvTrJl3EElSf2NLPckhwHvAV4BPAd4VZLnLBn2WuC7VfVM4F3AO8YdVJLUX5+Z+0nAnqr6WlX9GLgO2LpkzFbgqu7x9cCpSTK+mJKk5UhVPfKA5CzgtKp6Xbd8DvCiqjpv0Zg7uzHz3fK/dWPuX7KvbcC2bvHZwFfG9QdZ4kjg/pGjJmNas01rLpjebNOaC6Y327TmgunNtjTX06pqZtQXbeix42Ez8KWvCH3GUFXbge09nnNVksxV1exaP89KTGu2ac0F05ttWnPB9Gab1lwwvdlWmqvPYZl5YPOi5aOBew40JskG4EnAA8sNI0kajz7l/iXg2CRPT/Jo4Gxgx5IxO4Bzu8dnAZ+uUcd7JElrZuRhmaram+Q84BPAYcAVVbUrycXAXFXtAC4H3p9kD4MZ+9lrGbqHNT/0swrTmm1ac8H0ZpvWXDC92aY1F0xvthXlGvmGqiTp4OMVqpLUIMtdkhrUXLmPulXCJCTZnOQzSXYn2ZXk/ElnWizJYUn+JcnHJp1lsSRHJLk+yZe7792LJ51pnyS/2/1d3pnk2iSPnWCWK5Lc111vsm/dTye5KclXu89PnpJcf9r9fd6e5MNJjljvXAfKtmjb7yWpJEdOS64kb+x6bVeSP+mzr6bKveetEiZhL/CWqvoF4GTgDVOSa5/zgd2TDjHEXwB/X1XHAc9nSjIm2QS8CZitquMZnGgwyZMIrgROW7LuAuBTVXUs8Klueb1dyf65bgKOr6rnAf8KXLjeoTpXsn82kmwGXg58c70Dda5kSa4kv8rgLgDPq6rnAn/WZ0dNlTv9bpWw7qrq3qra2T3+PoOS2jTZVANJjgZeCbxv0lkWS/JTwK8wOBOLqvpxVX1vsqn+nw3A47rrOh7P/td+rJuq+iz7X1ey+JYgVwG/tq6hGJ6rqj5ZVXu7xS8wuG5m3R3gewaDe2P9AUMuwlwPB8j1euCSqvpRN+a+Pvtqrdw3AXcvWp5nSkp0n+6OmScAN082yU/8OYN/zP876SBLPANYAP66O2T0viSHTzoUQFV9i8Hs6ZvAvcCDVfXJyabaz89W1b0wmFwAT5lwnmF+B/j4pEPsk+RM4FtVdduksyzxLOCXuzvu/lOSX+zzRa2Ve6/bIExKkicAHwLeXFUPTUGeM4D7quqWSWcZYgNwIvDeqjoB+E8mc2hhP93x663A04GfAw5P8luTTXVwSfI2Bocrr5l0FoAkjwfeBlw06SxDbACezOCQ7u8DH+xzY8bWyr3PrRImIslGBsV+TVXdMOk8nZcAZyb5BoNDWC9N8jeTjfQT88B8Ve37Ded6BmU/DV4GfL2qFqrqv4EbgF+acKalvp3kKIDuc69f5ddDknOBM4DfnKIr2X+ewYv1bd3Pw9HAziRPnWiqgXnghhr4IoPfske+2dtaufe5VcK6615lLwd2V9U7J51nn6q6sKqOrqotDL5Xn66qqZiBVtV/AHcneXa36lTgrglGWuybwMlJHt/93Z7KlLzZu8jiW4KcC3x0gll+IslpwB8CZ1bVDyadZ5+quqOqnlJVW7qfh3ngxO7f4aR9BHgpQJJnAY+mx90rmyr37o2afbdK2A18sKp2TTYVMJghn8NgZnxr93H6pEMdBN4IXJPkduAFwB9POA8A3W8T1wM7gTsY/BxN7NL1JNcC/ww8O8l8ktcClwAvT/JVBmd/XDIluS4Fngjc1P0cXLbeuR4h28QdINcVwDO60yOvA87t8xuPtx+QpAY1NXOXJA1Y7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalB/wfAWUtgJJpIkwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from IPython import display\n",
    "\n",
    "d = dive( verbose = False )\n",
    "idx = np.arange( len( d.compartments ) )\n",
    "plt.xticks( idx, d.compartments, rotation = 30 )\n",
    "       \n",
    "\n",
    "PROFILE = [(24,\"1:40\"),(30,\"6:10\"),(29,\"10:20\"),(23,\"16:20\"),(15,\"21:30\"),(10,\"25:30\"),(3,\"40:0\")]\n",
    "\n",
    "for i in PROFILE : \n",
    "    display.clear_output( wait = True )\n",
    "    d.segment( newdepth = i[0], newtimestr = i[1] )\n",
    "    plt.ylim( 0, 3.0 )\n",
    "    plt.bar( idx, d.loadings )\n",
    "    plt.show()\n",
    "    display.display( plt.gcf() )\n",
    "    time.sleep( 1.5 )\n",
    "\n",
    "\n"
   ]
  },
  {
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}