Consider previous dives when calculating deco

This also initializes the N2 tissue saturations to correct numbers
(setting them to zero was clearly silly).

With this commit we walk back in the dive_table until we find a surface
intervall that's longer than 48h. Or a dive that comes after the last one
we looked at; that would indicate that this is a divelist that contains
dives from multiple divers or dives that for other reasons are not
ordered. In a sane environment one would assume that the dives that need
to be taken into account when doing deco calculations are organized as one
trip in the XML file and so this logic should work.

One major downside of the current implementation is that we recalculate
everything whenever the plot_info is recreated - which happens quite
frequently, for example when resizing the window or even when we go into
loup mode. While this isn't all that compute intensive, this is an utter
waste and we should at least cache the saturation inherited from previous
dives (and clear that number when the selected dive changes). We don't
want to cache all of it as the recreation of the plot_info may be
triggered by the user changing equipment (and most importantly, gasmix)
information. In that case the deco data for this dive does indeed have to
be recreated. But without changing the current dive the saturation after
the last surface intervall should stay the same.

Signed-off-by: Dirk Hohndel <dirk@hohndel.org>

5 files changed, 137 insertions, 18 deletions

diff --git a/deco.c b/deco.c
index 998d6637c..e18779ac8 100644
--- a/deco.c
+++ b/deco.c
@@ -7,10 +7,11 @@
  * The implementation below is (C) Dirk Hohndel 2012 and released under the GPLv2
  *
  * clear_deco()  - call to initialize for a new deco calculation
- * add_segment(pressure, gasmix) - add 1 second at the given pressure, breathing gasmix
+ * add_segment(pressure, gasmix, seconds) - add <seconds> at the given pressure, breathing gasmix
  * deco_allowed_depth(tissues_tolerance, surface_pressure, dive, smooth)
  *				- ceiling based on lead tissue, surface pressure, 3m increments or smooth
  */
+#include <math.h>
 #include "dive.h"
 
 //! Option structure for Buehlmann decompression.
@@ -76,6 +77,7 @@ const double buehlmann_He_factor_expositon_one_second[] = {
 	1.00198406028040E-004, 7.83611475491108E-005, 6.13689891868496E-005, 4.81280465299827E-005};
 
 #define WV_PRESSURE 0.0627 /* water vapor pressure */
+#define N2_IN_AIR 0.7902
 #define DIST_FROM_3_MTR 0.28
 #define PRESSURE_CHANGE_3M 0.3
 #define TOLERANCE 0.02
