Merge branch 'deco2'

Bringing in the first attempts to do our own deco calculations

9 files changed, 478 insertions, 6 deletions

diff --git a/Makefile b/Makefile
index 81e4014fd..dc3a66aa8 100644
--- a/Makefile
+++ b/Makefile
@@ -130,7 +130,7 @@ LIBS = $(LIBXML2) $(LIBXSLT) $(LIBGTK) $(LIBGCONF2) $(LIBDIVECOMPUTER) $(EXTRALI
 MSGLANGS=$(notdir $(wildcard po/*po))
 MSGOBJS=$(addprefix share/locale/,$(MSGLANGS:.po=.UTF-8/LC_MESSAGES/subsurface.mo))
 
-OBJS =	main.o dive.o time.o profile.o info.o equipment.o divelist.o \
+OBJS =	main.o dive.o time.o profile.o info.o equipment.o divelist.o deco.o \
 	parse-xml.o save-xml.o libdivecomputer.o print.o uemis.o uemis-downloader.o \
 	gtk-gui.o statistics.o file.o cochran.o $(OSSUPPORT).o $(RESFILE)
 
@@ -244,6 +244,9 @@ divelist.o: divelist.c dive.h display.h divelist.h
 print.o: print.c dive.h display.h display-gtk.h
 	$(CC) $(CFLAGS) $(GTK2CFLAGS) $(GLIB2CFLAGS) $(XML2CFLAGS) -c print.c
 
+deco.o: deco.c dive.h
+	$(CC) $(CFLAGS) $(GLIB2CFLAGS) -c deco.c
+
 libdivecomputer.o: libdivecomputer.c dive.h display.h display-gtk.h libdivecomputer.h
 	$(CC) $(CFLAGS) $(GTK2CFLAGS) $(GLIB2CFLAGS) $(XML2CFLAGS) \
 			$(LIBDIVECOMPUTERCFLAGS) \
diff --git a/color.h b/color.h
index 23065d306..5d88fb5f7 100644
--- a/color.h
+++ b/color.h
@@ -7,9 +7,11 @@
 // Greens
 #define CAMARONE1    { 0.0, 0.4, 0.0, 1 }
 #define FUNGREEN1    { 0.0, 0.4, 0.2, 1 }
+#define FUNGREEN1_HIGH_TRANS    { 0.0, 0.4, 0.2, 0.25 }
 #define KILLARNEY1   { 0.2, 0.4, 0.2, 1 }
 #define APPLE1       { 0.2, 0.6, 0.2, 1 }
 #define APPLE1_MED_TRANS { 0.2, 0.6, 0.2, 0.5 }
+#define APPLE1_HIGH_TRANS { 0.2, 0.6, 0.2, 0.25 }
 #define LIMENADE1    { 0.4, 0.8, 0.0, 1 }
 #define ATLANTIS1    { 0.4, 0.8, 0.2, 1 }
 #define ATLANTIS2    { 0.6, 0.8, 0.2, 1 }
diff --git a/deco.c b/deco.c
new file mode 100644
index 000000000..efcdfd39b
--- /dev/null
+++ b/deco.c
@@ -0,0 +1,255 @@
+/* calculate deco values
+ * based on Bühlmann ZHL-16b
+ * based on an implemention by heinrichs weikamp for the DR5
+ * the original file doesn't carry a license and is used here with
+ * the permission of Matthias Heinrichs
+ *
+ * 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, 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
+ * set_gf(gflow, gfhigh) - set Buehlmann gradient factors
+ */
+#include <math.h>
+#include "dive.h"
+
+//! Option structure for Buehlmann decompression.
+struct buehlmann_config {
+  double  satmult;		//! safety at inert gas accumulation as percentage of effect (more than 100).
+  double  desatmult;		//! safety at inert gas depletion as percentage of effect (less than 100).
+  double  safety_dist_deco_stop;//! assumed distance to official decompression where decompression takes places.
+  int  last_deco_stop_in_mtr;	//! depth of last_deco_stop.
+  double  gf_high;		//! gradient factor high (at surface).
+  double  gf_low;		//! gradient factor low (at bottom/start of deco calculation).
+  double  gf_low_position_min;	//! gf_low_position below surface_min_shallow.
+  double  gf_low_position_max;	//! gf_low_position below surface_max_depth.
+  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 dive_data {
+        double pressure;	//! pesent ambient pressure
+        double surface;		//! pressure at water surface
+        struct gasmix *gasmix;	//! current selected gas
+};
+
+const double buehlmann_N2_a[] = {1.1696, 1.0, 0.8618, 0.7562,
+				 0.62, 0.5043, 0.441, 0.4,
+				 0.375, 0.35, 0.3295, 0.3065,
+				 0.2835, 0.261, 0.248, 0.2327};
+
+const double buehlmann_N2_b[] = {0.5578, 0.6514, 0.7222, 0.7825,
+				 0.8126, 0.8434, 0.8693, 0.8910,
+				 0.9092, 0.9222, 0.9319, 0.9403,
+				 0.9477, 0.9544, 0.9602, 0.9653};
+
+const double buehlmann_N2_t_halflife[] = {5.0, 8.0, 12.5, 18.5,
+					 27.0, 38.3, 54.3, 77.0,
+					 109.0, 146.0, 187.0, 239.0,
+					 305.0, 390.0, 498.0, 635.0};
+
+const double buehlmann_N2_factor_expositon_one_second[] = {
+	2.30782347297664E-003, 1.44301447809736E-003, 9.23769302935806E-004, 6.24261986779007E-004,
+	4.27777107246730E-004, 3.01585140931371E-004, 2.12729727268379E-004, 1.50020603047807E-004,
+	1.05980191127841E-004, 7.91232600646508E-005, 6.17759153688224E-005, 4.83354552742732E-005,
+	3.78761777920511E-005, 2.96212356654113E-005, 2.31974277413727E-005, 1.81926738960225E-005};
+
+const double buehlmann_He_a[] = { 1.6189, 1.383 , 1.1919, 1.0458,
+				  0.922 , 0.8205, 0.7305, 0.6502,
+				  0.595 , 0.5545, 0.5333, 0.5189,
+				  0.5181, 0.5176, 0.5172, 0.5119};
+
+const double buehlmann_He_b[] = {0.4770, 0.5747, 0.6527, 0.7223,
+				 0.7582, 0.7957, 0.8279, 0.8553,
+				 0.8757, 0.8903, 0.8997, 0.9073,
+				 0.9122, 0.9171, 0.9217, 0.9267};
+
+const double buehlmann_He_t_halflife[] = {1.88, 3.02, 4.72, 6.99,
+					  10.21, 14.48, 20.53, 29.11,
+					  41.20, 55.19, 70.69, 90.34,
+					  115.29, 147.42, 188.24, 240.03};
+
+const double buehlmann_He_factor_expositon_one_second[] = {
+	6.12608039419837E-003, 3.81800836683133E-003, 2.44456078654209E-003, 1.65134647076792E-003,
+	1.13084424730725E-003, 7.97503165599123E-004, 5.62552521860549E-004, 3.96776399429366E-004,
+	2.80360036664540E-004, 2.09299583354805E-004, 1.63410794820518E-004, 1.27869320250551E-004,
+	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
+
+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;
+
+
+
+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)
+{
+	int ci = -1;
+	double tissue_inertgas_saturation, buehlmann_inertgas_a, buehlmann_inertgas_b;
+	double ret_tolerance_limit_ambient_pressure = -1.0;
+
+	for (ci = 0; ci < 16; ci++)
+	{
+		tissue_inertgas_saturation = tissue_n2_sat[ci] + tissue_he_sat[ci];
+		buehlmann_inertgas_a = ((buehlmann_N2_a[ci] * tissue_n2_sat[ci]) + (buehlmann_He_a[ci] * tissue_he_sat[ci])) / tissue_inertgas_saturation;
+		buehlmann_inertgas_b = ((buehlmann_N2_b[ci] * tissue_n2_sat[ci]) + (buehlmann_He_b[ci] * tissue_he_sat[ci])) / tissue_inertgas_saturation;
+
+		tissue_tolerated_ambient_pressure[ci] = (tissue_inertgas_saturation - buehlmann_inertgas_a) * buehlmann_inertgas_b;
+
+		if (tissue_tolerated_ambient_pressure[ci] > ret_tolerance_limit_ambient_pressure)
+		{
+			ci_pointing_to_guiding_tissue = ci;
+			ret_tolerance_limit_ambient_pressure = tissue_tolerated_ambient_pressure[ci];
+		}
+	}
+	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)
+{
+	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 */
+		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 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] = (surface_pressure - WV_PRESSURE) * N2_IN_AIR;
+		tissue_he_sat[ci] = 0.0;
+		tissue_tolerated_ambient_pressure[ci] = 0.0;
+	}
+	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;
+}
+
+void set_gf(double gflow, double gfhigh)
+{
+	if (gflow != -1.0)
+		buehlmann_config.gf_low = gflow;
+	if (gfhigh != -1.0)
+		buehlmann_config.gf_high = gfhigh;
+}
diff --git a/display-gtk.h b/display-gtk.h
index a7279fa9a..629972434 100644
--- a/display-gtk.h
+++ b/display-gtk.h
@@ -40,6 +40,10 @@ struct preferences {
 	visible_cols_t visible_cols;
 	partial_pressure_graphs_t pp_graphs;
 	gboolean profile_red_ceiling;
+	gboolean profile_calc_ceiling;
+	gboolean calc_ceiling_3m_incr;
+	double gflow;
+	double gfhigh;
 };
 
 extern struct preferences prefs;
diff --git a/dive.h b/dive.h
index 9b92baaa7..57ba280e0 100644
--- a/dive.h
+++ b/dive.h
@@ -572,6 +572,11 @@ 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 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);
+extern void set_gf(double gflow, double gfhigh);
 #ifdef DEBUGFILE
 extern char *debugfilename;
 extern FILE *debugfile;
diff --git a/divelist.c b/divelist.c
index e389d038f..f2e678181 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/gtk-gui.c b/gtk-gui.c
index 43e0d14e7..0540de266 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, FALSE, 30.0, 90.0
 };
 
 struct dcnicknamelist {
@@ -511,9 +511,29 @@ OPTIONCALLBACK(po2_toggle, prefs.pp_graphs.po2)
 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)
 
