\documentclass[a4paper,12pt]{article}
\usepackage[utf8]{inputenc}
\usepackage{mathtools}
\usepackage{amssymb}

\begin{document}
\author{Jan Darowski}
\title{VPM-B algorithm}
\maketitle


\begin{abstract}
This article is a description of the VPM-B decompression algorithm as implemented for the Subsurface project, during 
the Google Summer of Code 2015 program. Its main goal is to explain next phases of the algorithm and an intuitive 
meaning of calculations. It assumes some basic understanding of the Bubble Theory and does not provide full 
derivations of all the equations (however, they can be found in the pointed resources).
\end{abstract}

\section{The Algorithm}
First part of the VPM-B algorithm is a \textit{dive simulation}. It starts from the 
begining of the dive and is intended to gather all the information about a state of the diver's body in the moment of 
the start of ascent. This includes \textit{$He$ and $N_2$ saturation} in tissue compartments and the maximum 
\textit{crushing pressure} for each compartment. These two variables are the main factors, that decide how fast the 
diver can ascend at the end of the dive. 

After the dive simulation, a \textit{start gradient} is calculated. It is the maximum allowable difference between the 
\textit{ambient pressure} and the \textit{tissue saturation}. This gradient is calculated so that the maximum number 
of bubbles in each tissue is not exceeded after the ascent. Then, the \textit{Critical Volume Algorithm} starts. It 
simulates the ascent with the calculated gradient, checks if the gradient can be improved (increased) and starts again 
with a new gradient. It keeps the volume of the generated gas below some constant safe level. Next gradients converge 
pretty fast, so after a few iterations, the final deco profile is set.
For more conservative results, during the ascent simulation, Boyle's law compensation is added. It modifies the 
gradient on every deco stop, to represent the bubbles' expansion during the decompression.

\section{Critical radius}
	Basic VPM assumption is, that in every tissue compartment, some safe number of nuclei can grow into bubbles. 
	Ordering Theory says, that when applying some pressure schedule on a group of nuclei, their size order is kept. So 
	we can introduce some gas-specific critical radius of a nucleon, that is the biggest one that can grow into a 
	bubble. Then, we can simulate the pressure schedule on this nucleon and we will be able to say, what is its final 
	size, just prior the ascent. Knowing that, the only thing we need to do, is prevent it from growing into a 
	bubble. 
	
	For VPM-B, $0.45$ and $0.55$ microns are proposed as the basic values for \textit{$He$} and \textit{$N_2$} nuclei
	respectively. This values can be increased for more conservative deco. 

\section{Dive Simulation}
\subsection{Tissue saturation}
	In the Bubble Theory based algorithms, saturation of \textit{$He$} and \textit{$N_2$} inside 16 different 
	tissue compartment types is calculated. In this paper I will call it $P_{sat}$. This value can be easily calculated 
	using Schreiner's and Haldane's equations. They are quiet popular so I will not describe them here. Haldane 
	equation assumes, that the diver stays at a constant depth. Schreiner equation is its expansion, for a linear 
	ascent / descent.
	
	The original code used both of these equations, but as the Subsurface implementation already had the tissue 
	saturation calculation using only the Haldane's equation on 1s segments, I kept it this way, checking that it 
	does not have any influence on the results.
	