@@ -145,32 +147,65 @@ 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)
+double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds)
 {
 	int ci;
 	double ppn2 = (pressure - WV_PRESSURE) * (1000 - gasmix->o2.permille - gasmix->he.permille) / 1000.0;
 	double pphe = (pressure - WV_PRESSURE) * gasmix->he.permille / 1000.0;
 
-	divetime++;
 	/* right now we just do OC */
-	for (ci = 0; ci < 16; ci++) {
-		if (ppn2 - tissue_n2_sat[ci] > 0)
-			tissue_n2_sat[ci] += buehlmann_config.satmult * (ppn2 - tissue_n2_sat[ci]) * buehlmann_N2_factor_expositon_one_second[ci];
-		else
-			tissue_n2_sat[ci] += buehlmann_config.desatmult * (ppn2 - tissue_n2_sat[ci]) * buehlmann_N2_factor_expositon_one_second[ci];
-		if (pphe - tissue_he_sat[ci] > 0)
-			tissue_he_sat[ci] += buehlmann_config.satmult * (pphe - tissue_he_sat[ci]) * buehlmann_He_factor_expositon_one_second[ci];
-		else
-			tissue_he_sat[ci] += buehlmann_config.desatmult * (pphe - tissue_he_sat[ci]) * buehlmann_He_factor_expositon_one_second[ci];
+	if (period_in_seconds == 1) { /* that's what we do during the dive */
+		for (ci = 0; ci < 16; ci++) {
+			if (ppn2 - tissue_n2_sat[ci] > 0)
+				tissue_n2_sat[ci] += buehlmann_config.satmult * (ppn2 - tissue_n2_sat[ci]) *
+								buehlmann_N2_factor_expositon_one_second[ci];
+			else
+				tissue_n2_sat[ci] += buehlmann_config.desatmult * (ppn2 - tissue_n2_sat[ci]) *
+								buehlmann_N2_factor_expositon_one_second[ci];
+			if (pphe - tissue_he_sat[ci] > 0)
+				tissue_he_sat[ci] += buehlmann_config.satmult * (pphe - tissue_he_sat[ci]) *
+								buehlmann_He_factor_expositon_one_second[ci];
+			else
+				tissue_he_sat[ci] += buehlmann_config.desatmult * (pphe - tissue_he_sat[ci]) *
+								buehlmann_He_factor_expositon_one_second[ci];
+		}
+	} else { /* all other durations */
+		for (ci = 0; ci < 16; ci++)
+		{
+			if (ppn2 - tissue_n2_sat[ci] > 0)
+				tissue_n2_sat[ci] += buehlmann_config.satmult * (ppn2 - tissue_n2_sat[ci]) *
+					(1 - pow(2.0,(- period_in_seconds / (buehlmann_N2_t_halflife[ci] * 60))));
+			else
+				tissue_n2_sat[ci] += buehlmann_config.desatmult * (ppn2 - tissue_n2_sat[ci]) *
+					(1 - pow(2.0,(- period_in_seconds / (buehlmann_N2_t_halflife[ci] * 60))));
+			if (pphe - tissue_he_sat[ci] > 0)
+				tissue_he_sat[ci] += buehlmann_config.satmult * (pphe - tissue_he_sat[ci]) *
+					(1 - pow(2.0,(- period_in_seconds / (buehlmann_He_t_halflife[ci] * 60))));
+			else
+				tissue_he_sat[ci] += buehlmann_config.desatmult * (pphe - tissue_he_sat[ci]) *
+					(1 - pow(2.0,(- period_in_seconds / (buehlmann_He_t_halflife[ci] * 60))));
+		}
 	}
 	return tissue_tolerance_calc();
 }
 
-void clear_deco()
+void dump_tissues()
+{
+	int ci;
+	printf("N2 tissues:");
+	for (ci = 0; ci < 16; ci++)
+		printf(" %6.3e", tissue_n2_sat[ci]);
+	printf("\nHe tissues:");
+	for (ci = 0; ci < 16; ci++)
+		printf(" %6.3e", tissue_he_sat[ci]);
+	printf("\n");
+}
+
+void clear_deco(double surface_pressure)
 {
 	int ci;
 	for (ci = 0; ci < 16; ci++) {
-		tissue_n2_sat[ci] = 0.0;
+		tissue_n2_sat[ci] = (surface_pressure - WV_PRESSURE) * N2_IN_AIR;
 		tissue_he_sat[ci] = 0.0;
 		tissue_tolerated_ambient_pressure[ci] = 0.0;
 	}
diff --git a/dive.h b/dive.h
index e2af00d9d..0cee818fe 100644
--- a/dive.h
+++ b/dive.h
@@ -572,8 +572,9 @@ 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);
-extern void clear_deco(void);
+extern double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds);
+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);
 #ifdef DEBUGFILE
 extern char *debugfilename;
diff --git a/divelist.c b/divelist.c
index 5df8fbb57..32f5f44a3 100644
--- a/divelist.c
+++ b/divelist.c
@@ -818,6 +818,83 @@ static int calculate_sac(struct dive *dive, struct divecomputer *dc)
 	return sac * 1000;
 }
 