+static void gflow_edit(GtkWidget *w, gpointer _data)
+{
+	double gflow;
+	const char *buf;
+	buf = gtk_entry_get_text(GTK_ENTRY(w));
+	sscanf(buf, "%lf", &gflow);
+	set_gf(prefs.gflow, -1.0);
+}
+
+static void gfhigh_edit(GtkWidget *w, gpointer _data)
+{
+	double gfhigh;
+	const char *buf;
+	buf = gtk_entry_get_text(GTK_ENTRY(w));
+	sscanf(buf, "%lf", &gfhigh);
+	set_gf(-1.0, prefs.gfhigh);
+}
+
 static void event_toggle(GtkWidget *w, gpointer _data)
 {
 	gboolean *plot_ev = _data;
@@ -575,7 +595,7 @@ static void preferences_dialog(GtkWidget *w, gpointer data)
 {
 	int result;
 	GtkWidget *dialog, *notebook, *font, *frame, *box, *vbox, *button, *xmlfile_button;
-	GtkWidget *entry_po2, *entry_pn2, *entry_phe;
+	GtkWidget *entry_po2, *entry_pn2, *entry_phe, *entry_gflow, *entry_gfhigh;
 	const char *current_default, *new_default;
 	char threshold_text[10], utf8_buf[128];
 	struct preferences oldprefs = prefs;
@@ -768,15 +788,49 @@ static void preferences_dialog(GtkWidget *w, gpointer data)
 	box = gtk_hbox_new(FALSE, 6);
 	gtk_container_add(GTK_CONTAINER(vbox), box);
 
-	button = gtk_check_button_new_with_label(_("Show ceiling in red"));
+	button = gtk_check_button_new_with_label(_("Show dc reported ceiling in red"));
 	gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(button), prefs.profile_red_ceiling);
 	gtk_box_pack_start(GTK_BOX(box), button, FALSE, FALSE, 6);
 	g_signal_connect(G_OBJECT(button), "toggled", G_CALLBACK(red_ceiling_toggle), NULL);
 
+	box = gtk_hbox_new(FALSE, 6);
+	gtk_container_add(GTK_CONTAINER(vbox), box);
+
+	button = gtk_check_button_new_with_label(_("Show calculated ceiling"));
+	gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(button), prefs.profile_calc_ceiling);
+	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);
+
+	box = gtk_hbox_new(FALSE, 6);
+	gtk_container_add(GTK_CONTAINER(vbox), box);
+
+	frame = gtk_frame_new(_("GFlow"));
+	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);
+	gtk_entry_set_text(GTK_ENTRY(entry_gflow), threshold_text);
+	gtk_container_add(GTK_CONTAINER(frame), entry_gflow);
+	g_signal_connect(G_OBJECT(entry_gflow), "changed", G_CALLBACK(gflow_edit), NULL);
+
+	frame = gtk_frame_new(_("GFhigh"));
+	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);
+	gtk_entry_set_text(GTK_ENTRY(entry_gfhigh), threshold_text);
+	gtk_container_add(GTK_CONTAINER(frame), entry_gfhigh);
+	g_signal_connect(G_OBJECT(entry_gflow), "changed", G_CALLBACK(gfhigh_edit), NULL);
+
 	gtk_widget_show_all(dialog);
 	result = gtk_dialog_run(GTK_DIALOG(dialog));
 	if (result == GTK_RESPONSE_ACCEPT) {
-		const char *po2_threshold_text, *pn2_threshold_text, *phe_threshold_text;
+		const char *po2_threshold_text, *pn2_threshold_text, *phe_threshold_text, *gflow_text, *gfhigh_text;
 		/* Make sure to flush any modified old dive data with old units */
 		update_dive(NULL);
 
@@ -790,6 +844,11 @@ static void preferences_dialog(GtkWidget *w, gpointer data)
 		sscanf(pn2_threshold_text, "%lf", &prefs.pp_graphs.pn2_threshold);
 		phe_threshold_text = gtk_entry_get_text(GTK_ENTRY(entry_phe));
 		sscanf(phe_threshold_text, "%lf", &prefs.pp_graphs.phe_threshold);
+		gflow_text = gtk_entry_get_text(GTK_ENTRY(entry_gflow));
+		sscanf(gflow_text, "%lf", &prefs.gflow);
+		gfhigh_text = gtk_entry_get_text(GTK_ENTRY(entry_gfhigh));
+		sscanf(gfhigh_text, "%lf", &prefs.gfhigh);
+		set_gf(prefs.gflow, prefs.gfhigh);
 
 		update_screen();
 
@@ -817,6 +876,10 @@ static void preferences_dialog(GtkWidget *w, gpointer data)
 		subsurface_set_conf("pn2threshold", PREF_STRING, pn2_threshold_text);
 		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));
