Show effective gradient factors for VPMB-plans

For each stop, this computes an effective gradient factor
that gives the same ceiling. Then, it does linear regression
to find values for GFlow and GFhigh that give a similar deco
profile.

Note that this optimises the average gradient factor. The
runtime however depends strongly at the gradient factor at
the last depth. So we don't necessarily to get the runtime
right.

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

3 files changed, 78 insertions, 5 deletions

diff --git a/core/deco.c b/core/deco.c
index b138a9fbd..57acfe951 100644
--- a/core/deco.c
+++ b/core/deco.c
@@ -30,6 +30,7 @@
 
 
 extern bool in_planner();
+extern int plot_depth;
 
 extern pressure_t first_ceiling_pressure;
 
@@ -175,6 +176,10 @@ double bottom_he_gradient[16];
 double initial_n2_gradient[16];
 double initial_he_gradient[16];
 
+int sumx, sum1;
+long sumxx;
+double sumy, sumxy;
+
 double get_crit_radius_He()
 {
 	if (vpmb_config.conservatism <= 4)
@@ -308,6 +313,23 @@ double tissue_tolerance_calc(const struct dive *dive, double pressure)
 			}
 		// We are doing ok if the gradient was computed within ten centimeters of the ceiling.
 		} while (fabs(ret_tolerance_limit_ambient_pressure - reference_pressure) > 0.01);
+
+		if (plot_depth) {
+			++sum1;
+			sumx += plot_depth;
+			sumxx += plot_depth * plot_depth;
+			double n2_gradient, he_gradient, total_gradient;
+			n2_gradient = update_gradient(depth_to_bar(plot_depth, &displayed_dive), bottom_n2_gradient[ci_pointing_to_guiding_tissue]);
+			he_gradient = update_gradient(depth_to_bar(plot_depth, &displayed_dive), bottom_he_gradient[ci_pointing_to_guiding_tissue]);
+			total_gradient = ((n2_gradient * tissue_n2_sat[ci_pointing_to_guiding_tissue]) + (he_gradient * tissue_he_sat[ci_pointing_to_guiding_tissue]))
+					/ (tissue_n2_sat[ci_pointing_to_guiding_tissue] + tissue_he_sat[ci_pointing_to_guiding_tissue]);
+
+			double buehlmann_gradient = (1.0 / buehlmann_inertgas_b[ci_pointing_to_guiding_tissue] - 1.0) * depth_to_bar(plot_depth, &displayed_dive) + buehlmann_inertgas_a[ci_pointing_to_guiding_tissue];
+			double gf = (total_gradient - vpmb_config.other_gases_pressure) / buehlmann_gradient;
+			sumxy += gf * plot_depth;
+			sumy += gf;
+			plot_depth = 0;
+		}
 	}
 	return ret_tolerance_limit_ambient_pressure;
 }
@@ -632,3 +654,31 @@ double get_gf(double ambpressure_bar, const struct dive *dive)
 		gf = gf_low;
 	return gf;
 }
+
+double regressiona()
+{
+	if (sum1) {
+		double avxy = sumxy / sum1;
+		double avx = (double)sumx / sum1;
+		double avy = sumy / sum1;
+		double avxx = (double) sumxx / sum1;
+		return (avxy - avx * avy) / (avxx - avx*avx);
+	}
+	else
+		return 0.0;
+}
+
+double regressionb()
+{
+	if (sum1)
+		return sumy / sum1 - sumx * regressiona() / sum1;
+	else
+		return 0.0;
+}
+
+void reset_regression()
+{
+	sumx = sum1 = 0;
+	sumxx = 0L;
+	sumy = sumxy = 0.0;
+}
diff --git a/core/dive.h b/core/dive.h
index f1c144481..4a314dad5 100644
--- a/core/dive.h
+++ b/core/dive.h
@@ -859,6 +859,7 @@ struct diveplan {
 	short gfhigh;
 	short vpmb_conservatism;
 	struct divedatapoint *dp;
+	int eff_gflow, eff_gfhigh;
 };
 
 struct divedatapoint *plan_add_segment(struct diveplan *diveplan, int duration, int depth, int cylinderid, int po2, bool entered);