+/* for now we do this based on the first divecomputer */
+static void add_dive_to_deco(struct dive *dive)
+{
+	struct divecomputer *dc = &dive->dc;
+	int i;
+
+	if (!dc)
+		return;
+	for (i = 1; i < dive->dc.samples; i++) {
+		struct sample *psample = dc->sample + i - 1;
+		struct sample *sample = dc->sample + i;
+		int t0 = psample->time.seconds;
+		int t1 = sample->time.seconds;
+		int j;
+
+		for (j = t0; j < t1; j++) {
+			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);
+		}
+	}
+}
+
+static struct gasmix air = { .o2.permille = 209 };
+
+/* take into account previous dives until there is a 48h gap between dives */
+void init_decompression(struct dive *dive)
+{
+	int i, divenr = -1;
+	timestamp_t when;
+	gboolean deco_init = FALSE;
+
+	if (!dive)
+		return;
+	while (++divenr < dive_table.nr && get_dive(divenr) != dive)
+		;
+	when = dive->when;
+	i = divenr;
+	while (--i) {
+		struct dive* pdive = get_dive(i);
+		if (!pdive || pdive->when > when || pdive->when + pdive->duration.seconds + 48 * 60 * 60 < when)
+			break;
+		when = pdive->when;
+	}
+
+	while (++i < divenr) {
+		struct dive* pdive = get_dive(i);
+		double surface_pressure = pdive->surface_pressure.mbar ? pdive->surface_pressure.mbar / 1000.0 : 1.013;
+		unsigned int surface_time = get_dive(i+1)->when - pdive->when - pdive->duration.seconds;
+
+		if (!deco_init) {
+			clear_deco(surface_pressure);
+			deco_init = TRUE;
+#if DEBUG & 16
+			dump_tissues();
+#endif
+		}
+		add_dive_to_deco(pdive);
+#if DEBUG & 16
+		printf("added dive #%d\n", pdive->number);
+		dump_tissues();
+#endif
+		add_segment(surface_pressure, &air, surface_time);
+#if DEBUG & 16
+		printf("after surface intervall of %d:%02u\n", FRACTION(surface_time,60));
+		dump_tissues();
+#endif
+	}
+	if (!deco_init) {
+		double surface_pressure = dive->surface_pressure.mbar ? dive->surface_pressure.mbar / 1000.0 : 1.013;
+		clear_deco(surface_pressure);
+#if DEBUG & 16
+		printf("no previous dive\n");
+		dump_tissues();
+#endif
+	}
+}
+
 void update_cylinder_related_info(struct dive *dive)
 {
 	if (dive != NULL) {
diff --git a/divelist.h b/divelist.h
index b309e9a54..41a9f38ae 100644
--- a/divelist.h
+++ b/divelist.h
@@ -15,4 +15,5 @@ extern void remember_tree_state(void);
 extern void restore_tree_state(void);
 extern void select_next_dive(void);
 extern void select_prev_dive(void);
+extern void init_decompression(struct dive * dive);
 #endif
diff --git a/profile.c b/profile.c
index 656e286bd..32a91d2e6 100644
--- a/profile.c
+++ b/profile.c
@@ -1566,7 +1566,8 @@ static struct plot_info *create_plot_info(struct dive *dive, struct divecomputer
 	pi = &gc->pi;
 
 	/* reset deco information to start the calculation */
-	clear_deco();
+	init_decompression(dive);
+
 	/* we want to potentially add synthetic plot_info elements for the gas changes */
 	nr = dc->samples + 4 + 2 * count_gas_change_events(dc);
 	if (last_pi_entry)
@@ -1756,13 +1757,17 @@ static struct plot_info *create_plot_info(struct dive *dive, struct divecomputer
 			float ceiling_pressure = 0;
 			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);
+				double min_pressure = add_segment(depth_to_mbar(depth, dive) / 1000.0,
+								&dive->cylinder[cylinderindex].gasmix, 1);
 				if (min_pressure > ceiling_pressure)
 					ceiling_pressure = min_pressure;
 			}
 			entry->ceiling = deco_allowed_depth(ceiling_pressure, surface_pressure, dive, !prefs.calc_ceiling_3m_incr);
 		}
 	}
+#if DECO_CALC_DEBUG
+	dump_tissues();
+#endif
 
 	if (entry)
 		current->t_end = entry->sec;