Use boyle_compensation in profile

otherwise VPM-B planned profiles seem to violate the ceiling. This needs
the first_stop_pressure to be available also in the profile, so I made
it global in planner.c

Important lesson: If you want to use deco_allowed_depth on a tissue_tolerance
that comes from a VPM-B planned dive, you have to call boyles_law() before
add_segment()!

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

4 files changed, 17 insertions, 9 deletions

diff --git a/deco.c b/deco.c
index 21332e1e6..7b8d7c2be 100644
--- a/deco.c
+++ b/deco.c
@@ -23,6 +23,8 @@
 
 extern bool in_planner();
 
+extern int first_stop_pressure;
+
 //! Option structure for Buehlmann decompression.
 struct buehlmann_config {
 	double satmult;			    //! safety at inert gas accumulation as percentage of effect (more than 100).
@@ -344,24 +346,27 @@ double solve_cubic(double A, double B, double C)
 
 }
 
-double update_gradient(double first_stop_pressure, double next_stop_pressure, double first_gradient)
+double update_gradient(double next_stop_pressure, double first_gradient)
 {
 	double first_radius = 2.0 * vpmb_config.surface_tension_gamma / first_gradient;
 	double A = next_stop_pressure;
 	double B = -2.0 * vpmb_config.surface_tension_gamma;
-	double C = (first_stop_pressure + 2.0 * vpmb_config.surface_tension_gamma / first_radius) * cube(first_radius);
+	double C = (first_stop_pressure / 1000.0 + 2.0 * vpmb_config.surface_tension_gamma / first_radius) * cube(first_radius);
 
 	double next_radius = solve_cubic(A, B, C);
 
 	return 2.0 * vpmb_config.surface_tension_gamma / next_radius;
 }
 
-void boyles_law(double first_stop_pressure, double next_stop_pressure)
+void boyles_law(double next_stop_pressure)
 {
 	int ci;
+
+	if (!first_stop_pressure)
+		return;
 	for (ci = 0; ci < 16; ++ci) {
-		allowable_n2_gradient[ci] = update_gradient(first_stop_pressure, next_stop_pressure, bottom_n2_gradient[ci]);
-		allowable_he_gradient[ci] = update_gradient(first_stop_pressure, next_stop_pressure, bottom_he_gradient[ci]);
+		allowable_n2_gradient[ci] = update_gradient(next_stop_pressure, bottom_n2_gradient[ci]);
+		allowable_he_gradient[ci] = update_gradient(next_stop_pressure, bottom_he_gradient[ci]);
 
 		total_gradient[ci] = ((allowable_n2_gradient[ci] * tissue_n2_sat[ci]) + (allowable_he_gradient[ci] * tissue_he_sat[ci])) / (tissue_n2_sat[ci] + tissue_he_sat[ci]);
 	}
diff --git a/dive.h b/dive.h
index ee3f63701..0125b76a5 100644
--- a/dive.h
+++ b/dive.h
@@ -801,7 +801,7 @@ extern double restore_deco_state(char *data);
 extern void nuclear_regeneration(double time);
 extern void vpmb_start_gradient();
 extern void vpmb_next_gradient(double deco_time, double surface_pressure);
-extern void boyles_law(double first_stop_pressure, double next_stop_pressure);
+extern void boyles_law(double next_stop_pressure);
 
 /* this should be converted to use our types */
 struct divedatapoint {
diff --git a/planner.c b/planner.c
index f7b3a3f93..0c3d3dc4a 100644
--- a/planner.c
+++ b/planner.c
@@ -33,6 +33,8 @@ 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;
 
+int first_stop_pressure;
+
 const char *disclaimer;
 
 #if DEBUG_PLAN
@@ -921,7 +923,6 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 	int bottom_depth;
 	int bottom_gi;
 	int bottom_stopidx;
-	int first_stop_pressure;
 	bool is_final_plan = true;
 	int deco_time;
 	int previous_deco_time;
@@ -996,6 +997,7 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 		return(false);
 	}
 	tissue_tolerance = tissue_at_end(&displayed_dive, cached_datap);
+	displayed_dive.surface_pressure.mbar = diveplan->surface_pressure;
 
 #if DEBUG_PLAN & 4
 	printf("gas %s\n", gasname(&gas));
@@ -1160,7 +1162,7 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 				// Boyles Law compensation
 				if (first_stop_pressure == 0)
 					first_stop_pressure = depth_to_mbar(depth, &displayed_dive);
-				boyles_law(first_stop_pressure / 1000.0, depth_to_mbar(stoplevels[stopidx], &displayed_dive) / 1000.0);
+				boyles_law(depth_to_mbar(stoplevels[stopidx], &displayed_dive) / 1000.0);
 
 				/* Check we need to change cylinder.
 				 * We might not if the cylinder was chosen by the user
@@ -1215,7 +1217,7 @@ bool plan(struct diveplan *diveplan, char **cached_datap, bool is_planner, bool
 					// Boyles Law compensation
 					if (first_stop_pressure == 0)
 						first_stop_pressure = depth_to_mbar(depth, &displayed_dive);
-					boyles_law(first_stop_pressure / 1000.0, depth_to_mbar(stoplevels[stopidx], &displayed_dive) / 1000.0);
+					boyles_law(depth_to_mbar(stoplevels[stopidx], &displayed_dive) / 1000.0);
 				}
 
 				/* Are we waiting to switch gas?
diff --git a/profile.c b/profile.c
index 51f322cfd..67168547b 100644
--- a/profile.c
+++ b/profile.c
@@ -853,6 +853,7 @@ void calculate_deco_information(struct dive *dive, struct divecomputer *dc, stru
 			time_stepsize = t1 - t0;
 		for (j = t0 + time_stepsize; j <= t1; j += time_stepsize) {
 			int depth = interpolate(entry[-1].depth, entry[0].depth, j - t0, t1 - t0);
+			boyles_law(depth_to_mbar(entry->depth, &displayed_dive) / 1000.0);
 			double min_pressure = add_segment(depth_to_mbar(depth, dive) / 1000.0,
 							  &dive->cylinder[entry->cylinderindex].gasmix, time_stepsize, entry->o2pressure.mbar, dive, entry->sac);
 			tissue_tolerance = min_pressure;