Clean up the handling of surface pressure

There are two ways to look at surface pressure. One is to say "what was
the surface pressure during that dive?" - in that case we now return an
average over the pressure reported by the different divecomputers (or the
standard 1013mbar if none reported any).

Or you want to do specific calculations for a specific divecomputer - in
which case we access only the pressure reported by THAT divecomputer, if
present (and fall back to the previous case, otherwise).

We still have lots of places in Subsurface that only act on the first
divecomputer. As a side effect of this change we now make this more
obvious as we in those cases pass a pointer to the first divecomputer
explicitly to the calculations.

Either way, this commit should prevent us from ever mistakenly basing our
calculations on a surface pressure of 0 (which is the initial bug in
deco.c that triggered all this).

Similar changes need to be made for other elements that we currently only
use from the first divecomputer, i.e., salinity.

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

7 files changed, 51 insertions, 31 deletions

diff --git a/deco.c b/deco.c
index 4bddad576..87899b5ee 100644
--- a/deco.c
+++ b/deco.c
@@ -90,7 +90,7 @@ static double tissue_tolerance_calc(const struct dive *dive)
 	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->dc.surface_pressure.mbar / 1000.0;
+	double surface = get_surface_pressure_in_mbar(dive, TRUE) / 1000.0;
 	double lowest_ceiling;
 
 	for (ci = 0; ci < 16; ci++)
diff --git a/dive.c b/dive.c
index afe115376..7083e857c 100644
--- a/dive.c
+++ b/dive.c
@@ -230,6 +230,22 @@ int get_duration_in_sec(struct dive *dive)
 	return duration;
 }
 
+int get_surface_pressure_in_mbar(const struct dive *dive, gboolean non_null)
+{
+	int count = 0, pressure = 0;
+	const struct divecomputer *dc = &dive->dc;
+	do {
+		if (dc->surface_pressure.mbar) {
+			pressure = (double)(count * pressure + dc->surface_pressure.mbar) / (count + 1) + 0.5;
+			count++;
+		}
+		dc = dc->next;
+	} while (dc);
+	if (!pressure && non_null)
+		pressure = SURFACE_PRESSURE;
+	return pressure;
+}
+
 static void update_temperature(temperature_t *temperature, int new)
 {
 	if (new) {
diff --git a/dive.h b/dive.h
index fcbaca713..ddb9327c6 100644
--- a/dive.h
+++ b/dive.h
@@ -333,17 +333,21 @@ static inline void copy_gps_location(struct dive *from, struct dive *to)
 	}
 }
 
+extern int get_surface_pressure_in_mbar(const struct dive *dive, gboolean non_null);
+
 /* Pa = N/m^2 - so we determine the weight (in N) of the mass of 10m
  * of water (and use standard salt water at 1.03kg per liter if we don't know salinity)
  * and add that to the surface pressure (or to 1013 if that's unknown) */
-static inline int depth_to_mbar(int depth, struct dive *dive)
+static inline int depth_to_mbar(int depth, struct dive *dive, struct divecomputer *dc)
 {
 	double specific_weight = 1.03 * 0.981;
-	int surface_pressure = SURFACE_PRESSURE;
-	if (dive->dc.salinity)
-		specific_weight = dive->dc.salinity / 10000.0 * 0.981;
-	if (dive->dc.surface_pressure.mbar)
-		surface_pressure = dive->dc.surface_pressure.mbar;
+	int surface_pressure;
+	if (dc->salinity)
+		specific_weight = dc->salinity / 10000.0 * 0.981;
+	if (dc->surface_pressure.mbar)
+		surface_pressure = dc->surface_pressure.mbar;
+	else
+		surface_pressure = get_surface_pressure_in_mbar(dive, TRUE);
 	return depth / 10.0 * specific_weight + surface_pressure + 0.5;
 }
 
diff --git a/divelist.c b/divelist.c
index 825c372ed..97981a0f2 100644
--- a/divelist.c
+++ b/divelist.c
@@ -655,7 +655,7 @@ static int calculate_otu(struct dive *dive)
 			po2 = sample->po2;
 		} else {
 			int o2 = active_o2(dive, dc, sample->time);
-			po2 = o2 / 1000.0 * depth_to_mbar(sample->depth.mm, dive);
+			po2 = o2 / 1000.0 * depth_to_mbar(sample->depth.mm, dive, &dive->dc);
 		}
 		if (po2 >= 500)
 			otu += pow((po2 - 500) / 1000.0, 0.83) * t / 30.0;
