diyzhl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DIY ZHL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References etc.\n",
    "\n",
    "  * Mark Powell's \"Deco for Divers\" has all the background information in a succinct and easy to read form. The book should be available on interlibrary loan, ask your librarian.\n",
    "  \n",
    "  * Erik Baker's papers: \"Decolessons\", \"Understanding M-values\", \"Clearing up the confusion about Deep Stops\", are more technical and tend to omit much of the background info. PDFs are freely available on-line from several sources.\n",
    "  \n",
    "  * Buhlmann's \"Decompression - Decompression sickness\".\n",
    "  \n",
    "  * Several open-source implementations are out there, most notable Subsurface software, but you have to read C.\n",
    "  \n",
    "  * Plenty of other on-line sources, e.g. Stuart Morrison's \"DIY Decompresion\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Primer\n",
    "\n",
    "### Dissolved gas\n",
    "\n",
    "Gas is driven in and out of tissue by the difference in pressure: \"outside\" aka ambient $ P_{amb} $ &nbsp; and \"inside\" aka dissolved gas pressure, $ P_{t} $. The rate is assumed to be $ ln(2) $\n",
    "\n",
    "**Partial pressure** refers to pressures of individual gases in the breathing mix: oxygen, nitrogen, etc.\n",
    "\n",
    "### Tissue compartments\n",
    "\n",
    "To simulate real tissue the model uses hypothetical \"compartments\". They are defined by the time $ T_{ h } $ it takes for the pressure difference (above) to get to -- since we're using log 2 -- half (or double, depending on whether the tissues in on- or off-gassing) its intial value. That time is callled ``half time``. It is usually given in minutes.\n",
    "\n",
    "The body is modelled as a set of \"tissue compartments\" with different half-times. They have nothing to do with real tissues, but to illustrate the concept it is usually said that blood is \"fast\": gas gets in and out very quickly, whereas bone is very \"slow\".\n",
    "\n",
    "### M-values\n",
    "\n",
    "The greater the difference between $ P_{amb} $ &nbsp; and $ P_t $, the faster the gas flows. There is, however, such thing as too fast. To recreational divers \"too fast\" is only relevant on the way up.\n",
    "\n",
    "The model defines \"maximum allowed pressure difference\" value, ``M-value`` for short. The model does not attempt to quantify the risk, it only works in binary \"above the M-value\" (bad) or \"below the M-value\" (OK) terms.\n",
    "\n",
    "Haldane originally set his M-value as a fixed ratio (2:1). Workman replaced it with a function: a straight line $ M = \\Delta M + M_0 $ with greater overpressure tolerated at greater depths."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Water vapour etc. in alveolae\n",
    "\n",
    "Breathing gas in respiratory passages is diluted with water vapour and expelled carbon dioxide. The fraction of water vapour and carbon dioxide is subtracted from the partial pressure of inert gas.\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "P_{alv} = ( P_{amb} - P_{H_{2}O} + ( 1 - RQ ) / RQ * P_{CO_{2}} ) Q\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "where\n",
    "  * $ P_{ CO_2 } $ is pressure of carbon dioxide, 40 mm Hg or 0.0534 bar,\n",
    "  * $ RQ $ is respiratory quotient,\n",
    "  * $ P_{ H_{2}O } $ is pressure of water vapour, 47 mm Hg or 0.0627 bar,\n",
    "  * $ Q $ is the fraction of inert gas,\n",
    "  * $ P_{amb} $ is ambient pressure.\n",
    "\n",
    "$ RQ $ is a ratio of oxygen consumed to carbon dioxide produced in metabolism, it depends on nutrition and exertion and is in the range 0.7..1. Smaller number results in greater calculated gas loading. Buhlmann uses 1 whereas Shreiner uses 0.8 and USN uses 0.9. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Schreiner equation to calculate inert gas pressure in a tissue compartment\n",
    "\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "P_{t}(t) = P_{alv0} + R * (t - 1 / k) - (P_{alv0} - P_{t0} - R / k) * e^{-k * t}\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "*where*\n",
    "\n",
    "  * $ P_{t}(t) $ : partial pressure of inert gas in the tissue (bar) at time $ t $\n",
    "  \n",
    "  * $ P_{alv0} $ : initial alveolar partial pressure of the gas in the breathing mix (bar)\n",
    "  \n",
    "  * $ P_{t0} $ : initial partial pressure of the gas in the tissue (bar) (at time $ t0 $)\n",
    "  \n",
    "  * $ R $ is rate of change in inert gas pressure in the breathing mix, bar/min. This is rate of change in ambient pressure, $ \\Delta P $ times fraction of inert gas. $ \\Delta P $ is negative for ascent, positive for descent, zero for constant depth.\n",
    "\n",
    "  * $ t $ : time of exposure in minutes\n",
    "\n",
    "  * $ k $ : gas decay constant for a tissue compartment: $ k = \\frac  {\\ln(2)} {T_{h}} $\n",
    "\n",
    "    * $ \\quad T_{h} $ : gas \"half-life\" time for tissue compartment in minutes\n",
    "\n",
    "  * $ e $ : 2.7182818284590452353602874713527 and counting\n",
    "\n",
    "Note that for constant depth this reduces to a \"short form\" given in many sources, the code below just uses the full form with $ R = 0 $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## M-values\n",
    "\n",
    "### Workman\n",
    "\n",
    "Workman-style M-value $ M $ is defined in terms of $ \\Delta M $ and $ M_{ 0 } $, where\n",
    "\n",
    "  * $ M_{0} $ is the surfacing M-value and\n",
    "  * $ \\Delta M $ is the slope of M-value line: multiply by depth to get tha actual number at that depth.\n",
    "\n",
    "$ M = \\Delta M * Depth + M_{0} $\n",
    "\n",
    "\n",
    "$ Tolerated Depth = ( P - M_{0} ) / \\Delta M $\n",
    "\n",
    "where Depth is measured as pressure from sea level and Tolerated Depth is how high you can ascend without exceeding the M-value.\n",
    "\n",
    "### Buhlmann\n",
    "\n",
    "In Buhlmann model the M-values are expresed in terms of two coeficients: $ a $ and $ b $. In addition, it is using absolute instead od sea-level pressure for diving at altitude, and comes with separate numbers for different inert gases. Buhlmann's model calls \"Tolerated Depth\" an \"ascent ceiling\".\n",
    "\n",
    "### Conversion\n",
    "\n",
    "  * Workman to Buhlmann:\n",
    "  \n",
    "   * $ a = M_{0} - \\Delta M * P_{ \\text{ amb. (surface at sea level) }} $ \n",
    "   \n",
    "   * $ b = 1 / \\Delta M $  \n",
    "\n",
    "  * Buhlmann to Workman:\n",
    "\n",
    "   * $ \\Delta M = 1 / b $\n",
    "   \n",
    "   * $ M_{0} = a + \\frac { P_{ \\text{ amb. (surface at sea level) }}} { b } $\n",
    "   \n",
    "**Note** that this only really works with $ \\Delta M $ greater than 1. See Erik Baker's \"M-values\" for the why.\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## No-stop time\n",
    "\n",
    "### by (variant of) Schreiner equation\n",
    "\n",
    "This needs Workman-style M-values:\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "P_{no-deco} \\ = ( M_{0} - P_{alv0} - R * (t_{asc} - 1/k) ) * e ^ {k * t_{asc}} \\ + P_{alv0} - R / k\n",
    "\\end{align}\n",
    "$\n",
    "  * $ P_{no-deco} $ &nbsp; is the pressure at which ascent must be started (bar),\n",
    "  * $ M_{0} $ is the surfacing M-value,\n",
    "  * $ P_{alv0} $ is alveolar partial pressure (bar),\n",
    "  * $ t_{asc} $ is the time needed for ascent (minutes): depth divided by speed.\n",
    "  * $ R, k, e $ are the same as in Schreiner's equation above.\n",
    "  \n",
    "You have to figure out how long it'll take to get from the current inspired pressure $ P_{t} $ to $ P_{no-deco} $ &nbsp; : run `schreiner()` in a loop adding a minute of exposure on each iteration until you hit $ P_{no-deco} $ &nbsp; (or do a fancy binary search from, say, 99 minutes down, etc.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### by Buhlmann equation\n",
    "\n",
    "Buhlmann equation calculates Tolerated Depth in absolute pressure and calls it \"safe ascent ceiling\":\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "P_{amb.tol.} = ( P_{t} - a ) * b\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "where\n",
    "  * $ P_{amb.tol.} $ &nbsp; is the ascent ceiling,\n",
    "  * $ P_{t} $ is current inert gas pressure in the tisue, and\n",
    "  * $ a $ and $ b $ are Buhlmann's coefficients.\n",
    "\n",
    "When the calculated $ P_{amb.tol.} $ &nbsp; is above 1, i.e. under the water surface (assuming surface pressure is 1 bar of course), that means we have a deco ceiling and must stop and off-gas there.\n",
    "\n",
    "Buhlmann equation modified by Erik Baker's Gradient Factors is\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "P_{amb.tol.} = ( P_{t} - a * GF ) / ( GF / b + 1.0 - GF )\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "where $ GF $ is the \"current\" gradient factor. \n",
    "\n",
    "(Without the gradient factor: $ GF = 1 $, that reduces to $ ( P_{t} - a * 1 ) / ( 1/b + 1 - 1 ) => ( P_{t} - a ) / ( 1/b ) =>  ( P_{t} - a ) * b $: Buhlmann's original.)\n",
    "\n",
    "Gradient Factors plural are a pair of numbers: $ GFLo $ that applies at first decompression stop and $ GFHi $ that applies at the last stop: the surface. The \"current\" $ GF $ at any given depth between those two points is calculated from a simple proportion. **Note** that this only works when you have a relatively deep decompression stop.\n",
    "\n",
    "### You can calculate no-stop time\n",
    "by running `shreiner()` and `buhlmann()` in a loop, adding a minute of exposure at each iteration, until the latter produces a ceiling. **Note** however, that neither this method nor the equation above factor in off-gassing that will take place during ascent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Helium\n",
    "\n",
    "With multiple inert gases, each has to be tracked separately by Schreiner equation.\n",
    "\n",
    "For Buhlmann's formula you add up the partial pressures and proportionally adjust the $ a $ and $ b $:\n",
    "\n",
    "$ P_{t} = P_{iN_2} + P_{iHe} $\n",
    "\n",
    "$ a = (a_{N_2} * P_{N_2} + a_{He} * P_{He}) / P_{t} $\n",
    "\n",
    "$ b = (b_{N_2} * P_{N_2} + b_{He} * P_{He}) / P_{t} $\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decompression stop\n",
    "\n",
    "... is ascent ceiling rounded to next 3 metres down.\n",
    "\n",
    "The best I could find, in Buhlmann's \"Decompression\", is in open water it's hard to maintain stops at smaller intervals. Buhlmann suggests using smaller intervals or even \"ride the M-value line\" continuously for in-chamber decompression there. (**Note** this assumes no gas switching.)\n",
    "\n",
    "Once you're at deco stop, run `schreiner()` and `buhlmann()` in a loop adding a minute on each iteration, as usual, until your next ceiling is at your next stop: 3 metres above you current depth if using the 3 m spacing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tissue compartments\n",
    "\n",
    "Tables are coded as python dicts: keys for each compartment are ``t``, ``a``, and ``b`` for half-time and the coefficients resp. E.g. ZH-L12's $ b $ coefficient for nitrogen in TC 3 is ``ZHL12N[3][\"b\"]``\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ZH-L16\n",
    "\n",
    "**Note** that you have two options for the 1st (fastest) tissue compartment: 4 minutes or 5. Pick one. The 5-minute TC is usually called \"1b\" and is keyed as \"1.1\" in the table below.\n",
    "\n",
    "The ''a'' key in the nitrogen tables is a dict itself with three sub-keys: ''A'', ''B'', adn ''C'' for the three versions of ZH-L16 resp. Helium numbers seem to be unobtainable for ''A'' or ''C'', nor for the 4-minute tissue compartment. I.e. if you want to use this for trimix, ZH-L16B with 5-minute TC is your only option.\n",
    "\n",
    "M-values for \"-A\" were calculated as\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "a = 2 * T_{h} ^ {- \\frac{1}{3}}\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "b = 1.005 - T_{h} ^ {- \\frac{1}{2}}\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "\n",
    "  * $ T_{h} $ is tissue compartment half-time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workman\n",
    "\n",
    "Not useful as is, would need to be recalculated for $ \\Delta M $ instead of $ M $ first, then converted to Buhlmann's `a` and `b` coefficients."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DSAT\n",
    "\n",
    "Convert to Buhlmann's `a` and `b` using $ \\Delta M = 1 $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code\n",
    "\n",
    "Run the cell below to load python code (assuming you have the `diyzhl.py` in the same directory as the notebook file).\n",
    "\n",
    "The functions are written to closely follow the above formulae. That includes calculating constants on the fly, like e.g. $ K $ or $ P_{alv} $ &nbsp; -- you'd normally use pre-calculated numbers for your chosen parameters.\n",
    "\n",
    "About the only added features are \n",
    " * some of them save intermediate results in variables: mainly because I wanted to play with rounding at different steps. Note that although I round the numbers to 4 decimal digits, they're probably only meaningful to 2 at best.\n",
    " \n",
    " * Typecasts to `float` are there so that passing in not-a-number will throw an exception early.      \n",
    " \n",
    " * Similarly, `asserts` are used to catch division by zero early."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load diyzhl.py\n",
    "#!/usr/bin/python -u\n",
    "#\n",
    "# (K) Copy Rites Reversed: reuse what you like (but give credit)\n",
    "#\n",
    "# Credits:\n",
    "#\n",
    "# Mark Powell's \"Deco for Divers\"\n",
    "# Erik Baker's papers: \"Decolessons\", \"Understanding M-values\", and \"Clearing up the confusion about Deep Stops\" in particular\n",
    "# Buhlmann's \"Decompression - Decompression Sickness\", English edition\n",
    "# Several open-source implementations, most notably Subsurface software (and people, Robert in particular)\n",
    "# Plenty of other on-line sources, e.g. Stuart Morrison's \"DIY Decompresion\"\n",
    "#\n",
    "# The goal here is \"by the book\" implementation to use for learning this stuff.\n",
    "# There's a lot of things that can be done differently, more efficeint, and so on... and\n",
    "#  that's not what this code is for.\n",
    "#\n",
    "\n",
    "import sys\n",
    "import math\n",
    "\n",
    "# ZH-L12 from \"Decompression\"\n",
    "#\n",
    "ZHL12N = {\n",
    "    1 :  { \"t\" : 2.65,  \"a\" : 2.2,   \"b\" : 0.82 },\n",
    "    2 :  { \"t\" : 7.94,  \"a\" : 1.5,   \"b\" : 0.82 },\n",
    "    3 :  { \"t\" : 12.2,  \"a\" : 1.08,  \"b\" : 0.825 },\n",
    "    4 :  { \"t\" : 18.5,  \"a\" : 0.9,   \"b\" : 0.835 },\n",
    "    5 :  { \"t\" : 26.5,  \"a\" : 0.75,  \"b\" : 0.845 },\n",
    "    6 :  { \"t\" : 37.0,  \"a\" : 0.58,  \"b\" : 0.86 },\n",
    "    7 :  { \"t\" : 53.0,  \"a\" : 0.47,  \"b\" : 0.87 },\n",
    "    8 :  { \"t\" : 79.0,  \"a\" : 0.45,  \"b\" : 0.89 },\n",
    "    9 :  { \"t\" : 114.0, \"a\" : 0.45,  \"b\" : 0.89 },\n",
    "    10 : { \"t\" : 146.0, \"a\" : 0.455, \"b\" : 0.934 },\n",
    "    11 : { \"t\" : 185.0, \"a\" : 0.455, \"b\" : 0.934 },\n",
    "    12 : { \"t\" : 238.0, \"a\" : 0.38,  \"b\" : 0.944 },\n",
    "    13 : { \"t\" : 304.0, \"a\" : 0.255, \"b\" : 0.962 },\n",
    "    14 : { \"t\" : 397.0, \"a\" : 0.255, \"b\" : 0.962 },\n",
    "    15 : { \"t\" : 503.0, \"a\" : 0.255, \"b\" : 0.962 },\n",
    "    16 : { \"t\" : 635.0, \"a\" : 0.255, \"b\" : 0.962 }\n",
    "}\n",
    "\n",
    "ZHL12He = {\n",
    "    1 :  { \"t\" : 1.0,   \"a\" : 2.2,   \"b\" : 0.82 },\n",
    "    2 :  { \"t\" : 3.0,   \"a\" : 1.5,   \"b\" : 0.82 },\n",
    "    3 :  { \"t\" : 4.6,   \"a\" : 1.08,  \"b\" : 0.825 },\n",
    "    4 :  { \"t\" : 7.0,   \"a\" : 0.9,   \"b\" : 0.835 },\n",
    "    5 :  { \"t\" : 10.0,  \"a\" : 0.75,  \"b\" : 0.845 },\n",
    "    6 :  { \"t\" : 14.0,  \"a\" : 0.58,  \"b\" : 0.86 },\n",
    "    7 :  { \"t\" : 20.0,  \"a\" : 0.47,  \"b\" : 0.87 },\n",
    "    8 :  { \"t\" : 30.0,  \"a\" : 0.45,  \"b\" : 0.89 },\n",
    "    9 :  { \"t\" : 43.0,  \"a\" : 0.45,  \"b\" : 0.89 },\n",
    "    10 : { \"t\" : 55.0,  \"a\" : 0.515, \"b\" : 0.926 },\n",
    "    11 : { \"t\" : 70.0,  \"a\" : 0.515, \"b\" : 0.926 },\n",
    "    12 : { \"t\" : 90.0,  \"a\" : 0.515, \"b\" : 0.926 },\n",
    "    13 : { \"t\" : 115.0, \"a\" : 0.515, \"b\" : 0.926 },\n",
    "    14 : { \"t\" : 150.0, \"a\" : 0.515, \"b\" : 0.926 },\n",
    "    15 : { \"t\" : 190.0, \"a\" : 0.515, \"b\" : 0.926 },\n",
    "    16 : { \"t\" : 240.0, \"a\" : 0.515, \"b\" : 0.926 },\n",
    "}\n",
    "\n",
    "# ZH_L16: several sources incl. a photo of a page from Tauchmedizin @\n",
    "# http://www.nigelhewitt.co.uk/stuff/aab.jpg\n",
    "# It appears nobody has Helium numbers for \"-A\" and \"-C\", nor for the 4-minute TC\n",
    "# 5-minute TC is keyed as 1.1\n",
    "#\n",
    "ZHL16N = {\n",
    "    1   : { \"t\" : 4.0,   \"b\" : 0.505,  \"a\" : { \"A\" : 1.2599, \"B\" : 1.2599, \"C\" : 1.2599 } },\n",
    "    1.1 : { \"t\" : 5.0,   \"b\" : 0.5578, \"a\" : { \"A\" : 1.1696, \"B\" : 1.1696, \"C\" : 1.1696 } },\n",
    "    2   : { \"t\" : 8.0,   \"b\" : 0.6514, \"a\" : { \"A\" : 1.0,    \"B\" : 1.0,    \"C\" : 1.0 } },\n",
    "    3   : { \"t\" : 12.5,  \"b\" : 0.7222, \"a\" : { \"A\" : 0.8618, \"B\" : 0.8618, \"C\" : 0.8618 } },\n",
    "    4   : { \"t\" : 18.5,  \"b\" : 0.7825, \"a\" : { \"A\" : 0.7562, \"B\" : 0.7562, \"C\" : 0.7562 } },\n",
    "    5   : { \"t\" : 27.0,  \"b\" : 0.8126, \"a\" : { \"A\" : 0.6667, \"B\" : 0.6667, \"C\" : 0.62 } },\n",
    "    6   : { \"t\" : 38.3,  \"b\" : 0.8434, \"a\" : { \"A\" : 0.5933, \"B\" : 0.56,   \"C\" : 0.5043 } },\n",
    "    7   : { \"t\" : 54.3,  \"b\" : 0.8693, \"a\" : { \"A\" : 0.5282, \"B\" : 0.4947, \"C\" : 0.441 } },\n",
    "    8   : { \"t\" : 77.0,  \"b\" : 0.891,  \"a\" : { \"A\" : 0.4701, \"B\" : 0.45,   \"C\" : 0.4 } },\n",
    "    9   : { \"t\" : 109.0, \"b\" : 0.9092, \"a\" : { \"A\" : 0.4187, \"B\" : 0.4187, \"C\" : 0.375 } },\n",
    "    10  : { \"t\" : 146.0, \"b\" : 0.9222, \"a\" : { \"A\" : 0.3798, \"B\" : 0.3798, \"C\" : 0.35 } },\n",
    "    11  : { \"t\" : 187.0, \"b\" : 0.9319, \"a\" : { \"A\" : 0.3497, \"B\" : 0.3497, \"C\" : 0.3295 } },\n",
    "    12  : { \"t\" : 239.0, \"b\" : 0.9403, \"a\" : { \"A\" : 0.3223, \"B\" : 0.3223, \"C\" : 0.3065 } },\n",
    "    13  : { \"t\" : 305.0, \"b\" : 0.9477, \"a\" : { \"A\" : 0.2971, \"B\" : 0.285,  \"C\" : 0.2835 } },\n",
    "    14  : { \"t\" : 390.0, \"b\" : 0.9544, \"a\" : { \"A\" : 0.2737, \"B\" : 0.2737, \"C\" : 0.261 } },\n",
    "    15  : { \"t\" : 498.0, \"b\" : 0.9602, \"a\" : { \"A\" : 0.2523, \"B\" : 0.2523, \"C\" : 0.248 } },\n",
    "    16  : { \"t\" : 635.0, \"b\" : 0.9653, \"a\" : { \"A\" : 0.2327, \"B\" : 0.2327, \"C\" : 0.2327 } }\n",
    "}\n",
    "\n",
    "ZHL16He = {\n",
    "    1.1 : { \"t\" : 1.88,   \"a\" : { \"B\" : 1.6189 }, \"b\" : 0.477 },\n",
    "    2 :   { \"t\" : 3.02,   \"a\" : { \"B\" : 1.383 },  \"b\" : 0.5747 },\n",
    "    3 :   { \"t\" : 4.72,   \"a\" : { \"B\" : 1.1919 }, \"b\" : 0.6527 },\n",
    "    4 :   { \"t\" : 6.99,   \"a\" : { \"B\" : 1.0458 }, \"b\" : 0.7223 },\n",
    "    5 :   { \"t\" : 10.21,  \"a\" : { \"B\" : 0.922 },  \"b\" : 0.7582 },\n",
    "    6 :   { \"t\" : 14.48,  \"a\" : { \"B\" : 0.8205 }, \"b\" : 0.7957 },\n",
    "    7 :   { \"t\" : 20.53,  \"a\" : { \"B\" : 0.7305 }, \"b\" : 0.8279 },\n",
    "    8 :   { \"t\" : 29.11,  \"a\" : { \"B\" : 0.6502 }, \"b\" : 0.8553 },\n",
    "    9 :   { \"t\" : 41.20,  \"a\" : { \"B\" : 0.595 },  \"b\" : 0.8757 },\n",
    "    10 :  { \"t\" : 55.19,  \"a\" : { \"B\" : 0.5545 }, \"b\" : 0.8903 },\n",
    "    11 :  { \"t\" : 70.69,  \"a\" : { \"B\" : 0.5333 }, \"b\" : 0.8997 },\n",
    "    12 :  { \"t\" : 90.34,  \"a\" : { \"B\" : 0.5189 }, \"b\" : 0.9073 },\n",
    "    13 :  { \"t\" : 115.29, \"a\" : { \"B\" : 0.5181 }, \"b\" : 0.9122 },\n",
    "    14 :  { \"t\" : 147.42, \"a\" : { \"B\" : 0.5176 }, \"b\" : 0.9171 },\n",
    "    15 :  { \"t\" : 188.24, \"a\" : { \"B\" : 0.5172 }, \"b\" : 0.9217 },\n",
    "    16 :  { \"t\" : 240.03, \"a\" : { \"B\" : 0.5119 }, \"b\" : 0.9267 }\n",
    "}\n",
    "\n",
    "# From Deco for Divers\n",
    "# Do not use: \n",
    "#  need to convert to bar and delta-M to use,\n",
    "#  M0 is in msw and M calculated as Delta M * Depth + M0,\n",
    "#  Depth is unknown\n",
    "#\n",
    "WORKMAN = {\n",
    "    1 : { \"t\" : 5.0, \"M0\" : 31.5, \"M\" : 1.8 },\n",
    "    2 : { \"t\" : 10.0, \"M0\" : 26.8, \"M\" : 1.6 },\n",
    "    3 : { \"t\" : 20.0, \"M0\" : 21.9, \"M\" : 1.5 },\n",
    "    4 : { \"t\" : 40.0, \"M0\" : 17.0, \"M\" : 1.4 },\n",
    "    5 : { \"t\" : 80.0, \"M0\" : 16.4, \"M\" : 1.3 },\n",
    "    6 : { \"t\" : 120.0, \"M0\" : 15.8, \"M\" : 1.2 },\n",
    "    7 : { \"t\" : 160.0, \"M0\" : 15.5, \"M\" : 1.15 },\n",
    "    8 : { \"t\" : 200.0, \"M0\" : 15.5, \"M\" : 1.1 },\n",
    "    9 : { \"t\" : 240.0, \"M0\" : 15.2, \"M\" : 1.1 }\n",
    "}\n",
    "\n",
    "# Also from Deco for Divers\n",
    "#\n",
    "# Since DSAT's primary concern is no-stop diving, it only uses M0\n",
    "#  -- there is no Delta M i.e. Delta M = 1\n",
    "# values are in msw\n",
    "#\n",
    "# Conveert to Buhlmann with\n",
    "#  m_w2b( M0 = NNN / 10, dM = 1, P = 1 )\n",
    "# and run ZHL with DSAT compartments and M-values\n",
    "#\n",
    "DSAT = {\n",
    "    1 : { \"t\" : 5.0, \"M0\" : 30.42 },\n",
    "    2 : { \"t\" : 10.0, \"M0\" : 25.37 },\n",
    "    3 : { \"t\" : 20.0, \"M0\" : 20.54 },\n",
    "    4 : { \"t\" : 30.0, \"M0\" : 18.34 },\n",
    "    5 : { \"t\" : 40.0, \"M0\" : 17.11 },\n",
    "    6 : { \"t\" : 60.0, \"M0\" : 15.79 },\n",
    "    7 : { \"t\" : 80.0, \"M0\" : 15.11 },\n",
    "    8 : { \"t\" : 100.0, \"M0\" : 14.69 },\n",
    "    9 : { \"t\" : 120.0, \"M0\" : 14.41 },\n",
    "    10 : { \"t\" : 160.0, \"M0\" : 14.06 },\n",
    "    11 : { \"t\" : 200.0, \"M0\" : 13.84 },\n",
    "    12 : { \"t\" : 240.0, \"M0\" : 13.69 },\n",
    "    13 : { \"t\" : 360.0, \"M0\" : 13.45 },\n",
    "    14 : { \"t\" : 480.0, \"M0\" : 13.33 }\n",
    "}\n",
    "\n",
    "# return alveolar inert gas pressure\n",
    "# with P amb = 1 bar, fraction of inert gas = 0.79, and RQ = 0.9\n",
    "# this should return 0.79 - 0.0567 = 0.7451 or 0.7452 dep. on where you round it\n",
    "#\n",
    "def palv( Pamb = 1, Q = 0.79, RQ = 0.9 ) :\n",
    "    assert float( RQ ) != 0.0\n",
    "    vw = float( Pamb ) - 0.0627 + (1.0 - float( RQ )) / float( RQ ) * 0.0534\n",
    "    return round( vw * float( Q ), 4 )\n",
    "\n",
    "# return k: constant for tissue compartment (min^-1)\n",
    "# Th : tissue compartment half-time in minutes\n",
    "# for 5-minute compartment it's 0.8452\n",
    "#\n",
    "def kay( Th = 5 ) :\n",
    "    assert float( Th ) > 0.0\n",
    "    return round( math.log( 2 ) / float( Th ), 4 )\n",
    "\n",
    "# return rate of pressure change in bar/min\n",
    "# d0 : start pressure, bar\n",
    "# dt : end pressure, bar\n",
    "# t : time, min\n",
    "# Q : fraction of inert gas (same Q as in palv()\n",
    "#\n",
    "def arr( d0 = 1.0, dt = 1.0, t = 1, Q = 0.79 ) :\n",
    "    assert float( t ) > 0.0\n",
    "    dP = (float( dt ) - float( d0 )) / float( t )\n",
    "    rc = dP * float( Q )\n",
    "    return round( rc, 4 )\n",
    "\n",
    "# Schreiner equation\n",
    "# Palv + R * (t - 1/k) - (Palv - Pi - R/k) * e^(-k * t)\n",
    "#\n",
    "# returns pressure in tissue compartment after time t at depth Pa & dP\n",
    "#\n",
    "# Pi: initial pressure of inert gas in tissue (bar)\n",
    "# Palv: initial pressure of inert gas in the lungs (bar, output of palv())\n",
    "# t: time (minutes)\\n\",\n",
    "# R: rate of pressure change (output of arr()),\n",
    "# k: gas decay constant (output of kay()\n",
    "#\n",
    "# (Intermediate variables b/c I was playing with rounding)\n",
    "#\n",
    "\n",
    "def schreiner( Pi = 0.7451, Palv = 0.7451, t = 1, R = 0, k = 0.1386, verbose = False ) :\n",
    "\n",
    "    assert float( k ) != 0.0\n",
    "    x1 = float( R ) * (float( t ) - 1.0 / float( k ))\n",
    "    x2 = float( Palv ) - float( Pi ) - float( R ) / float( k )\n",
    "    x3 = math.e ** (float( -k ) * float( t ))\n",
    "    rc = round( float( Palv ) + x1 - x2 * x3, 4 )\n",
    "    if verbose : sys.stdout.write( \"x1: %f, x2: %f, x3: %f, rc: %f\\n\" % (x1, x2, x3, rc,) )\n",
    "    return round( rc, 4 )\n",
    "\n",
    "# M-value: workman to buhlmann\n",
    "# P is ambient pressure in bar\n",
    "# returns pair ( a, b )\n",
    "#\n",
    "# TODO: add GF\n",
    "#\n",
    "def m_w2b( M0 = 2.9624, dM = 1.7928, P = 1 ) :\n",
    "    assert float( dM ) >= 1.0\n",
    "    a = float( M0 ) - float( dM ) * float( P )\n",
    "    b = 1.0 / float( dM )\n",
    "    return (round( a, 4 ), round( b, 4 ))\n",
    "\n",
    "# M-value: buhlmann to workman\n",
    "# returns pair ( M0, dM )\n",
    "#\n",
    "def m_b2w( a = 1.1696, b = 0.5578, P = 1 ) :\n",
    "    assert float( b ) > 0.0\n",
    "    M0 = float( a ) + float( P ) / float( b )\n",
    "    dM = 1.0 / float( b )\n",
    "    return (round( M0, 4 ), round( dM, 4 ))\n",
    "\n",
    "# no-stop time by Schreiner\n",
    "#\n",
    "# Palv: initial pressure of inert gas in the lungs (bar, output of palv())\n",
    "# t: time (minutes)\\n\",\n",
    "# R: rate of pressure change (output of arr()),\n",
    "# k: gas decay constant (output of kay()\n",
    "#  -- same as schreiner()\n",
    "# M0: surfacing M-value as per Workman\n",
    "#\n",
    "def ndl( Palv = 0.7451, M0 = 2.9624, t = 0, R = 0, k = 0.1386, verbose = False ) :\n",
    "\n",
    "# (M0 - Palv - R * (t - 1/k)) * math.e ** (k * t) + Palv - R / k\n",
    "    assert float( k ) != 0.0\n",
    "    x1 = float( M0 ) - float( Palv ) - float( R ) * (float( t ) - 1.0 / float( k ))\n",
    "    x2 = math.e ** (float( k ) * float( t ))\n",
    "    rc = x1 * x2 + float( Palv ) - float( R ) / float( k )\n",
    "    if verbose : sys.stdout.write( \"x1: %f, x2: %f, rc: %f\\n\" % (x1, x2, rc,) )\n",
    "    return round( rc, 4 )\n",
    "\n",
    "# Buhlman formula with GF and Helium\n",
    "# returns safe ascent ceiling\n",
    "#\n",
    "# Pn is N2 pressure in tissue compartment\n",
    "# Phe is He pressure in tissue compartment\n",
    "# an is N2 a coefficient\n",
    "# bn is N2 b coefficient\n",
    "# ahe is He a coefficient\n",
    "# bhe is He b coefficient\n",
    "# gf is current gradient factor\n",
    "#\n",
    "def buhlmann( Pn, an, bn, Phe = 0, ahe = 0, bhe = 0, gf = 1 ) :\n",
    "\n",
    "    P = float( Pn ) + float( Phe )\n",
    "    assert float( P ) != 0.0\n",
    "    a = (float( an ) * float( Pn ) + float( ahe ) * float( Phe )) / P\n",
    "    b = (float( bn ) * float( Pn ) + float( bhe ) * float( Phe )) / P\n",
    "    num = float( P ) - float( a ) * float( gf )\n",
    "    den = float( gf ) / float( b ) + 1.0 - float( gf )\n",
    "    assert den != 0.0\n",
    "    rc = num / den\n",
    "    return round( rc, 4 )\n",
    "\n",
    "# eof\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
}