Use gradient factors in deco calculation

Usually dive computers show the ceiling in terms of the next deco stop -
and those are in 3m increments. This commit also adds the ability to chose
either the typical 3m increments or the smooth ceiling that the Bühlmann
algorithm actually calculates.

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

5 files changed, 86 insertions, 13 deletions

diff --git a/deco.c b/deco.c
index 6901e6c6a..998d6637c 100644
--- a/deco.c
+++ b/deco.c
@@ -8,6 +8,8 @@
  *
  * clear_deco()  - call to initialize for a new deco calculation
  * add_segment(pressure, gasmix) - add 1 second 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 "dive.h"
 
@@ -24,9 +26,8 @@ struct buehlmann_config {
   double  gf_high_emergency;	//! emergency gf factors
   double  gf_low_emergency;	//! gradient factor low (at bottom/start of deco calculation).
 };
-
-struct dive_data
-{
+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 dive_data {
         double pressure;	//! pesent ambient pressure
         double surface;		//! pressure at water surface
         struct gasmix *gasmix;	//! current selected gas
@@ -75,14 +76,50 @@ 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 DIST_FROM_3_MTR 0.28
+#define PRESSURE_CHANGE_3M 0.3
+#define TOLERANCE 0.02
 
 double tissue_n2_sat[16];
 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;
 
-struct buehlmann_config buehlmann_config = { 1.0, 1.01, 0.5, 3, 95.0, 95.0, 10.0, 30.0, 95.0, 95.0 };
+
+
+static double actual_gradient_limit(const struct dive_data *data)
+{
+	double pressure_diff, limit_at_position;
+	double gf_high = buehlmann_config.gf_high;
+	double gf_low = buehlmann_config.gf_low;
+
+	pressure_diff = data->pressure - data->surface;
+
+	if (pressure_diff > TOLERANCE) {
+		if (pressure_diff < gf_low_position_this_dive)
+			limit_at_position = gf_high - ((gf_high - gf_low) * pressure_diff / gf_low_position_this_dive);
+		else
+			limit_at_position = gf_low;
+	} else {
+		limit_at_position = gf_high;
+	}
+	return limit_at_position;
+}
+
+static double gradient_factor_calculation(const struct dive_data *data)
+{
+	double tissue_inertgas_saturation;
+
+	tissue_inertgas_saturation = tissue_n2_sat[ci_pointing_to_guiding_tissue] +
+					tissue_he_sat[ci_pointing_to_guiding_tissue];
+	if (tissue_inertgas_saturation < data->pressure)
+		return 0.0;
+	else
+		return (tissue_inertgas_saturation - data->pressure) /
+			(tissue_inertgas_saturation - tissue_tolerated_ambient_pressure[ci_pointing_to_guiding_tissue]);
+}
 
 static double tissue_tolerance_calc(void)
 {
@@ -104,7 +141,6 @@ static double tissue_tolerance_calc(void)
 			ret_tolerance_limit_ambient_pressure = tissue_tolerated_ambient_pressure[ci];
 		}
 	}
-	printf("%d:%02u %lf\n",FRACTION(divetime, 60), ret_tolerance_limit_ambient_pressure);
 	return (ret_tolerance_limit_ambient_pressure);
 }
 
@@ -116,7 +152,6 @@ double add_segment(double pressure, struct gasmix *gasmix)
 	double pphe = (pressure - WV_PRESSURE) * gasmix->he.permille / 1000.0;
 
 	divetime++;
-	printf("%2d:%02u N2 %2.3lf He %2.3lf",FRACTION(divetime, 60), ppn2, pphe);
 	/* right now we just do OC */
 	for (ci = 0; ci < 16; ci++) {
 		if (ppn2 - tissue_n2_sat[ci] > 0)
@@ -141,3 +176,36 @@ void clear_deco()
 	}
 	divetime = 0;
 }
+
+unsigned int deco_allowed_depth(double tissues_tolerance, double surface_pressure, struct dive *dive, gboolean smooth)
+{
+	unsigned int depth, multiples_of_3m;
+	gboolean below_gradient_limit;
+	double new_gradient_factor;
+	double pressure_delta = tissues_tolerance - surface_pressure;
+	struct dive_data mydata;
+
+	if (pressure_delta > 0) {
+		if (!smooth) {
+			multiples_of_3m = (pressure_delta + DIST_FROM_3_MTR) / 0.3;
+			depth = 3000 * multiples_of_3m;
+		} else {
+			depth = rel_mbar_to_depth(pressure_delta * 1000, dive);
+		}
+	} else {
+		depth = 0;
+	}
+	mydata.pressure = surface_pressure + depth / 10000.0;
+	mydata.surface = surface_pressure;
+
+	new_gradient_factor = gradient_factor_calculation(&mydata);
+	below_gradient_limit = (new_gradient_factor < actual_gradient_limit(&mydata));
+	while(!below_gradient_limit)
+	{
+		mydata.pressure += PRESSURE_CHANGE_3M;
+		new_gradient_factor = gradient_factor_calculation(&mydata);
+		below_gradient_limit = (new_gradient_factor < actual_gradient_limit(&mydata));
+	}
+
+	return depth;
+}
diff --git a/display-gtk.h b/display-gtk.h
index 74d86ef6d..0523dc0aa 100644
--- a/display-gtk.h
+++ b/display-gtk.h
@@ -41,6 +41,7 @@ struct preferences {
 	partial_pressure_graphs_t pp_graphs;
 	gboolean profile_red_ceiling;
 	gboolean profile_calc_ceiling;
+	gboolean calc_ceiling_3m_incr;
 };
 
 extern struct preferences prefs;
