Meaning of GF settings

Robert C. Helling helling at lmu.de
Tue Jan 8 06:37:41 PST 2013


Hi,


On Sun, 6 Jan 2013, Dirk Hohndel wrote:

>> Having your code in place will incredibly helpful since with this I
>> don't have to think about how to integrate this with the rest of the
>> program. If you have some patience with me for a few days, I will try
>> to come up with my version of tissue_loadings_to_ceiling.
>
> I will stop messing with the deco routines and wait for your rewrite.

ok, here we go. This does look somewhat reasonable to me. Please check 
that it gives sensible results for your test cases!

Here is what I did:

o) Instead of using gradient factors as means of comparison, I now use 
pressure (as in: maximal ambient pressure).

o) tissue_tolerance_calc() now computes the maximal ambient pressure now 
respecting gradient factors. For this, it needs to know about the surface 
pressure (as refernce for GF_high), thus gets *dive as an argument. It is 
called from add_segment() which this also needs *dive as an additional 
argument.

o) This implies deco_allowed_depth is now mainly a ambient-pressure to 
depth conversion with decorations to avoid negative depth (i.e. no deco 
obliation), implementation of quantization (!smooth => multiples of 3m) 
and explicit setting of last deco depth (e.g. 6m for O2 deco).

o) gf_low_pressure_this_dive (slight change of name), the max depth in 
pressure units is updated in add_segment. I set the minimal value in 
buehlmann_config to the equivalent of 20m as otherwise good values of 
GF_low add a lot of deco to shallow dives which do not need deep stops in 
the first place.

o) The bogus loop is gone as well as actual_gradient_limit() and 
gradient_factor_calculation() and large parts of deco_allowed_depth() 
although I did not delete the code but put it in comments.

o) The meat is in the formula in lines 147-154 of deco.c. Here is the 
rationale:

Without gradient factors, the M-value (i.e the maximal tissue pressure) 
at a given depth is given by ambient_pressure / buehlmann_b + a.

According to "Clearing Up The Confusion About "Deep Stops" by Erik C. 
Baker (as found via google) the effect of the gradient factors is no 
replace this by a reduced affine relation (i.e. another line) such that at 
the surface the difference between M-value and ambient pressure is reduced 
by a factor GF_high and at the maximal depth by a factor GF_low.

That is, we are looking for parameters alpha and beta such that

alpha surface + beta = surface + gf_high * (surface/b + a - surface)

and

alpha max_p + beta = max_p + gf_low * (max_p/b + a - max_p)

This can be solved for alpha and beta (I did this in Mathematica, notebook 
attached as well) and then inverted to obtain the max ambient pressure 
given tissue loadings. The result is the above mentioned formula.

Have fun!
Robert

PS: And now I can start my paid work for today...

-- 
.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 --------------
diff --git a/deco.c b/deco.c
index 841ca31..3197bd8 100644
--- a/deco.c
+++ b/deco.c
@@ -31,7 +31,7 @@ struct buehlmann_config {
   double  gf_high_emergency;	//! emergency gf factors
   double  gf_low_emergency;	//! gradient factor low (at bottom/start of deco calculation).
 };