diff --git a/core/planner.c b/core/planner.c
index 3e67b4674..488fcb9c4 100644
--- a/core/planner.c
+++ b/core/planner.c
@@ -34,10 +34,14 @@ int decostoplevels_imperial[] = { 0, 3048, 6096, 9144, 12192, 15240, 18288, 2133
 double plangflow, plangfhigh;
 bool plan_verbatim, plan_display_runtime, plan_display_duration, plan_display_transitions;
 
+extern double regressiona();
+extern double regressionb();
+extern void reset_regression();
+
 pressure_t first_ceiling_pressure, max_bottom_ceiling_pressure = {};
 
 const char *disclaimer;
-
+int plot_depth = 0;
 #if DEBUG_PLAN
 void dump_plan(struct diveplan *diveplan)
 {
@@ -577,10 +581,15 @@ static void add_plan_to_notes(struct diveplan *diveplan, struct dive *dive, bool
 		snprintf(temp, sz_temp, translate("gettextFromC", "based on Bühlmann ZHL-16C with GFlow = %d and GFhigh = %d"),
 			diveplan->gflow, diveplan->gfhigh);
 	} else if (prefs.deco_mode == VPMB){
+		int temp_len;
 		if (diveplan->vpmb_conservatism == 0)
-			snprintf(temp, sz_temp, "%s", translate("gettextFromC", "based on VPM-B at nominal conservatism"));
+			temp_len = snprintf(temp, sz_temp, "%s", translate("gettextFromC", "based on VPM-B at nominal conservatism"));
 		else
-			snprintf(temp, sz_temp, translate("gettextFromC", "based on VPM-B at +%d conservatism"), diveplan->vpmb_conservatism);
+			temp_len = snprintf(temp, sz_temp, translate("gettextFromC", "based on VPM-B at +%d conservatism"), diveplan->vpmb_conservatism);
+		if(diveplan->eff_gflow)
+			temp_len += snprintf(temp + temp_len, sz_temp - temp_len,  translate("gettextFromC", ", effective GF=%d/%d"), diveplan->eff_gflow
+					     , diveplan->eff_gfhigh);
+
 	} else if (prefs.deco_mode == RECREATIONAL){
 		snprintf(temp, sz_temp, translate("gettextFromC", "recreational mode based on Bühlmann ZHL-16B with GFlow = %d and GFhigh = %d"),
 			diveplan->gflow, diveplan->gfhigh);
@@ -996,6 +1005,7 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 	int breakcylinder = 0;
 	int error = 0;
 	bool decodive = false;
+	int first_stop_depth = 0;
 
 	set_gf(diveplan->gflow, diveplan->gfhigh, prefs.gf_low_at_maxdepth);
 	set_vpmb_conservatism(diveplan->vpmb_conservatism);
@@ -1165,6 +1175,7 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 		gas = bottom_gas;
 		stopping = false;
 		decodive = false;
+		first_stop_depth = 0;
 		stopidx = bottom_stopidx;
 		breaktime = -1;
 		breakcylinder = 0;
@@ -1183,7 +1194,7 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 			report_error(translate("gettextFromC", "Can't find gas %s"), gasname(&gas));
 			current_cylinder = 0;
 		}
-
+		reset_regression();
 		while (1) {
 			/* We will break out when we hit the surface */
 			do {
@@ -1204,6 +1215,9 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 								TIMESTEP, po2, &displayed_dive, prefs.decosac);
 				clock += TIMESTEP;
 				depth -= deltad;
+				/* Print VPM-Gradient as gradient factor, this has to be done from within deco.c */
+				if (decodive)
+					plot_depth = depth;
 			} while (depth > 0 && depth > stoplevels[stopidx]);
 
 			if (depth <= 0)
@@ -1258,7 +1272,10 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 					break; /* We did not hit the ceiling */
 
 				/* Add a minute of deco time and then try again */
-				decodive = true;
+				if (!decodive) {
+					decodive = true;
+					first_stop_depth = depth;
+				}
 				if (!stopping) {
 					/* The last segment was an ascend segment.
 					 * Add a waypoint for start of this deco stop */
@@ -1349,6 +1366,11 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 	add_plan_to_notes(diveplan, &displayed_dive, show_disclaimer, error);
 	fixup_dc_duration(&displayed_dive.dc);
 
+	if(prefs.deco_mode == VPMB) {
+		diveplan->eff_gfhigh = rint(100.0 * regressionb());
+		diveplan->eff_gflow = rint(100*(regressiona() * first_stop_depth + regressionb()));
+	}
+
 	free(stoplevels);
 	free(gaschanges);
 	free(bottom_cache);