+		subsurface_set_conf("gflow", PREF_STRING, gflow_text);
+		subsurface_set_conf("gfhigh", PREF_STRING, gfhigh_text);
 
 		new_default = strdup(gtk_button_get_label(GTK_BUTTON(xmlfile_button)));
 
@@ -1236,6 +1299,20 @@ void init_ui(int *argcp, char ***argvp)
 		free((void *)conf_value);
 	}
 	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));
+	conf_value = subsurface_get_conf("gflow", PREF_STRING);
+	if (conf_value) {
+		sscanf(conf_value, "%lf", &prefs.gflow);
+		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);
+		set_gf(-1.0, prefs.gfhigh);
+		free((void *)conf_value);
+	}
 	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 dea1ad6ee..32a91d2e6 100644
--- a/profile.c
+++ b/profile.c
@@ -40,6 +40,7 @@ struct plot_data {
 	int temperature;
 	/* Depth info */
 	int depth;
+	int ceiling;
 	int ndl;
 	int stoptime;
 	int stopdepth;
@@ -77,7 +78,7 @@ typedef enum {
 	TEXT_BACKGROUND, ALERT_BG, ALERT_FG, EVENTS, SAMPLE_DEEP, SAMPLE_SHALLOW,
 	SMOOTHED, MINUTE, TIME_GRID, TIME_TEXT, DEPTH_GRID, MEAN_DEPTH, DEPTH_TOP,
 	DEPTH_BOTTOM, TEMP_TEXT, TEMP_PLOT, SAC_DEFAULT, BOUNDING_BOX, PRESSURE_TEXT, BACKGROUND,
-	CEILING_SHALLOW, CEILING_DEEP
+	CEILING_SHALLOW, CEILING_DEEP, CALC_CEILING_SHALLOW, CALC_CEILING_DEEP
 } color_indice_t;
 
 typedef struct {
@@ -135,6 +136,8 @@ static const color_t profile_color[] = {
 	[BACKGROUND]      = {{SPRINGWOOD1, BLACK1_LOW_TRANS}},
 	[CEILING_SHALLOW] = {{REDORANGE1_HIGH_TRANS, REDORANGE1_HIGH_TRANS}},
 	[CEILING_DEEP]    = {{RED1_MED_TRANS, RED1_MED_TRANS}},
+	[CALC_CEILING_SHALLOW] = {{FUNGREEN1_HIGH_TRANS, FUNGREEN1_HIGH_TRANS}},
+	[CALC_CEILING_DEEP]    = {{APPLE1_HIGH_TRANS, APPLE1_HIGH_TRANS}},
 
 };
 
@@ -806,6 +809,24 @@ static void plot_depth_profile(struct graphics_context *gc, struct plot_info *pi
 		cairo_close_path(gc->cr);
 		cairo_fill(gc->cr);
 	}
+	/* finally, plot the calculated ceiling over all this */
+	if (prefs.profile_calc_ceiling) {
+		pat = cairo_pattern_create_linear (0.0, 0.0,  0.0, 256.0 * plot_scale);
+		pattern_add_color_stop_rgba (gc, pat, 0, CALC_CEILING_SHALLOW);
+		pattern_add_color_stop_rgba (gc, pat, 1, CALC_CEILING_DEEP);
+		cairo_set_source(gc->cr, pat);
+		cairo_pattern_destroy(pat);
+		entry = pi->entry;
+		move_to(gc, 0, 0);
+		for (i = 0; i < pi->nr; i++, entry++) {
+			if (entry->ceiling)
+				line_to(gc, entry->sec, entry->ceiling);
+			else
+				line_to(gc, entry->sec, 0);
+		}
+		cairo_close_path(gc->cr);
+		cairo_fill(gc->cr);
+	}
 	/* next show where we have been bad and crossed the ceiling */
 	pat = cairo_pattern_create_linear (0.0, 0.0,  0.0, 256.0 * plot_scale);
 	pattern_add_color_stop_rgba (gc, pat, 0, CEILING_SHALLOW);
@@ -1539,10 +1560,14 @@ static struct plot_info *create_plot_info(struct dive *dive, struct divecomputer
 	struct plot_data *entry = NULL;
 	struct event *ev;
 	double amb_pressure, po2;
+	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;
 
+	/* reset deco information to start the calculation */
+	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)
@@ -1724,7 +1749,25 @@ static struct plot_info *create_plot_info(struct dive *dive, struct divecomputer
 		current->pressure_time += (entry->sec - (entry-1)->sec) *
 			depth_to_mbar((entry->depth + (entry-1)->depth) / 2, dive) / 1000.0;
 		missing_pr |= !SENSOR_PRESSURE(entry);
+		/* and now let's try to do some deco calculations */
+		if (i > 0) {
+			int j;
+			int t0 = (entry - 1)->sec;
+			int t1 = entry->sec;
+			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, 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;
@@ -1965,6 +2008,11 @@ static void plot_string(struct plot_data *entry, char *buf, size_t bufsize,
 		memcpy(buf2, buf, bufsize);
 		snprintf(buf, bufsize, "%s\nT:%.1f %s", buf2, tempvalue, temp_unit);
 	}
+	if (entry->ceiling) {
+		depthvalue = get_depth_units(entry->ceiling, NULL, &depth_unit);
+		memcpy(buf2, buf, bufsize);
+		snprintf(buf, bufsize, "%s\nCalculated ceiling %.0f %s", buf2, depthvalue, depth_unit);
+	}
 	if (entry->stopdepth) {
 		depthvalue = get_depth_units(entry->stopdepth, NULL, &depth_unit);
 		memcpy(buf2, buf, bufsize);