\subsection{Crushing pressure}
	In normal, permeable situation, there is a direct relation between the ambient pressure ($P_{amb}$) and the 
	pressure of gas in a nucleon ($P_{in}$), inside the tissue. In such a situation, we define the crushing pressure 
	as the difference between the ambient pressure and the total tissue saturation:
	\begin{align*}
		H = P_{crush} = P_{amb} - (P_{sat_{He}} + P_{sat_{N_2}} + P_{sat_{other}})
	\end{align*}
	And in this situation, crushing pressure for both gases is the same.
	
	But the nucleon is not always permeable! The $H$ value can be understood as a pressure gradient on the walls of 
	the nucleon. If it's too big, the bubble stops being permeable, gases cannot freely flow in and out the nucleon 
	and $P_{crush} = H$ is no longer kept. VPM-B authors assume some constant value for all the compartments, that 
	defines the critical gradient of permeability $H_{max}$. So what if $H> H_{max}$? Then, for the given nucleon,
	knowing its radius at the beginning of the dive, we can calculate its radius in the moment of becoming 
	impermeable. Next, using the Bubble Theory equations, we can calculate its radius in any moment of the impermeable 
	phase, knowing the ambient pressure. From there, we can calculate the pressure inside the nucleon using the Boyle's 
	law.
	\begin{align*}
		& r_{onset} = \frac{1}{ \frac{H_{max}}{ 2 \cdot (\gamma_c - \gamma ) } + \frac{1}{ r_{start} }} \\
		& A = P_{amb} - H_{max} + \frac{2 \cdot (\gamma_c - \gamma ) }{r_{onset}} \\
		& B = 2 \cdot (\gamma_c - \gamma ) \\
		& C = P_{sat_{onset}} \cdot {r_{onset}}^3 \\
		& \\
		& A \cdot {r_{curr}}^3 - B \cdot {r_{curr}}^2 - C = 0 \\
		& P_{in} = \frac{P_{sat_{onset}} \cdot {r_{onset}}^3}{{r_{curr}}^3} \\
		& P_{crush} = P_{amb} - P_{in} \\
		& \\
	\end{align*}
	\begin{tabular}{l l l}
		& where & \\
		& $r_{onset}$ & nucleon's~radius~in~the~last~permeable~moment \\
		& $r_{start}$ & nucleon's~radius~at~the~beginning~of~the~dive \\
		& $\gamma$ & surface~tension~constant \\
		& $\gamma_c$ & skin~compression~constant \\
		& $P_{amb}$ & current~ambient~pressure \\
		& $P_{sat_{onset}}$ & compartment's~saturation~in~the~last~permeable~moment \\
		& $r_{curr}$ & current~nucleon's~radius
	\end{tabular}
	

\subsection{Bottom time}
	When the profile reaches the end of the bottom time, we can calculate radius of the nuclei starting with the 
	specified \textit{Critical radius}, after applying the crushing pressure on them. Then, authors of the algorithm 
	suggest nuclear regeneration calculations, but their influence on the results is minimal.
	\begin{align*}
		& r = \frac{1}{\frac{P_{crush}}{2 \cdot (\gamma_c - \gamma)} + \frac{1}{r_{start}} } \\
		& r_{regen} = r + ( r_{start} - r) \cdot (1 - e^{\frac{-t}{T}}) \\ 
		& \\
	\end{align*}
	\begin{tabular}{l l l}
		& where  &\\
		& $r$ & nucleon's~radius~after~the~crushing~pressure \\
		& $r_{start}$ & initial~nucleon's~radius \\
		& $r_{regen}$ & radius~after~the~nuclear~regeneration \\
		& $t$ & time~from~the~start~of~the~dive \\
		& $T$ & time~regeneration~constant
	\end{tabular}
	

\section{Start Gradient}
	When we are at the end of the bottom time, we need to find out what conditions would make our tracked critical 
	nucleon grow into a bubble. Having the radius of the nucleon, we can calculate, what gradient between the tissue 
	saturation and ambient pressure would make it grow into bubble ($G_{start}$). Then, we can simulate the ascent, 
	keeping an eye on the current gradient ($G$), stopping when the calculated max gradient would be violated. For the 
	single nucleon, the $G_{start}$ can be formulated as follows:
	
	\begin{align*}
		&G_{max} = G_{start} = 2 \cdot \frac{\gamma}{\gamma_c} \cdot \frac{\gamma_c - \gamma}{r} \\
		& \\
	\end{align*}
	\begin{tabular}{l l l}
		& where & \\
		& $r$ & radius~of~the~nucleon~prior~the~ascent
	\end{tabular}
	\\
	
	We need to calculate it for both \textit{$He$} and \textit{$N_2$} separately. Then, the total gradient for a 
	compartment is given by:
	\begin{align*}
		&G_{start} = \frac{G_{start_{He}} \cdot P_{sat_{He}} + G_{start_{N_2}} \cdot P_{sat_{N_2}} } {P_{sat_{He}} + 
		\cdot P_{sat_{N_2}}}
	\end{align*}
	For the given compartment, $G$ is defined as follows:
	\begin{align*}
		&G = P_{sat_{He}} + P_{sat_{N_2}} + P_{sat_{other}} - P_{amb} \\
		& \\
		\end{align*}
	\begin{tabular}{l l l}
		& where &\\
		& $P_{sat_{other}}$ & other~tissue~gases~constant
	\end{tabular}
	
	