-struct buehlmann_config buehlmann_config = { 1.0, 1.01, 0.5, 3, 0.75, 0.35, 1.0, 6.0, 0.95, 0.95 };
+struct buehlmann_config buehlmann_config = { 1.0, 1.01, 0.5, 0, 0.75, 0.35, 3.0, 6.0, 0.95, 0.95 };
 struct dive_data {
         double pressure;	//! pesent ambient pressure
         double surface;		//! pressure at water surface
@@ -86,50 +86,55 @@ const double buehlmann_He_factor_expositon_one_second[] = {
 #define PRESSURE_CHANGE_3M 0.3
 #define TOLERANCE 0.02
 
+#define DECO_STOPS_MULTIPLIER_MM 3000.0
+
 double tissue_n2_sat[16];
 double tissue_he_sat[16];
 double tissue_tolerated_ambient_pressure[16];
 int ci_pointing_to_guiding_tissue;
-double gf_low_position_this_dive;
+double gf_low_pressure_this_dive;
 #define TISSUE_ARRAY_SZ sizeof(tissue_n2_sat)
 
-static double actual_gradient_limit(const struct dive_data *data)
-{
-	double pressure_diff, limit_at_position;
-	double gf_high = buehlmann_config.gf_high;
-	double gf_low = buehlmann_config.gf_low;
-
-	pressure_diff = data->pressure - data->surface;
-
-	if (pressure_diff > TOLERANCE) {
-		if (pressure_diff < gf_low_position_this_dive)
-			limit_at_position = gf_high - ((gf_high - gf_low) * pressure_diff / gf_low_position_this_dive);
-		else
-			limit_at_position = gf_low;
-	} else {
-		limit_at_position = gf_high;
-	}
-	return limit_at_position;
-}
-
-static double gradient_factor_calculation(const struct dive_data *data)
-{
-	double tissue_inertgas_saturation;
-
-	tissue_inertgas_saturation = tissue_n2_sat[ci_pointing_to_guiding_tissue] +
-					tissue_he_sat[ci_pointing_to_guiding_tissue];
-	if (tissue_inertgas_saturation < data->pressure)
-		return 0.0;
-	else
-		return (tissue_inertgas_saturation - data->pressure) /
-			(tissue_inertgas_saturation - tissue_tolerated_ambient_pressure[ci_pointing_to_guiding_tissue]);
-}
-
-static double tissue_tolerance_calc(void)
+/* static double actual_gradient_limit_obsolete(const struct dive_data *data) */
+/* { */
+/* 	double pressure_diff, limit_at_position; */
+/* 	double gf_high = buehlmann_config.gf_high; */
+/* 	double gf_low = buehlmann_config.gf_low; */
+
+/* 	pressure_diff = data->pressure - data->surface; */
+
+/* 	if (pressure_diff > TOLERANCE) { */
+/* 		if (pressure_diff < gf_low_position_this_dive) */
+/* 			limit_at_position = gf_high - ((gf_high - gf_low) * pressure_diff / gf_low_position_this_dive); */
+/* 		else */
+/* 			limit_at_position = gf_low; */
+/* 	} else { */
+/* 		limit_at_position = gf_high; */
+/* 	} */
+/* 	return limit_at_position; */
+/* } */
+
+/* static double gradient_factor_calculation_obsolete(const struct dive_data *data) */
+/* { */
+/* 	double tissue_inertgas_saturation; */
+
+/* 	tissue_inertgas_saturation = tissue_n2_sat[ci_pointing_to_guiding_tissue] + */
+/* 					tissue_he_sat[ci_pointing_to_guiding_tissue]; */
+/* 	if (tissue_inertgas_saturation < data->pressure) */
+/* 		return 0.0; */
+/* 	else */
+/* 		return (tissue_inertgas_saturation - data->pressure) / */
+/* 			(tissue_inertgas_saturation - tissue_tolerated_ambient_pressure[ci_pointing_to_guiding_tissue]); */
+/* } */
+
+static double tissue_tolerance_calc(const struct dive *dive)
 {
 	int ci = -1;
 	double tissue_inertgas_saturation, buehlmann_inertgas_a, buehlmann_inertgas_b;
 	double ret_tolerance_limit_ambient_pressure = 0.0;
+	double gf_high = buehlmann_config.gf_high;
+	double gf_low = buehlmann_config.gf_low;
+	double surface = dive->surface_pressure.mbar / 1000.0;
 
 	for (ci = 0; ci < 16; ci++)
 	{
@@ -137,8 +142,16 @@ static double tissue_tolerance_calc(void)
 		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;
 
-		tissue_tolerated_ambient_pressure[ci] = (tissue_inertgas_saturation - buehlmann_inertgas_a) * buehlmann_inertgas_b;
+		/* tissue_tolerated_ambient_pressure[ci] = (tissue_inertgas_saturation - buehlmann_inertgas_a) * buehlmann_inertgas_b; */
 
+		tissue_tolerated_ambient_pressure[ci] = (-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));
+		
 		if (tissue_tolerated_ambient_pressure[ci] > ret_tolerance_limit_ambient_pressure)
 		{
 			ci_pointing_to_guiding_tissue = ci;
@@ -149,13 +162,16 @@ static double tissue_tolerance_calc(void)
 }
 
 /* add a second at the given pressure and gas to the deco calculation */
-double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds, double ccpo2)
+double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds, double ccpo2, const struct dive *dive)
 {
 	int ci;
 	int fo2 = gasmix->o2.permille ? gasmix->o2.permille : 209;
 	double ppn2 = (pressure - WV_PRESSURE) * (1000 - fo2 - gasmix->he.permille) / 1000.0;
 	double pphe = (pressure - WV_PRESSURE) * gasmix->he.permille / 1000.0;
 
+	if(pressure > gf_low_pressure_this_dive)
+	        gf_low_pressure_this_dive = pressure;
+
 	if (ccpo2 > 0.0) { /* CC */
 		double rel_o2_amb, f_dilutent;
 		rel_o2_amb = ccpo2 / pressure;
@@ -200,7 +216,7 @@ double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds
 					(1 - pow(2.0,(- period_in_seconds / (buehlmann_He_t_halflife[ci] * 60))));
 		}
 	}
-	return tissue_tolerance_calc();
+	return tissue_tolerance_calc(dive);
 }
 
 void dump_tissues()
