Deco artefacts with low GFlow

In a dive, when you choose a very low GFlow (like 5 or 9) and a trimix
with quite some He (12/48 in the example) and descend fast, the ceiling
seems to do strange things in the first minutes of the dive (very very
deep for example or jumping around).

To understand what is going on we have to recall what gradient factors do
in detail: Plain Buehlmann gives you for each tissue a maximal inert gas
pressure that is a straight line when plotted against the ambient
pressure. So for each depth (=ambient pressure) there is a maximally
allowed over-pressure.

The idea of gradient factors is that one does not use all the possible
over-pressure that Buehlmann gives us but only a depth dependent fraction.
GFhigh is the fraction of the possible over-pressure at the surface while
GFlow is the fraction at the first deco stop. In between, the fraction is
linearly interpolated. As the Buehlmann over-pressure is increasing with
depth and typically also the allowed overpressure after applications of
gradient factors increases with depth or said differently: the tissue
saturation has to be lower if the diver wants to ascent.

The main problem is: What is the first stop (where to apply GFlow)? In a
planned dive, we could take the first deco stop, but in a real dive from a
dive computer download it is impossible to say what constitutes a stop and
what is only a slow ascent?

What I have used so far is not exactly the first stop but rather the first
theoretical stop: During all of the dive, I have calculated the ceiling
under the assumption that GFlow applies everywhere (and not just at a
single depth). The deepest of these ceilings I have used as the “first
stop depth”, the depth at which GFlow applies.

Even more, I only wanted to use the information that a diver has during
the dive, so I actually only considered the ceilings in the past (and not
in the future of a given sample).

But this brings with it the problem that early in the dive, in particular
during the descent the lowest ceiling so far is very shallow (as not much
gas has built up in the body so far).

This problem now interferes with a second one: If at the start of the dive
when the all compartments have 790mbar N2 the diver starts breathing a
He-heavy mix (like 12/48) and descents fast the He builds up in the
tissues before the N2 can diffuse out. So right at the start, we already
encounter high tissue loadings.

If now we have a large difference between GFhigh and GFlow but they apply
at very similar depth (the surface and a very shallow depth of the deepest
ceiling (which for a non-decompression dive would be theoretically at
negative depth) so far) it can happen that the linear interpolation as
opposite slope then in the typical case above: The allowed over-pressure
is degreasing with depth, shallower depth do not require lower gas loading
in the tissue (i.e. can be reached after further off-gasing) but but
tolerate higher loadings. In that situation the ceiling disappears (or is
rather a floor).

So far, I got rid of that problem, by stating that the minimum depth for
GFlow was 20m (after all, GFlow is about deep stops, so it should better
not be too shallow). Now the dive reported in ticket #549 takes values to
an extreme in such away that 20m (which is determined by
buehlmann_config.gf_low_position_min in deco.c) was not enough to prevent
this inversion problem (or in a milder form that the interpolation of
gradient factors is in fact an extrapolation with quite extreme values).

This patch that gets rid of the problem for the dive described above but
still it is possible to find (more extreme) parameter choices that lead to
non-realistic ceilings.

Let me close by pointing out that all this is only about the descent, as
it is about too shallow depth for GFlow. So no real deco (i.e. later part
of the dive) is inflicted. This is only about a theoretical ceiling
displayed possibly in the first minutes of a dive. So this is more an
aesthetically than a practical problem.

Fixes #549

Signed-off-by: Robert C. Helling <helling@atdotde.de>
Signed-off-by: Dirk Hohndel <dirk@hohndel.org>

1 files changed, 29 insertions, 16 deletions

diff --git a/deco.c b/deco.c
index f5ffb8958..ffa8ed4ba 100644
--- a/deco.c
+++ b/deco.c
@@ -29,7 +29,7 @@ struct buehlmann_config {
 	double gf_low_position_min;	 //! gf_low_position below surface_min_shallow.
 	bool gf_low_at_maxdepth;	    //! if true, gf_low applies at max depth instead of at deepest ceiling.
 };
