From 404205c8f1aa37be374c2ba17acceb3541620bce Mon Sep 17 00:00:00 2001 From: "Robert C. Helling" Date: Wed, 23 Nov 2016 11:50:50 +0100 Subject: [PATCH 2/2] Show effective gradient factors for VPMB-plans To: subsurface@subsurface-divelog.org 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 --- core/deco.c | 50 ++++++++++++++++++++++++++++++++++++++++++++++++++ core/dive.h | 1 + core/planner.c | 32 +++++++++++++++++++++++++++----- 3 files changed, 78 insertions(+), 5 deletions(-) diff --git a/core/deco.c b/core/deco.c index b138a9f..57acfe9 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 f1c1444..4a314da 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 3e67b46..488fcb9 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); -- 2.8.4 (Apple Git-73)