Graph documentation [was: Re: Meaning of GF settings]

[Robert C. Helling]

On Thu, 10 Jan 2013, Jan Schubert wrote:

> Robert, I'm not arguing that your calculation is wrong, it's just a
> little bit hard to understand. I also do not have access to Mathematica,
> does it make sense to post a screenshot from the solving part?
> I'd like to check the outcome of the deco calculation in total (not
> specific to the function above) is correct, especially as we see such
> notable differences to dive computers.

OK, let me be a bit more explicit.

We use (and I did not change anything there as it is good as far as I can 
see) tissue loadings according to Buehlmann as calculated in 
add_segment(). You can view this as a black box giving 
tissue_inert_gas_saturation for each compartment (actually it gives this 
separately for N and He and they are combined in the beginning of the loop 
in tissue_tolerance_calc()). For some background, there are many 
explanations on the web. Of course, I would recommend my own (in German)

http://euve10195.vserver.de/~robert/dekotheorie.pdf 
("Dekotheorie verstehen")

>From this saturation, you want to calculate the ceiling, or equivalently 
the minimal ambient pressure. Buehlmann claims there is an affine relation 
between the two (i.e. the graph of the function is a straight line):

min_amb_pressure = b * (tissue_saturation - a)

for some constants that Buehlmann gives separately for each compartment. 
Again, actually, he gives those separately for N and He, and you are 
supposed to average them weighted by the amounts of the respective gases 
in the compartement. This is done (and again I did not tough this) in the 
loop in tissue_tolerance_calc() where buehlmann_inertgas_a and 
buehlmann_inertgas_b are computed. This is the classical Buehlmann model.

To understand gradient factors, according to Eric Baker, see

ftp://gordonhamilton.com/Diving/Technical/deepstop.pdf

it's best to first invert the question: What is the maximal tissue 
saturation at a given ambient pressure? According to Buehlmann (see 
above), this is

tissue_saturation = ambient_pressure / b + a .

For historical reasons, this tissue saturation is called "M-value". At 
this point you should have a look at Fig. 4 of the Baker paper. This plot 
has the ambient pressure as the x-axis and the tissue_saturation as 
y-axis. You can see that the M-value line is always above the line for the 
ambient pressure since a and b are positive constants.

The idea of gradent factors is now not to use the full tolerated 
overpressure

M_value - ambient pressure

but only a fraction GF of it, i.e.

M_value_corrected = GF * (M_value - ambient_pressure) + ambient_pressure

                   = GF * (ambient_pressure / b + a - ambient_pressure) + ambient_pressure .

Let's call this equation (1). Let's further make GF depth dependent, such 
that in the deeper part (the early part of the ascent) it is much smaller 
than close to the surface.

Again, it is assumed that GF as a function of depth is such that, when 
M_value_corrected plotted as a function of ambient pressure, it is a 
streight line, the dashed line in Baker's Fig. 4. In formulas, this means

M_value_corrected = alpha * ambient_pressure + beta

for some constants alpha and beta. Call this equation (2). Now, let's 
determine these constants.

First, at the surface, the gradient factor should be GF_high, i.e. alpha 
and beta have to be such that (equating (1) and (2))

alpha * surface_pressure + beta = GF_high (surface_p / b + a - surface_p)
 	                                   + surface_p

Similarly, at the maximal depth (*) the gradient factor should be GF_low, 
i.e. it should hold

alpha * max_amb_p + beta = GF_low (max_amb_p / b + a - max_amb_p)
                                            + max_amb_p

We have now two points that the gradient factor line should go through. 
This determines the line, that is we can solve these last two equations 
for alpha and beta. This you could do by hand but I used Mathematica, see 
the PDF that I attached of that session. The result for alpha and beta is 
Out[95].

Finally, recall that in (2) the M_value_corrected is not really the 
unknown, this is our tissue_saturation, but we want to determine the 
corresponding ambient pressure. So we invert (2) to

ambient_pressure = (M_value_corrected - beta) / alpha

where we plug in the values of alpha and beta that we just found. This is 
Out[96] in the Mathematica file or the six-line-formula in 
tissue_tolerance_calc.

Of course, that can be broken into smaller parts, like first computing 
temporary variables for nominator and denominator. But I don't see that 
that would be any clearer since at least I cannot see an indepentent 
meaning of the nominator or denominator separately.