-struct buehlmann_config buehlmann_config = { 1.0, 1.01, 0, 0.75, 0.35, 2.0, false };
+struct buehlmann_config buehlmann_config = { 1.0, 1.01, 0, 0.75, 0.35, 1.0, false };
 
 const double buehlmann_N2_a[] = { 1.1696, 1.0, 0.8618, 0.7562,
 				  0.62, 0.5043, 0.441, 0.4,
@@ -90,38 +90,49 @@ double tolerated_by_tissue[16];
 static double tissue_tolerance_calc(const struct dive *dive)
 {
 	int ci = -1;
-	double tissue_inertgas_saturation, buehlmann_inertgas_a, buehlmann_inertgas_b;
+	double tissue_inertgas_saturation[16], buehlmann_inertgas_a[16], buehlmann_inertgas_b[16];
 	double ret_tolerance_limit_ambient_pressure = 0.0;
 	double gf_high = buehlmann_config.gf_high;
 	double gf_low = buehlmann_config.gf_low;
 	double surface = get_surface_pressure_in_mbar(dive, true) / 1000.0;
+	double lowest_ceiling = 0.0;
+	double tissue_lowest_ceiling[16];
 
 	for (ci = 0; ci < 16; ci++) {
-		double tolerated;
+		tissue_inertgas_saturation[ci] = tissue_n2_sat[ci] + tissue_he_sat[ci];
+		buehlmann_inertgas_a[ci] = ((buehlmann_N2_a[ci] * tissue_n2_sat[ci]) + (buehlmann_He_a[ci] * tissue_he_sat[ci])) / tissue_inertgas_saturation[ci];
+		buehlmann_inertgas_b[ci] = ((buehlmann_N2_b[ci] * tissue_n2_sat[ci]) + (buehlmann_He_b[ci] * tissue_he_sat[ci])) / tissue_inertgas_saturation[ci];
 
-		tissue_inertgas_saturation = tissue_n2_sat[ci] + tissue_he_sat[ci];
-		buehlmann_inertgas_a = ((buehlmann_N2_a[ci] * tissue_n2_sat[ci]) + (buehlmann_He_a[ci] * tissue_he_sat[ci])) / tissue_inertgas_saturation;
-		buehlmann_inertgas_b = ((buehlmann_N2_b[ci] * tissue_n2_sat[ci]) + (buehlmann_He_b[ci] * tissue_he_sat[ci])) / tissue_inertgas_saturation;
 
 		/* tolerated = (tissue_inertgas_saturation - buehlmann_inertgas_a) * buehlmann_inertgas_b; */
 
+		tissue_lowest_ceiling[ci] = (buehlmann_inertgas_b[ci] * tissue_inertgas_saturation[ci] - gf_low * buehlmann_inertgas_a[ci] * buehlmann_inertgas_b[ci]) /
+					     ((1.0 - buehlmann_inertgas_b[ci]) * gf_low + buehlmann_inertgas_b[ci]);
+		if (tissue_lowest_ceiling[ci] > lowest_ceiling)
+			lowest_ceiling = tissue_lowest_ceiling[ci];
 		if (!buehlmann_config.gf_low_at_maxdepth) {
-			double lowest_ceiling = (buehlmann_inertgas_b * tissue_inertgas_saturation - gf_low * buehlmann_inertgas_a * buehlmann_inertgas_b) /
-						((1.0 - buehlmann_inertgas_b) * gf_low + buehlmann_inertgas_b);
 			if (lowest_ceiling > gf_low_pressure_this_dive)
 				gf_low_pressure_this_dive = lowest_ceiling;
 		}
+	}
+	for (ci = 0; ci <16; ci++) {
+		double tolerated;
+
+		if ((surface / buehlmann_inertgas_b[ci] + buehlmann_inertgas_a[ci] - surface) * gf_high + surface <
+		    (gf_low_pressure_this_dive / buehlmann_inertgas_b[ci] + buehlmann_inertgas_a[ci] - gf_low_pressure_this_dive) * gf_low + gf_low_pressure_this_dive)
+			tolerated = (-buehlmann_inertgas_a[ci] * buehlmann_inertgas_b[ci] * (gf_high * gf_low_pressure_this_dive - gf_low * surface) -
+				     (1.0 - buehlmann_inertgas_b[ci]) * (gf_high - gf_low) * gf_low_pressure_this_dive * surface +
+				     buehlmann_inertgas_b[ci] * (gf_low_pressure_this_dive - surface) * tissue_inertgas_saturation[ci]) /
+				    (-buehlmann_inertgas_a[ci] * buehlmann_inertgas_b[ci] * (gf_high - gf_low) +
+				     (1.0 - buehlmann_inertgas_b[ci]) * (gf_low * gf_low_pressure_this_dive - gf_high * surface) +
+				     buehlmann_inertgas_b[ci] * (gf_low_pressure_this_dive - surface));
+		else
+			tolerated = ret_tolerance_limit_ambient_pressure;
 
-		tolerated = (-buehlmann_inertgas_a * buehlmann_inertgas_b * (gf_high * gf_low_pressure_this_dive - gf_low * surface) -
-			     (1.0 - buehlmann_inertgas_b) * (gf_high - gf_low) * gf_low_pressure_this_dive * surface +
-			     buehlmann_inertgas_b * (gf_low_pressure_this_dive - surface) * tissue_inertgas_saturation) /
-			    (-buehlmann_inertgas_a * buehlmann_inertgas_b * (gf_high - gf_low) +
-			     (1.0 - buehlmann_inertgas_b) * (gf_low * gf_low_pressure_this_dive - gf_high * surface) +
-			     buehlmann_inertgas_b * (gf_low_pressure_this_dive - surface));
 
 		tolerated_by_tissue[ci] = tolerated;
 
-		if (tolerated > ret_tolerance_limit_ambient_pressure) {
+		if (tolerated >= ret_tolerance_limit_ambient_pressure) {
 			ci_pointing_to_guiding_tissue = ci;
 			ret_tolerance_limit_ambient_pressure = tolerated;
 		}
@@ -231,7 +242,9 @@ void clear_deco(double surface_pressure)
 		tissue_n2_sat[ci] = (surface_pressure - WV_PRESSURE) * N2_IN_AIR / 1000;
 		tissue_he_sat[ci] = 0.0;
 	}
-	gf_low_pressure_this_dive = surface_pressure + buehlmann_config.gf_low_position_min;
+	gf_low_pressure_this_dive = surface_pressure;
+	if (!buehlmann_config.gf_low_at_maxdepth)
+		gf_low_pressure_this_dive += buehlmann_config.gf_low_position_min;
 }
 
 void cache_deco_state(double tissue_tolerance, char **cached_datap)