Fix deco calculations to correctly use GF values and add CC support

The old implementation was broken in several ways.

For one thing the GF values are percentages, so they should normally be
0 < GF < 1 (well, some crazy people like to go above that).

With this most of the Bühlmann config constants were wrong.

Furthermore, after we adjust the pressure tolerance based on the gradient
factors, we need to convert this back into a depth (instead of passing
back the unmodified depth - oops).

Finally, this commit adds closed circuit support to the deco calculations.

Major progress and much more useful at this stage.

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

5 files changed, 34 insertions, 20 deletions

diff --git a/deco.c b/deco.c
index efcdfd39b..8d8fd686c 100644
--- a/deco.c
+++ b/deco.c
@@ -28,7 +28,7 @@ struct buehlmann_config {
   double  gf_high_emergency;	//! emergency gf factors
   double  gf_low_emergency;	//! gradient factor low (at bottom/start of deco calculation).
 };
-struct buehlmann_config buehlmann_config = { 1.0, 1.01, 0.5, 3, 75.0, 35.0, 10.0, 30.0, 95.0, 95.0 };
+struct buehlmann_config buehlmann_config = { 1.0, 1.01, 0.5, 3, 0.75, 0.35, 1.0, 6.0, 0.95, 0.95 };
 struct dive_data {
         double pressure;	//! pesent ambient pressure
         double surface;		//! pressure at water surface
@@ -88,9 +88,6 @@ double tissue_he_sat[16];
 double tissue_tolerated_ambient_pressure[16];
 int ci_pointing_to_guiding_tissue;
 double gf_low_position_this_dive;
-int divetime;
-
-
 
 static double actual_gradient_limit(const struct dive_data *data)
 {
@@ -128,7 +125,7 @@ static double tissue_tolerance_calc(void)
 {
 	int ci = -1;
 	double tissue_inertgas_saturation, buehlmann_inertgas_a, buehlmann_inertgas_b;
-	double ret_tolerance_limit_ambient_pressure = -1.0;
+	double ret_tolerance_limit_ambient_pressure = 0.0;
 
 	for (ci = 0; ci < 16; ci++)
 	{
@@ -144,18 +141,29 @@ static double tissue_tolerance_calc(void)
 			ret_tolerance_limit_ambient_pressure = tissue_tolerated_ambient_pressure[ci];
 		}
 	}
-	return (ret_tolerance_limit_ambient_pressure);
+	return ret_tolerance_limit_ambient_pressure;
 }
 
 /* add a second at the given pressure and gas to the deco calculation */
-double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds)
+double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds, double ccpo2)
 {
 	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;
 
-	/* right now we just do OC */
-	if (period_in_seconds == 1) { /* that's what we do during the dive */
+	if (ccpo2 > 0.0) { /* CC */
+		double rel_o2_amb, f_dilutent;
+		rel_o2_amb = ccpo2 / pressure;
+		f_dilutent = (1 - rel_o2_amb) / (1 - gasmix->o2.permille / 1000.0);
+		if (f_dilutent < 0) { /* setpoint is higher than ambient pressure -> pure O2 */
+			ppn2 = 0.0;
+			pphe = 0.0;
+		} else if (f_dilutent < 1.0) {
+			ppn2 *= f_dilutent;
+			pphe *= f_dilutent;
+		}
+	}
+	if (period_in_seconds == 1) { /* one second interval during 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]) *
@@ -210,7 +218,7 @@ void clear_deco(double surface_pressure)
 		tissue_he_sat[ci] = 0.0;
 		tissue_tolerated_ambient_pressure[ci] = 0.0;
 	}
-	divetime = 0;
+	gf_low_position_this_dive = buehlmann_config.gf_low_position_min;
 }
 
 unsigned int deco_allowed_depth(double tissues_tolerance, double surface_pressure, struct dive *dive, gboolean smooth)
@@ -231,7 +239,7 @@ unsigned int deco_allowed_depth(double tissues_tolerance, double surface_pressur
 	} else {
 		depth = 0;
 	}
-	mydata.pressure = surface_pressure + depth / 10000.0;
+	mydata.pressure = depth_to_mbar(depth, dive) / 1000.0;
 	mydata.surface = surface_pressure;
 
 	new_gradient_factor = gradient_factor_calculation(&mydata);
@@ -242,7 +250,7 @@ unsigned int deco_allowed_depth(double tissues_tolerance, double surface_pressur
 		new_gradient_factor = gradient_factor_calculation(&mydata);
 		below_gradient_limit = (new_gradient_factor < actual_gradient_limit(&mydata));
 	}
-
+	depth = rel_mbar_to_depth((mydata.pressure - surface_pressure) * 1000, dive);
 	return depth;
 }
 
diff --git a/dive.h b/dive.h
index 57ba280e0..2bc7c8bf8 100644
--- a/dive.h
+++ b/dive.h
@@ -572,7 +572,7 @@ 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, int period_in_seconds);
+extern double add_segment(double pressure, struct gasmix *gasmix, int period_in_seconds, double setpoint);
 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);
diff --git a/divelist.c b/divelist.c
index f2e678181..88aff9a49 100644
--- a/divelist.c
+++ b/divelist.c
@@ -834,8 +834,10 @@ static void add_dive_to_deco(struct dive *dive)
 		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);
+			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, sample->po2 / 1000.0);
 		}
 	}
 }
@@ -879,7 +881,7 @@ void init_decompression(struct dive *dive)
 		printf("added dive #%d\n", pdive->number);
 		dump_tissues();
 #endif
-		add_segment(surface_pressure, &air, surface_time);
+		add_segment(surface_pressure, &air, surface_time, 0.0);
 #if DEBUG & 16
 		printf("after surface intervall of %d:%02u\n", FRACTION(surface_time,60));
 		dump_tissues();
diff --git a/gtk-gui.c b/gtk-gui.c
index aad2c6f34..42bcb39e6 100644
--- a/gtk-gui.c
+++ b/gtk-gui.c
@@ -36,7 +36,7 @@ struct preferences prefs = {
 	SI_UNITS,
 	{ TRUE, FALSE, },
 	{ FALSE, FALSE, FALSE, 1.6, 4.0, 13.0},
-	FALSE, FALSE, FALSE, 30.0, 90.0
+	FALSE, FALSE, FALSE, 0.30, 0.75
 };
 
 struct dcnicknamelist {
@@ -821,7 +821,7 @@ static void preferences_dialog(GtkWidget *w, gpointer data)
 	gtk_box_pack_start(GTK_BOX(box), frame, FALSE, FALSE, 6);
 	entry_gflow = gtk_entry_new();
 	gtk_entry_set_max_length(GTK_ENTRY(entry_gflow), 4);
-	snprintf(threshold_text, sizeof(threshold_text), "%.0f", prefs.gflow);
+	snprintf(threshold_text, sizeof(threshold_text), "%.0f", prefs.gflow * 100);
 	gtk_entry_set_text(GTK_ENTRY(entry_gflow), threshold_text);
 	gtk_container_add(GTK_CONTAINER(frame), entry_gflow);
 	gtk_widget_add_events(entry_gflow, GDK_FOCUS_CHANGE_MASK);
@@ -831,7 +831,7 @@ static void preferences_dialog(GtkWidget *w, gpointer data)
 	gtk_box_pack_start(GTK_BOX(box), frame, FALSE, FALSE, 6);
 	entry_gfhigh = gtk_entry_new();
 	gtk_entry_set_max_length(GTK_ENTRY(entry_gfhigh), 4);
-	snprintf(threshold_text, sizeof(threshold_text), "%.0f", prefs.gfhigh);
+	snprintf(threshold_text, sizeof(threshold_text), "%.0f", prefs.gfhigh * 100);
 	gtk_entry_set_text(GTK_ENTRY(entry_gfhigh), threshold_text);
 	gtk_container_add(GTK_CONTAINER(frame), entry_gfhigh);
 	gtk_widget_add_events(entry_gfhigh, GDK_FOCUS_CHANGE_MASK);
@@ -858,6 +858,8 @@ static void preferences_dialog(GtkWidget *w, gpointer data)
 		sscanf(gflow_text, "%lf", &prefs.gflow);
 		gfhigh_text = gtk_entry_get_text(GTK_ENTRY(entry_gfhigh));
 		sscanf(gfhigh_text, "%lf", &prefs.gfhigh);
+		prefs.gflow /= 100.0;
+		prefs.gfhigh /= 100.0;
 		set_gf(prefs.gflow, prefs.gfhigh);
 
 		update_screen();
@@ -1316,12 +1318,14 @@ void init_ui(int *argcp, char ***argvp)
 	conf_value = subsurface_get_conf("gflow", PREF_STRING);
 	if (conf_value) {
 		sscanf(conf_value, "%lf", &prefs.gflow);
+		prefs.gflow /= 100.0;
 		set_gf(prefs.gflow, -1.0);
 		free((void *)conf_value);
 	}
 	conf_value = subsurface_get_conf("gfhigh", PREF_STRING);
 	if (conf_value) {
 		sscanf(conf_value, "%lf", &prefs.gfhigh);
+		prefs.gfhigh /= 100.0;
 		set_gf(-1.0, prefs.gfhigh);
 		free((void *)conf_value);
 	}
diff --git a/profile.c b/profile.c
index 32a91d2e6..6de888607 100644
--- a/profile.c
+++ b/profile.c
@@ -1758,7 +1758,7 @@ static struct plot_info *create_plot_info(struct dive *dive, struct divecomputer
 			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, 1);
+								&dive->cylinder[cylinderindex].gasmix, 1, entry->po2);
 				if (min_pressure > ceiling_pressure)
 					ceiling_pressure = min_pressure;
 			}