> Is there a way to split the tissue calculation from them GF
> calculation/interpolation?

I don't see any since the tissue calculation involves the GF's. Note, that 
by directly computing the inverted (2), we never have to compute an 
interpolated GF at any depth.

> Do you agree when I state GFlow should not have any (notable) impact on
> the shallow part/deco on existing profiles?

Not directly. Strictly speaking it only has no influence at the surface. 
In Baker's Fig. 4, you can see this as how much the dashed line moves when 
you move the "GF low generates first stop" point left and right. It moves 
the line everywhere, of course the less the further up you come.

> Can we check our understanding of GFlow? I remember you saying GFlow
> refers to deepest point in dive!? I thought GFlow defines the 1st deco
> stop (current GF should be <=GFlow).

OK, you are right, I had remembered this wrong but this is what Baker 
states. Note that this moves the point where GF_low applies up which 
implies that at all shallower depth (i.e. during deco) a lower GF results 
which makes the deco longer rather than shorter.

Of course, "GF_low" applies at first deco stop is a but tricky since the 
depth of the first deco stop again depends on GF_low, i.e. there is 
another equation to solve. You can do this by inverting equation (1) for 
the ambient pressure and use GF_low as the gradient factor. This yields

amb = (b * M_value_corrected - GF_low * a * b) / ((1-b) * GF_low + b)

I also attached a patch that implements this and you can switch between 
the to versions by setting the macro GF_LOW_AT_MAXDEPTH to true or false 
for you to play around with this. (It has the effect of changing the 
meaning of gf_low_pressure_this_dive from max_ambient_pressure to 
max_ceiling_pressure).

I found that for me this had little effect in the deco part of my 
reference dive. I also hat to change buehlmann_config.gf_low_position_min 
for it to have any noticible effect but for deeper dives this does not 
matter.

> How is the interpolation from GFlow to GFhigh to surface? In my
> understanding GFhigh refers to the max. GF when surfacing.

See above.

Let me add a final remark: Deco is an art, not a science. Any algorithm 
will only give you a rough guide how to run your decompression and 
hopefully be safe most of the time while not wasting too much time with 
unneccessary stops. But it is always a very rough approximation of very 
complex physiological processes in your body that even in principle cannot 
be perfectly simulated since we do not now the initial conditions of many 
if the relevant variables (like the exact composition of your blood and 
the shape of your blood vesses etc).

Thus I would not worry at all if two different models (or implementations) 
vary in their prediction for your stops by one or even two meters. That 
said, I find the precision at which the smooth ceiling matches your actual 
profile and what your dive computer says very very impressive (say from a 
depth of 50m). I expected much worse. Acutally, I would call that spot on, 
in particular given that your actual depth during deco stages also has 
fluctuations of the size of at least one meter. Which, as was mentioned in 
an earlier email is also the order of magnitude of surface pressure 
variations due to weather. And those you also do not take into account.

The mouse over value ("1m") is actually rounded to integer multiples of 
meters. If you want more precise values, find the word "Calculated" in 
profile.c and in the snprintf replace depthvalue by entry->ceiling. Then 
the ceiling is reported in mm and you will see that at the end of your 
dive it is even less than a meter I would expect (just don't commit that 
change).

Was this the explanation you were looking for?

Best
Robert

PS: For the patch

Signed-off-by: Robert C. Helling helling at atdotde.de

-- 
.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oOo.oO
Robert C. Helling     Elite Master Course Theoretical and Mathematical Physics
                       Scientific Coordinator
 		      Ludwig Maximilians Universitaet Muenchen, Dept. Physik
print "Just another   Phone: +49 89 2180-4523  Theresienstr. 39, rm. B339
     stupid .sig\n";   http://www.atdotde.de
-------------- next part --------------
A non-text attachment was scrubbed...
Name: deco.pdf
Type: application/pdf
Size: 57187 bytes
Desc: 
URL: <http://lists.hohndel.org/pipermail/subsurface/attachments/20130110/39f43f75/attachment-0001.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: gf_low_at_ceiling.patch
Type: text/x-diff
Size: 2594 bytes
Desc: 
URL: <http://lists.hohndel.org/pipermail/subsurface/attachments/20130110/39f43f75/attachment-0001.patch>