@@ -223,7 +239,7 @@ void clear_deco(double surface_pressure)
 		tissue_he_sat[ci] = 0.0;
 		tissue_tolerated_ambient_pressure[ci] = 0.0;
 	}
-	gf_low_position_this_dive = buehlmann_config.gf_low_position_min;
+	gf_low_pressure_this_dive = surface_pressure+ buehlmann_config.gf_low_position_min;
 }
 
 void cache_deco_state(double tissue_tolerance, char **cached_datap)
@@ -240,7 +256,7 @@ void cache_deco_state(double tissue_tolerance, char **cached_datap)
 	data += TISSUE_ARRAY_SZ;
 	memcpy(data, tissue_tolerated_ambient_pressure, TISSUE_ARRAY_SZ);
 	data += TISSUE_ARRAY_SZ;
-	memcpy(data, &gf_low_position_this_dive, sizeof(double));
+	memcpy(data, &gf_low_pressure_this_dive, sizeof(double));
 	data += sizeof(double);
 	memcpy(data, &tissue_tolerance, sizeof(double));
 	data += sizeof(double);
@@ -257,7 +273,7 @@ double restore_deco_state(char *data)
 	data += TISSUE_ARRAY_SZ;
 	memcpy(tissue_tolerated_ambient_pressure, data, TISSUE_ARRAY_SZ);
 	data += TISSUE_ARRAY_SZ;
-	memcpy(&gf_low_position_this_dive, data, sizeof(double));
+	memcpy(&gf_low_pressure_this_dive, data, sizeof(double));
 	data += sizeof(double);
 	memcpy(&tissue_tolerance, data, sizeof(double));
 	data += sizeof(double);