@@ -720,7 +720,7 @@ static int calculate_sac(struct dive *dive)
 		}
 	}
 	/* Mean pressure in bar (SAC calculations are in bar*l/min) */
-	pressure = depth_to_mbar(dive->dc.meandepth.mm, dive) / 1000.0;
+	pressure = depth_to_mbar(dive->dc.meandepth.mm, dive, &dive->dc) / 1000.0;
 	sac = airuse / pressure * 60 / duration;
 
 	/* milliliters per minute.. */
@@ -744,7 +744,7 @@ static void add_dive_to_deco(struct dive *dive)
 
 		for (j = t0; j < t1; j++) {
 			int depth = interpolate(psample->depth.mm, sample->depth.mm, j - t0, t1 - t0);
-			(void) add_segment(depth_to_mbar(depth, dive) / 1000.0,
+			(void) add_segment(depth_to_mbar(depth, dive, &dive->dc) / 1000.0,
 					   &dive->cylinder[sample->sensor].gasmix, 1, sample->po2, dive);
 		}
 	}
@@ -790,7 +790,7 @@ double init_decompression(struct dive *dive)
 		/* again skip dives from different trips */
 		if (dive->divetrip && dive->divetrip != pdive->divetrip)
 			continue;
-		surface_pressure = (pdive->dc.surface_pressure.mbar ? pdive->dc.surface_pressure.mbar : SURFACE_PRESSURE) / 1000;
+		surface_pressure = get_surface_pressure_in_mbar(pdive, TRUE) / 1000.0;
 		if (!deco_init) {
 			clear_deco(surface_pressure);
 			deco_init = TRUE;
@@ -816,7 +816,7 @@ double init_decompression(struct dive *dive)
 	/* add the final surface time */
 	if (lasttime && dive->when > lasttime) {
 		surface_time = dive->when - lasttime;
-		surface_pressure = (dive->dc.surface_pressure.mbar ? dive->dc.surface_pressure.mbar : SURFACE_PRESSURE) / 1000;
+		surface_pressure = get_surface_pressure_in_mbar(dive, TRUE) / 1000.0;
 		tissue_tolerance = add_segment(surface_pressure, &air, surface_time, 0, dive);
 #if DECO_CALC_DEBUG & 2
 		printf("after surface intervall of %d:%02u\n", FRACTION(surface_time,60));
@@ -824,7 +824,7 @@ double init_decompression(struct dive *dive)
 #endif
 	}
 	if (!deco_init) {
-		double surface_pressure = (dive->dc.surface_pressure.mbar ? dive->dc.surface_pressure.mbar : SURFACE_PRESSURE) / 1000;
+		double surface_pressure = get_surface_pressure_in_mbar(dive, TRUE) / 1000.0;
 		clear_deco(surface_pressure);
 #if DECO_CALC_DEBUG & 2
 		printf("no previous dive\n");
diff --git a/planner.c b/planner.c
index 847128866..5bf87fb1e 100644
--- a/planner.c
+++ b/planner.c
@@ -112,7 +112,7 @@ double tissue_at_end(struct dive *dive, char **cached_datap)
 			lastdepth = psample->depth.mm;
 		for (j = t0; j < t1; j++) {
 			int depth = interpolate(lastdepth, sample->depth.mm, j - t0, t1 - t0);
-			tissue_tolerance = add_segment(depth_to_mbar(depth, dive) / 1000.0,
+			tissue_tolerance = add_segment(depth_to_mbar(depth, dive, &dive->dc) / 1000.0,
 						       &dive->cylinder[gasidx].gasmix, 1, sample->po2, dive);
 		}
 		psample = sample;
@@ -139,7 +139,7 @@ int time_at_last_depth(struct dive *dive, int next_stop, char **cached_data_p)
 	gasidx = get_gasidx(dive, o2, he);
 	while (deco_allowed_depth(tissue_tolerance, surface_pressure, dive, 1) > next_stop) {
 		wait++;
-		tissue_tolerance = add_segment(depth_to_mbar(depth, dive) / 1000.0,
+		tissue_tolerance = add_segment(depth_to_mbar(depth, dive, &dive->dc) / 1000.0,
 					       &dive->cylinder[gasidx].gasmix, 1, sample->po2, dive);
 	}
 	return wait;
@@ -504,7 +504,7 @@ static void add_plan_to_notes(struct diveplan *diveplan, struct dive *dive)
 		gasidx = get_gasidx(dive, o2, he);
 		len = strlen(buffer);
 		if (dp->depth != lastdepth) {
-			used = diveplan->bottomsac / 1000.0 * depth_to_mbar((dp->depth + lastdepth) / 2, dive) *
+			used = diveplan->bottomsac / 1000.0 * depth_to_mbar((dp->depth + lastdepth) / 2, dive, &dive->dc) *
 						(dp->time - lasttime) / 60;
 			snprintf(buffer + len, sizeof(buffer) - len, _("Transition to %.*f %s in %d:%02d min - runtime %d:%02u on %s\n"),
 							decimals, depthvalue, depth_unit,
@@ -514,7 +514,7 @@ static void add_plan_to_notes(struct diveplan *diveplan, struct dive *dive)
 		} else {
 			/* we use deco SAC rate during the calculated deco stops, bottom SAC rate everywhere else */
 			int sac = dp->entered ? diveplan->bottomsac : diveplan->decosac;
-			used = sac / 1000.0 * depth_to_mbar(dp->depth, dive) * (dp->time - lasttime) / 60;
+			used = sac / 1000.0 * depth_to_mbar(dp->depth, dive, &dive->dc) * (dp->time - lasttime) / 60;
 			snprintf(buffer + len, sizeof(buffer) - len, _("Stay at %.*f %s for %d:%02d min - runtime %d:%02u on %s\n"),
 							decimals, depthvalue, depth_unit,
 							FRACTION(dp->time - lasttime, 60),
diff --git a/profile.c b/profile.c
index 86f7537af..3f24966e8 100644
--- a/profile.c
+++ b/profile.c
@@ -1011,16 +1011,16 @@ static void set_sac_color(struct graphics_context *gc, int sac, int avg_sac)
 }
 
 /* calculate the current SAC in ml/min and convert to int */
-#define GET_LOCAL_SAC(_entry1, _entry2, _dive)	(int)				\
+#define GET_LOCAL_SAC(_entry1, _entry2, _dive, _dc)	(int)		\
 	((GET_PRESSURE((_entry1)) - GET_PRESSURE((_entry2))) *			\
 		(_dive)->cylinder[(_entry1)->cylinderindex].type.size.mliter /	\
 		(((_entry2)->sec - (_entry1)->sec) / 60.0) /			\
-		depth_to_mbar(((_entry1)->depth + (_entry2)->depth) / 2.0, (_dive)))
+		depth_to_mbar(((_entry1)->depth + (_entry2)->depth) / 2.0, (_dive), (_dc)))
 
 #define SAC_WINDOW 45	/* sliding window in seconds for current SAC calculation */
 
 static void plot_cylinder_pressure(struct graphics_context *gc, struct plot_info *pi,
-				struct dive *dive)
+				struct dive *dive, struct divecomputer *dc)
 {
 	int i;
 	int last = -1, last_index = -1;
@@ -1050,12 +1050,12 @@ static void plot_cylinder_pressure(struct graphics_context *gc, struct plot_info
 		if (!last_entry) {
 			last = i;
 			last_entry = entry;
-			sac = GET_LOCAL_SAC(entry, pi->entry + i + 1, dive);
+			sac = GET_LOCAL_SAC(entry, pi->entry + i + 1, dive, dc);
 		} else {
 			int j;
 			sac = 0;
 			for (j = last; j < i; j++)
-				sac += GET_LOCAL_SAC(pi->entry + j, pi->entry + j + 1, dive);
+				sac += GET_LOCAL_SAC(pi->entry + j, pi->entry + j + 1, dive, dc);
 			sac /= (i - last);
 			if (entry->sec - last_entry->sec >= SAC_WINDOW) {
 				last++;
@@ -1400,12 +1400,12 @@ static void fill_missing_segment_pressures(pr_track_t *list)
  * scale pressures, so it ends up being a unitless scaling
  * factor.
  */
-static inline int pressure_time(struct dive *dive, struct plot_data *a, struct plot_data *b)
+static inline int pressure_time(struct dive *dive, struct divecomputer *dc, struct plot_data *a, struct plot_data *b)
 {
 	int time = b->sec - a->sec;
 	int depth = (a->depth + b->depth)/2;
 
-	return depth_to_mbar(depth, dive) * time;
+	return depth_to_mbar(depth, dive, dc) * time;
 }
 
 static void fill_missing_tank_pressures(struct dive *dive, struct plot_info *pi, pr_track_t **track_pr)
@@ -1453,7 +1453,7 @@ static void fill_missing_tank_pressures(struct dive *dive, struct plot_info *pi,
 		magic = (segment->end - segment->start) / (double) segment->pressure_time;
 
 		/* Use that overall pressure change to update the current pressure */
-		cur_pt = pressure_time(dive, entry-1, entry);
+		cur_pt = pressure_time(dive, &dive->dc, entry-1, entry);
 		pressure = cur_pr[cyl] + cur_pt * magic + 0.5;
 		INTERPOLATED_PRESSURE(entry) = pressure;
 		cur_pr[cyl] = pressure;
@@ -1768,7 +1768,7 @@ static void populate_pressure_information(struct dive *dive, struct divecomputer
 		struct plot_data *entry = pi->entry + i;
 
 		/* discrete integration of pressure over time to get the SAC rate equivalent */
-		current->pressure_time += pressure_time(dive, entry-1, entry);
+		current->pressure_time += pressure_time(dive, dc, entry-1, entry);
 
 		/* track the segments per cylinder and their pressure/time integral */
 		if (entry->cylinderindex != cylinderindex) {
@@ -1814,7 +1814,7 @@ static void calculate_deco_information(struct dive *dive, struct divecomputer *d
 {
 	int i;
 	double amb_pressure;
-	double surface_pressure = (dive->dc.surface_pressure.mbar ? dive->dc.surface_pressure.mbar : SURFACE_PRESSURE) / 1000.0;
+	double surface_pressure = (dc->surface_pressure.mbar ? dc->surface_pressure.mbar : get_surface_pressure_in_mbar(dive, TRUE)) / 1000.0;
 
 	for (i = 1; i < pi->nr; i++) {
 		int fo2, fhe, j, t0, t1;
@@ -1822,7 +1822,7 @@ static void calculate_deco_information(struct dive *dive, struct divecomputer *d
 		struct plot_data *entry = pi->entry + i;
 		int cylinderindex = entry->cylinderindex;
 
-		amb_pressure = depth_to_mbar(entry->depth, dive) / 1000.0;
+		amb_pressure = depth_to_mbar(entry->depth, dive, dc) / 1000.0;
 		fo2 = dive->cylinder[cylinderindex].gasmix.o2.permille ? : O2_IN_AIR;
 		fhe = dive->cylinder[cylinderindex].gasmix.he.permille;
 		double ratio = (double)fhe / (1000.0 - fo2);
@@ -1876,7 +1876,7 @@ static void calculate_deco_information(struct dive *dive, struct divecomputer *d
 		tissue_tolerance = 0;
 		for (j = t0+1; j <= t1; j++) {
 			int depth = interpolate(entry[-1].depth, entry[0].depth, j - t0, t1 - t0);
-			double min_pressure = add_segment(depth_to_mbar(depth, dive) / 1000.0,
+			double min_pressure = add_segment(depth_to_mbar(depth, dive, dc) / 1000.0,
 							&dive->cylinder[cylinderindex].gasmix, 1, entry->po2 * 1000, dive);
 			tissue_tolerance = min_pressure;
 		}
@@ -2039,7 +2039,7 @@ void plot(struct graphics_context *gc, struct dive *dive, scale_mode_t scale)
 	plot_temperature_profile(gc, pi);
 
 	/* Cylinder pressure plot */
-	plot_cylinder_pressure(gc, pi, dive);
+	plot_cylinder_pressure(gc, pi, dive, dc);
 
 	/* Text on top of all graphs.. */
 	plot_temperature_text(gc, pi);
diff --git a/statistics.c b/statistics.c
index e0df9c405..82070a937 100644
--- a/statistics.c
+++ b/statistics.c
@@ -579,7 +579,7 @@ static void show_single_dive_stats(struct dive *dive)
 	} else {
 		set_label(single_w.air_temp, "");
 	}
-	if (dive->dc.surface_pressure.mbar) {
+	if (get_surface_pressure_in_mbar(dive, FALSE)) {
 		set_label(single_w.air_press, "%d mbar", dive->dc.surface_pressure.mbar);
 	} else {
 		set_label(single_w.air_press, _("unknown"));