\section{Ascent Simulation}
	During the ascent simulation, we just need to keep track of the current $G$ value and make a deco stop when it gets 
	too close to $G_{max}$. That's the theory. Original code by Baker makes something slightly different. At the end of 
	each deco stop, it waits, until it could reach the next stop in $0$ time, without violating $G_{max}$. It makes the 
	schedule a little bit more conservative than if we take the ascent time into account. Especially for dives with 
	a short bottom time. In Subsurface, the original idea has been implemented.
	
\subsection{Boyle's Law Compensation}
	The letter B in VPM-B states for the Boyle's Law Compensation, which is a small modification of the original 
	algorithm, that greatly increases deco time of deeper dives. The idea is, that because of a depth difference 
	during the ascent, nuclei grow. So the gradient needed to turn the initial nucleon into bubble decreases. To 
	implement this idea, we assume, that the $G_{start}$ was the gradient working fine only to the first deco stop. 
	Then, on every stop, we update the current $G_{max}$, so it works fine for the next deco stop.\\
	What is important from the practical point of view: also gas change stops should be treated as deco stops. If the 
	first gas change is before any other deco stop, all the compensation should be done in relation to it.
	\begin{align*}
		& r_{first} =  \frac{2 \cdot \gamma }{G_{bottom}}\\
		& A = P_{amb_{next} } \\
		& B = -2 \cdot \gamma \\
		& C = (P_{amb_{first} } + \frac{2 \cdot \gamma }{r_{first} } ) \cdot {r_{first} }^3 \\
		& A \cdot {r_{next}}^3 - B \cdot {r_{next}}^2 - C = 0
		& \\
		& G_{max_{new}} = \frac{2 \cdot \gamma}{r_{next}}
		& \\
	\end{align*}
	\begin{tabular}{l l l}
		& where &\\
		& $r_{first}$ & radius~of~nucleon~at~the~first~stop \\
		& $G_{bottom}$ & gradient~for~the~first~stop~($G_{start}$~or~from~CVA) \\
		& $P_{amb_{next}}$ & ambient~pressure~at~the~next~stop \\
		& $P_{amb_{first}}$ & ambient~pressure~at~the~first~stop
	\end{tabular}
	

\section{Critical Volume Algorithm}
	The basic VPM algorithm, without CVA, assumes that the creation of bubbles from nuclei is instant and body can 
	tolerate some number ($N_{safe}$) of them forever. But in real dives, it is possible to accept some more bubbles, 
	as long as they do not exceed the critical volume ($V_{crit}$) of free gas generated during the desaturation. The 
	free gas volume is assumed to be proportional to the integral of $G$ in time, during the desaturation. To simplify 
	calculations, We assume, that through the whole ascent $G = G_{max}$. And so we can easily calculate $G \cdot 
	t_{deco}$. But the free gas is also generated at the surface, where the gradient decays exponentially.

	CVA's goal is to generate a new $G_{max_{new}}$, that will be closer to generating $V_{crit}$ free gas. It is not 
	precise and depends on the deco time, so we must simulate the whole ascent several times, recalculating $G_{max}$, 
	recording $t_{deco}$  and checking if it still differs significantly from the previous time. The formula for the 
	next $G_{max_{new}}$ is as follows:
	\begin{align*}
		& B = G_{start} + \frac{\lambda \cdot \gamma}{ \gamma_c \cdot t_{desat}} \\
		& C = \frac{ \gamma^2 \cdot \lambda \cdot P_{crush} }{{\gamma_c}^2 \cdot t_{desat} } \\
		& G_{max} = \frac{B + \sqrt{B^2 - 4C}}{2}\\
		& \\
	\end{align*}
	\begin{tabular}{l l l}
		& where & \\
		& $\lambda$ & critical~volume $\times$ time~constant
	\end{tabular}
	\\
	
	Of course it needs to be calculated for all 16 compartments, for two gases and the final result is a weighted 
	average of $G_{max_{He}}$ and $G_{max_{N_2}}$. Just like in the case of $G_{start}$. What can be unclear here is 
	the $t_{desat}$. It is the total time of desaturation, assuming that it is done with $G = G_{max_{old}}$, also on 
	the surface!
	\begin{align*}
		& t_{surf1} = \frac{ \frac{P_{sat_{He}}}{ T_{He}} + \frac{P_{sat_{N_{2}}} - P_{N_2 Ins}}{ 
		T_{N_{2}}}} {P_{sat_{He}} + P_{sat_{N_{2}}} - P_{N_2 Ins}} \\
		& \\
		& T_G = \frac{ln(\frac{P_{N_2 Ins} - P_{sat_{N_2}}} {P_{sat_{He}}} )}    { T_{N_2} - T_{He} } \\
		& VT_G = \frac{P_{sat_{He}}}{T_{He}} \cdot (1 - \exp(-T_{He} \cdot T_G)) +
				\frac{P_{sat_{N_2}} - P_{N_2 Ins}}{T_{N_2}} \cdot (1 - \exp(-T_{N_2} \cdot T_G)) \\
		& t_{surf2} = \frac{VT_G} {P_{sat_{He}} +P_{sat_{N_2}} - P_{N_2 Ins}} \\
		& \\
		& t_{desat} = t_{deco} + t_{surf} \\
		& \\
	\end{align*}
	\begin{tabular}{l l l}
		& where &\\
		& $t_{surf1}$ & when~$P_{sat_{N_2}} >  P_{N_2 Ins}$ \\
		& $t_{surf2}$ & otherwise \\
		& $P_{N_2 Ins}$ & inspired~$N_2$~pressure~at~surface \\
		& $T_x$ & half\-life~of~x 
	\end{tabular}
	
	
\section{Conservatism levels}
	As the basic parameterization is prepared for very aggressive divers, additional conservatism levels are 
	implemented. Next levels scale the starting critical radii by $1.05, 1.12, 1.22, 1.35$. These values were proposed 
	by the V-Planner's authors and provide enough spread for most users.

\section{Parameters}
	\begin{tabular}{l l l}
		$r_{start}$ &$N_2$~critical~radius &- 0.55~microns \\
		$ r_{start}$ &$He$~critical~radius &- 0.45~microns \\
		$ \lambda$  &Critical~volume $\times$ time~ &- 199.58~bar $\cdot$ min \\
		$ H_{max}$ &Max~gradient~of~impermeability~ &- 8.2~bar \\
		$ \gamma$ &Surface~tension~ &- 0.0179~N/m \\
		$ \gamma_c$ &Skin~compression~ &- 0.257~N/m\\
		$ T$ &Nucleon~regeneration~time~ &- 20160~min\\
		$ P_{sat_{other}}$ &Other~gases~saturation~ &- 0.1359888~bar\\
	\end{tabular}


\section{Resources}
\begin{enumerate}
	\item \textbf{http://www.deepocean.net} - Very good explanation of the Bubble Theory and the breathing equations. 
	It is a good place to start, not necessarily to base implementation on, as there are a few mistakes.
	\item \textbf{https://github.com/bwaite/vpmb} - Rewritten original implementation of VPM-B in C and Python. Very 
	handful for debugging and comparing partial results. 
	\item \textbf{http://www.decompression.org/maiken/VPM/VPM\_Program\_Site\_Map.htm} - Site of the authors of VPM-B. 
	Links to the ftp server containing priceless materials, like CVA equation's derivation or the original code.
	
\end{enumerate}
\end{document}