@@ -268,47 +284,62 @@ double restore_deco_state(char *data)
 
 unsigned int deco_allowed_depth(double tissues_tolerance, double surface_pressure, struct dive *dive, gboolean smooth)
 {
-	unsigned int depth, multiples_of_3m;
-	gboolean below_gradient_limit;
-	double new_gradient_factor;
-	double pressure_delta = tissues_tolerance - surface_pressure;
-	struct dive_data mydata;
-	int bail = 1000;
-
-	if (pressure_delta > 0) {
-		if (!smooth) {
-			multiples_of_3m = (pressure_delta + DIST_FROM_3_MTR) / 0.3;
-			depth = 3000 * multiples_of_3m;
-		} else {
-			depth = rel_mbar_to_depth(pressure_delta * 1000, dive);
-		}
-	} else {
-		depth = 0;
-	}
-	mydata.pressure = depth_to_mbar(depth, dive) / 1000.0;
-	mydata.surface = surface_pressure;
-
-	new_gradient_factor = gradient_factor_calculation(&mydata);
-	below_gradient_limit = (new_gradient_factor < actual_gradient_limit(&mydata));
-	while(!below_gradient_limit)
-	{
-		/* we run into bugs where this turns into an infinite loop; so add
-		 * some bailout code that prints a warning but prevents the code from hanging */
-		if (--bail == 0) {
-			printf("WARNING!!!\n==========\nThe deco_allowed_depth() loop appears to hang.\nBailing out.\n");
-			break;
-		}
-		if (!smooth)
-			mydata.pressure += PRESSURE_CHANGE_3M;
-		else
-			mydata.pressure += PRESSURE_CHANGE_3M / 30; /* 4in / 10cm instead */
-		new_gradient_factor = gradient_factor_calculation(&mydata);
-		below_gradient_limit = (new_gradient_factor < actual_gradient_limit(&mydata));
-	}
-	depth = rel_mbar_to_depth((mydata.pressure - surface_pressure) * 1000, dive);
+	unsigned int depth;
+	double pressure_delta;
+	
+	pressure_delta = tissues_tolerance > surface_pressure ? tissues_tolerance - surface_pressure : 0.0;  /* Avoid negative depths */
+	
+	
+	/* gboolean below_gradient_limit; */
+	/* double new_gradient_factor; */
+	/* double pressure_delta = tissues_tolerance - surface_pressure; */
+	/* struct dive_data mydata; */
+	/* int bail = 1000; */
+
+	/* if (pressure_delta > 0) { */
+	/* 	if (!smooth) { */
+	/* 		multiples_of_3m = (pressure_delta + DIST_FROM_3_MTR) / 0.3; */
+	/* 		depth = 3000 * multiples_of_3m; */
+	/* 	} else { */
+	/* 		depth = rel_mbar_to_depth(pressure_delta * 1000, dive); */
+	/* 	} */
+	/* } else { */
+	/* 	depth = 0; */
+	/* } */
+	/* mydata.pressure = depth_to_mbar(depth, dive) / 1000.0; */
+	/* mydata.surface = surface_pressure; */
+
+	/* new_gradient_factor = gradient_factor_calculation(&mydata); */
+	/* below_gradient_limit = (new_gradient_factor < actual_gradient_limit(&mydata)); */
+	/* while(!below_gradient_limit) */
+	/* { */
+	/* 	/\* we run into bugs where this turns into an infinite loop; so add */
+	/* 	 * some bailout code that prints a warning but prevents the code from hanging *\/ */
+	/* 	if (--bail == 0) { */
+	/* 		printf("WARNING!!!\n==========\nThe deco_allowed_depth() loop appears to hang.\nBailing out.\n"); */
+	/* 		break; */
+	/* 	} */
+	/* 	if (!smooth) */
+	/* 		mydata.pressure += PRESSURE_CHANGE_3M; */
+	/* 	else */
+	/* 		mydata.pressure += PRESSURE_CHANGE_3M / 30; /\* 4in / 10cm instead *\/ */
+	/* 	new_gradient_factor = gradient_factor_calculation(&mydata); */
+	/* 	below_gradient_limit = (new_gradient_factor < actual_gradient_limit(&mydata)); */
+	/* } */
+	/* depth = rel_mbar_to_depth((mydata.pressure - surface_pressure) * 1000, dive); */
+
+	depth = rel_mbar_to_depth(pressure_delta * 1000, dive);
+	
+	if(!smooth)
+	  depth = ceil(depth / DECO_STOPS_MULTIPLIER_MM) * DECO_STOPS_MULTIPLIER_MM;
+
+	if(depth > 0 && depth < buehlmann_config.last_deco_stop_in_mtr * 1000)
+	  depth = buehlmann_config.last_deco_stop_in_mtr * 1000;
+	
 	return depth;
 }
 