diff --git a/dive.h b/dive.h
index b5d3368a5..e2af00d9d 100644
--- a/dive.h
+++ b/dive.h
@@ -574,7 +574,7 @@ extern void subsurface_command_line_exit(gint *, gchar ***);
 
 extern double add_segment(double pressure, struct gasmix *gasmix);
 extern void clear_deco(void);
-
+extern unsigned int deco_allowed_depth(double tissues_tolerance, double surface_pressure, struct dive *dive, gboolean smooth);
 #ifdef DEBUGFILE
 extern char *debugfilename;
 extern FILE *debugfile;
diff --git a/gtk-gui.c b/gtk-gui.c
index 9637ba532..79003ac2f 100644
--- a/gtk-gui.c
+++ b/gtk-gui.c
@@ -512,6 +512,7 @@ OPTIONCALLBACK(pn2_toggle, prefs.pp_graphs.pn2)
 OPTIONCALLBACK(phe_toggle, prefs.pp_graphs.phe)
 OPTIONCALLBACK(red_ceiling_toggle, prefs.profile_red_ceiling)
 OPTIONCALLBACK(calc_ceiling_toggle, prefs.profile_calc_ceiling)
+OPTIONCALLBACK(calc_ceiling_3m_toggle, prefs.calc_ceiling_3m_incr)
 OPTIONCALLBACK(force_toggle, force_download)
 OPTIONCALLBACK(prefer_dl_toggle, prefer_downloaded)
 
@@ -779,6 +780,11 @@ static void preferences_dialog(GtkWidget *w, gpointer data)
 	gtk_box_pack_start(GTK_BOX(box), button, FALSE, FALSE, 6);
 	g_signal_connect(G_OBJECT(button), "toggled", G_CALLBACK(calc_ceiling_toggle), NULL);
 
+	button = gtk_check_button_new_with_label(_("3m increments for calculated ceiling"));
+	gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(button), prefs.calc_ceiling_3m_incr);
+	gtk_box_pack_start(GTK_BOX(box), button, FALSE, FALSE, 6);
+	g_signal_connect(G_OBJECT(button), "toggled", G_CALLBACK(calc_ceiling_3m_toggle), NULL);
+
 	gtk_widget_show_all(dialog);
 	result = gtk_dialog_run(GTK_DIALOG(dialog));
 	if (result == GTK_RESPONSE_ACCEPT) {
@@ -824,6 +830,7 @@ static void preferences_dialog(GtkWidget *w, gpointer data)
 		subsurface_set_conf("phethreshold", PREF_STRING, phe_threshold_text);
 		subsurface_set_conf("redceiling", PREF_BOOL, BOOL_TO_PTR(prefs.profile_red_ceiling));
 		subsurface_set_conf("calcceiling", PREF_BOOL, BOOL_TO_PTR(prefs.profile_calc_ceiling));
+		subsurface_set_conf("calcceiling3m", PREF_BOOL, BOOL_TO_PTR(prefs.calc_ceiling_3m_incr));
 
 		new_default = strdup(gtk_button_get_label(GTK_BUTTON(xmlfile_button)));
 
@@ -1244,6 +1251,7 @@ void init_ui(int *argcp, char ***argvp)
 	}
 	prefs.profile_red_ceiling = PTR_TO_BOOL(subsurface_get_conf("redceiling", PREF_BOOL));
 	prefs.profile_calc_ceiling = PTR_TO_BOOL(subsurface_get_conf("calcceiling", PREF_BOOL));
+	prefs.calc_ceiling_3m_incr = PTR_TO_BOOL(subsurface_get_conf("calcceiling3m", PREF_BOOL));
 	divelist_font = subsurface_get_conf("divelist_font", PREF_STRING);
 
 	default_filename = subsurface_get_conf("default_filename", PREF_STRING);
diff --git a/profile.c b/profile.c
index 6cc8ce074..656e286bd 100644
--- a/profile.c
+++ b/profile.c
@@ -1560,7 +1560,7 @@ static struct plot_info *create_plot_info(struct dive *dive, struct divecomputer
 	struct plot_data *entry = NULL;
 	struct event *ev;
 	double amb_pressure, po2;
-	int surface_pressure = dive->surface_pressure.mbar ? dive->surface_pressure.mbar : 1013;
+	double surface_pressure = (dive->surface_pressure.mbar ? dive->surface_pressure.mbar : 1013) / 1000.0;
 
 	/* The plot-info is embedded in the graphics context */
 	pi = &gc->pi;
@@ -1760,11 +1760,7 @@ static struct plot_info *create_plot_info(struct dive *dive, struct divecomputer
 				if (min_pressure > ceiling_pressure)
 					ceiling_pressure = min_pressure;
 			}
-			ceiling_pressure = ceiling_pressure * 1000.0 + 0.5;
-			if (ceiling_pressure > surface_pressure)
-				entry->ceiling = rel_mbar_to_depth(ceiling_pressure - surface_pressure, dive);
-			else
-				entry->ceiling = 0;
+			entry->ceiling = deco_allowed_depth(ceiling_pressure, surface_pressure, dive, !prefs.calc_ceiling_3m_incr);
 		}
 	}