+
 void set_gf(double gflow, double gfhigh)
 {
 	if (gflow != -1.0)
diff --git a/dive.h b/dive.h
index 45f0709..ec03bf5 100644
--- a/dive.h
+++ b/dive.h
@@ -572,7 +572,7 @@ extern void subsurface_command_line_exit(gint *, gchar ***);
 
 #define FRACTION(n,x) ((unsigned)(n)/(x)),((unsigned)(n)%(x))
 
-extern double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds, double setpoint);
+extern double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds, double setpoint, const struct dive *dive);
 extern void clear_deco(double surface_pressure);
 extern void dump_tissues(void);
 extern unsigned int deco_allowed_depth(double tissues_tolerance, double surface_pressure, struct dive *dive, gboolean smooth);
diff --git a/divelist.c b/divelist.c
index 6b09b7d..aaf7d2f 100644
--- a/divelist.c
+++ b/divelist.c
@@ -837,7 +837,7 @@ static void add_dive_to_deco(struct dive *dive)
 			int depth = 0.5 + psample->depth.mm + (j - t0) *
 					(sample->depth.mm - psample->depth.mm) / (t1 - t0);
 			(void) add_segment(depth_to_mbar(depth, dive) / 1000.0,
-					&dive->cylinder[sample->sensor].gasmix, 1, sample->po2 / 1000.0);
+					   &dive->cylinder[sample->sensor].gasmix, 1, sample->po2 / 1000.0, dive);
 		}
 	}
 }
@@ -882,7 +882,7 @@ double init_decompression(struct dive *dive)
 		printf("added dive #%d\n", pdive->number);
 		dump_tissues();
 #endif
-		tissue_tolerance = add_segment(surface_pressure, &air, surface_time, 0.0);
+		tissue_tolerance = add_segment(surface_pressure, &air, surface_time, 0.0, dive);
 #if DECO_CALC_DEBUG & 2
 		printf("after surface intervall of %d:%02u\n", FRACTION(surface_time,60));
 		dump_tissues();
diff --git a/planner.c b/planner.c
index 683fbef..6c9015c 100644
--- a/planner.c
+++ b/planner.c
@@ -65,7 +65,7 @@ double tissue_at_end(struct dive *dive, char **cached_datap)
 		for (j = t0; j < t1; j++) {
 			int depth = psample->depth.mm + (j - t0) * (sample->depth.mm - psample->depth.mm) / (t1 - t0);
 			tissue_tolerance = add_segment(depth_to_mbar(depth, dive) / 1000.0,
-						&dive->cylinder[sample->sensor].gasmix, 1, sample->po2);
+						       &dive->cylinder[sample->sensor].gasmix, 1, sample->po2, dive);
 		}
 		psample = sample;
 		t0 = t1;
@@ -90,7 +90,7 @@ int time_at_last_depth(struct dive *dive, int next_stop, char **cached_data_p)
 	while (deco_allowed_depth(tissue_tolerance, surface_pressure, dive, 1) > next_stop) {
 		wait++;
 		tissue_tolerance = add_segment(depth_to_mbar(depth, dive) / 1000.0,
-					&dive->cylinder[sample->sensor].gasmix, 1, sample->po2);
+					       &dive->cylinder[sample->sensor].gasmix, 1, sample->po2, dive);
 	}
 	return wait;
 }
diff --git a/profile.c b/profile.c
index c1603c3..1448b17 100644
--- a/profile.c
+++ b/profile.c
@@ -1870,7 +1870,7 @@ static void calculate_deco_information(struct dive *dive, struct divecomputer *d
 			for (j = t0; j < t1; j++) {
 				int depth = 0.5 + (entry - 1)->depth + (j - t0) * (entry->depth - (entry - 1)->depth) / (t1 - t0);
 				double min_pressure = add_segment(depth_to_mbar(depth, dive) / 1000.0,
-								&dive->cylinder[cylinderindex].gasmix, 1, entry->po2);
+								  &dive->cylinder[cylinderindex].gasmix, 1, entry->po2, dive);
 				if (min_pressure > tissue_tolerance)
 					tissue_tolerance = min_pressure;
 			}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: deco.nb
Type: application/mathematica
Size: 5467 bytes
Desc: 
URL: <http://lists.hohndel.org/pipermail/subsurface/attachments/20130108/3c97e9c5/attachment-0001.nb>


More information about the subsurface mailing list