First stab at simplistic dive planning

This comes with absolutely no gui - so the plan literally needs to be
compiled into Subsurface. Not exactly a feature, but this allowed me to
focus on the planning part instead of spending time on tedious UI work.

A new menu "Planner" with entry "Test Planner" calls into the hard-coded
function in planner.c. There a simple dive plan can be constructed with
calls to plan_add_segment(&diveplan, duration, depth at the end, fO2, pO2)

Calling plan(&diveplan) does the deco calculations and creates deco stops
that keep us below the ceiling (with the GFlow/high values currently
configured). The stop levels used are defined at the top of planner.c in
the stoplevels array - there is no need to do the traditional multiples of
3m or anything like that.

The dive including the ascents and deco stops all the way to the surface
is completed and then added as simulated dive to the end of the divelist
(I guess we could automatically select it later) and can be viewed.

This is crude but shows the direction we can go with this. Envision a nice
UI that allows you to simply enter the segments and pick the desired
stops.

What is missing is the ability to give the algorithm additional gases that
it can use during the deco phase - right now it simply keeps using the
last gas used in the diveplan.

All that said, there are clear bugs here - and sadly they seem to be in
the deco calculations, as with the example given the ceiling that is
calculated makes no sense. When displayed in smooth mode it has very
strange jumps up and down that I wouldn't expect. For example with GF
35/75 (the default) the deco ceiling when looking at the simulated dive
jumps from 16m back up to 13m around 14:10 into the dive. That seems very
odd.

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

7 files changed, 269 insertions, 9 deletions

diff --git a/Makefile b/Makefile
index dc3a66aa8..0f6bee13b 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 deco.o \
+OBJS =	main.o dive.o time.o profile.o info.o equipment.o divelist.o deco.o planner.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)
 
@@ -247,6 +247,9 @@ print.o: print.c dive.h display.h display-gtk.h
 deco.o: deco.c dive.h
 	$(CC) $(CFLAGS) $(GLIB2CFLAGS) -c deco.c
 
+planner.o: planner.c dive.h
+	$(CC) $(CFLAGS) $(GLIB2CFLAGS) -c planner.c
+
 libdivecomputer.o: libdivecomputer.c dive.h display.h display-gtk.h libdivecomputer.h
 	$(CC) $(CFLAGS) $(GTK2CFLAGS) $(GLIB2CFLAGS) $(XML2CFLAGS) \
 			$(LIBDIVECOMPUTERCFLAGS) \
diff --git a/dive.h b/dive.h
index 2bc7c8bf8..d58a0146e 100644
--- a/dive.h
+++ b/dive.h
@@ -577,6 +577,9 @@ 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);
+
+extern void test_planner(void);
+
 #ifdef DEBUGFILE
 extern char *debugfilename;
 extern FILE *debugfile;
diff --git a/divelist.c b/divelist.c
index f13f0c48f..fd141ba28 100644
--- a/divelist.c
+++ b/divelist.c
@@ -845,14 +845,15 @@ static void add_dive_to_deco(struct dive *dive)
 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)
+double init_decompression(struct dive *dive)
 {
 	int i, divenr = -1;
 	timestamp_t when;
 	gboolean deco_init = FALSE;
+	double tissue_tolerance;
 
 	if (!dive)
-		return;
+		return 0.0;
 	while (++divenr < dive_table.nr && get_dive(divenr) != dive)
 		;
 	when = dive->when;
@@ -881,7 +882,7 @@ void init_decompression(struct dive *dive)
 		printf("added dive #%d\n", pdive->number);
 		dump_tissues();
 #endif
-		add_segment(surface_pressure, &air, surface_time, 0.0);
+		tissue_tolerance = add_segment(surface_pressure, &air, surface_time, 0.0);
 #if DECO_CALC_DEBUG & 2
 		printf("after surface intervall of %d:%02u\n", FRACTION(surface_time,60));
 		dump_tissues();
@@ -895,6 +896,7 @@ void init_decompression(struct dive *dive)
 		dump_tissues();
 #endif
 	}
+	return tissue_tolerance;
 }
 
 void update_cylinder_related_info(struct dive *dive)
diff --git a/divelist.h b/divelist.h
index 41a9f38ae..1690ec657 100644
--- a/divelist.h
+++ b/divelist.h
@@ -15,5 +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);
+extern double init_decompression(struct dive * dive);
 #endif
diff --git a/gtk-gui.c b/gtk-gui.c
index 42bcb39e6..18c293598 100644
--- a/gtk-gui.c
+++ b/gtk-gui.c
@@ -1117,11 +1117,17 @@ static void next_dc(GtkWidget *w, gpointer data)
 	repaint_dive();
 }
 
+static void test_planner_cb(GtkWidget *w, gpointer data)
+{
+	test_planner();
+}
+
 static GtkActionEntry menu_items[] = {
 	{ "FileMenuAction", NULL, N_("File"), NULL, NULL, NULL},
 	{ "LogMenuAction",  NULL, N_("Log"), NULL, NULL, NULL},
 	{ "ViewMenuAction",  NULL, N_("View"), NULL, NULL, NULL},
 	{ "FilterMenuAction",  NULL, N_("Filter"), NULL, NULL, NULL},
+	{ "PlannerMenuAction",  NULL, N_("Planner"), NULL, NULL, NULL},
 	{ "HelpMenuAction", NULL, N_("Help"), NULL, NULL, NULL},
 	{ "NewFile",        GTK_STOCK_NEW, N_("New"),   CTRLCHAR "N", NULL, G_CALLBACK(file_close) },
 	{ "OpenFile",       GTK_STOCK_OPEN, N_("Open..."),   CTRLCHAR "O", NULL, G_CALLBACK(file_open) },
@@ -1144,6 +1150,7 @@ static GtkActionEntry menu_items[] = {
 	{ "ViewThree",      NULL, N_("Three"), CTRLCHAR "4", NULL, G_CALLBACK(view_three) },
 	{ "PrevDC",         NULL, N_("Prev DC"), NULL, NULL, G_CALLBACK(prev_dc) },
 	{ "NextDC",         NULL, N_("Next DC"), NULL, NULL, G_CALLBACK(next_dc) },
+	{ "TestPlan", NULL, N_("Test Planner"), NULL, NULL, G_CALLBACK(test_planner_cb) }
 };
 static gint nmenu_items = sizeof (menu_items) / sizeof (menu_items[0]);
 
@@ -1192,6 +1199,9 @@ static const gchar* ui_string = " \
 			<menu name=\"FilterMenu\" action=\"FilterMenuAction\"> \
 				<menuitem name=\"SelectEvents\" action=\"SelectEvents\" /> \
 			</menu> \
+			<menu name=\"PlannerMenu\" action=\"PlannerMenuAction\"> \
+				<menuitem name=\"TestPlan\" action=\"TestPlan\" /> \
+			</menu> \
 			<menu name=\"Help\" action=\"HelpMenuAction\"> \
 				<menuitem name=\"About\" action=\"About\" /> \
 			</menu> \
diff --git a/planner.c b/planner.c
new file mode 100644
index 000000000..9d407a7c4
--- /dev/null
+++ b/planner.c
@@ -0,0 +1,242 @@
+/* planner.c
+ *
+ * code that allows us to plan future dives
+ *
+ * (c) Dirk Hohndel 2013
+ */
+
+#include "dive.h"
+#include "divelist.h"
+
+int stoplevels[] = { 3000, 6000, 9000, 12000, 15000, 21000, 30000, 42000, 60000, 90000 };
+
+struct divedatapoint {
+	int time;
+	int depth;
+	int o2;
+	int he;
+	struct divedatapoint *next;
+};
+
+struct diveplan {
+	timestamp_t when;
+	int surface_pressure;
+	struct divedatapoint *dp;
+};
+
+/* returns the tissue tolerance at the end of this (partial) dive */
+double tissue_at_end(struct dive *dive)
+{
+	struct divecomputer *dc;
+	struct sample *sample, *psample;
+	int i, j, t0, t1;
+	double tissue_tolerance;
+
+	if (!dive)
+		return 0.0;
+	tissue_tolerance = init_decompression(dive);
+
+	dc = &dive->dc;
+	if (!dc->samples)
+		return 0.0;
+	psample = sample = dc->sample;
+	t0 = 0;
+	for (i = 0; i < dc->samples; i++, sample++) {
+		t1 = sample->time.seconds;
+		for (j = t0; j < t1; j++) {
+			int depth = psample->depth.mm + (j - t0) * (sample->depth.mm - psample->depth.mm) / (t1 - t0);
+			tissue_tolerance = add_segment(depth_to_mbar(depth, dive) / 1000.0,
+						&dive->cylinder[sample->sensor].gasmix, 1, sample->po2);
+		}
+		psample = sample;
+		t0 = t1;
+	}
+	return tissue_tolerance;
+}
+
+/* how many seconds until we can ascend to the next stop? */
+int time_at_last_depth(struct dive *dive, int next_stop)
+{
+	int depth;
+	double surface_pressure, tissue_tolerance;
+	int wait = 0;
+	struct sample *sample;
+
+	if (!dive)
+		return 0;
+	surface_pressure = dive->surface_pressure.mbar / 1000.0;
+	tissue_tolerance = tissue_at_end(dive);
+	sample = &dive->dc.sample[dive->dc.samples - 1];
+	depth = sample->depth.mm;
+	while (deco_allowed_depth(tissue_tolerance, surface_pressure, dive, 1) > next_stop) {
+		wait++;
+		tissue_tolerance = add_segment(depth_to_mbar(depth, dive) / 1000.0,
+					&dive->cylinder[sample->sensor].gasmix, 1, sample->po2);
+	}
+	return wait;
+}
+
+int add_gas(struct dive *dive, int o2, int he)
+{
+	int i;
+	struct gasmix *mix;
+	cylinder_t *cyl;
+
+	for (i = 0; i < MAX_CYLINDERS; i++) {
+		cyl = dive->cylinder + i;
+		if (cylinder_nodata(cyl))
+			break;
+		mix = &cyl->gasmix;
+		if (o2 == mix->o2.permille && he == mix->he.permille)
+			return i;
+	}
+	if (i == MAX_CYLINDERS) {
+		printf("too many cylinders\n");
+		return -1;
+	}
+	mix = &cyl->gasmix;
+	mix->o2.permille = o2;
+	mix->he.permille = he;
+	return i;
+}
+
+struct dive *create_dive_from_plan(struct diveplan *diveplan)
+{
+	struct dive *dive;
+	struct divedatapoint *dp;
+	struct divecomputer *dc;
+	struct sample *sample;
+	int gasused = 0;
+	int t = 0;
+	int lastdepth = 0;
+
+	if (!diveplan || !diveplan->dp)
+		return NULL;
+	dive = alloc_dive();
+	dive->when = diveplan->when;
+	dive->surface_pressure.mbar = diveplan->surface_pressure;
+	dc = &dive->dc;
+	dc->model = "Simulated Dive";
+	dp = diveplan->dp;
+	while (dp) {
+		int i, depth;
+
+		if (dp->o2 != dive->cylinder[gasused].gasmix.o2.permille ||
+		    dp->he != dive->cylinder[gasused].gasmix.he.permille)
+			gasused = add_gas(dive, dp->o2, dp->he);
+
+		for (i = t; i < dp->time; i += 10) {
+			depth = lastdepth + (i - t) * (dp->depth - lastdepth) / (dp->time - t);
+			sample = prepare_sample(dc);
+			sample->time.seconds = i;
+			sample->depth.mm = depth;
+			sample->sensor = gasused;
+			dc->samples++;
+		}
+		sample = prepare_sample(dc);
+		sample->time.seconds = dp->time;
+		sample->depth.mm = dp->depth;
+		sample->sensor = gasused;
+		lastdepth = dp->depth;
+		t = dp->time;
+		dp = dp->next;
+		dc->samples++;
+	}
+	return dive;
+}
+
+struct divedatapoint *create_dp(int time_incr, int depth, int o2, int he)
+{
+	struct divedatapoint *dp;
+
+	dp = malloc(sizeof(struct divedatapoint));
+	dp->time = time_incr;
+	dp->depth = depth;
+	dp->o2 = o2;
+	dp->he = he;
+	dp->next = NULL;
+	return dp;
+}
+
+void add_to_end_of_diveplan(struct diveplan *diveplan, struct divedatapoint *dp)
+{
+	struct divedatapoint **lastdp = &diveplan->dp;
+	struct divedatapoint *ldp = *lastdp;
+	while(*lastdp) {
+		ldp = *lastdp;
+		lastdp = &(*lastdp)->next;
+	}
+	*lastdp = dp;
+	if (ldp)
+		dp->time += ldp->time;
+}
+
+void plan_add_segment(struct diveplan *diveplan, int duration, int depth, int o2, int he)
+{
+	struct divedatapoint *dp = create_dp(duration, depth, o2, he);
+	add_to_end_of_diveplan(diveplan, dp);
+}
+
+void plan(struct diveplan *diveplan)
+{
+	struct dive *dive;
+	struct sample *sample;
+	int wait_time, o2, he;
+	int ceiling, depth, transitiontime;
+	int stopidx;
+	double tissue_tolerance;
+
+	if (!diveplan->surface_pressure)
+		diveplan->surface_pressure = 1013;
+	dive = create_dive_from_plan(diveplan);
+	record_dive(dive);
+
+	sample = &dive->dc.sample[dive->dc.samples - 1];
+	o2 = dive->cylinder[sample->sensor].gasmix.o2.permille;
+	he = dive->cylinder[sample->sensor].gasmix.he.permille;
+
+	tissue_tolerance = tissue_at_end(dive);
+	ceiling = deco_allowed_depth(tissue_tolerance, diveplan->surface_pressure / 1000.0, dive, 1);
+
+	for (stopidx = 0; stopidx < sizeof(stoplevels) / sizeof(int); stopidx++)
+		if (stoplevels[stopidx] >= ceiling)
+			break;
+
+	while (stopidx >= 0) {
+		depth = dive->dc.sample[dive->dc.samples - 1].depth.mm;
+		if (depth > stoplevels[stopidx]) {
+			transitiontime = (depth - stoplevels[stopidx]) / 150;
+			plan_add_segment(diveplan, transitiontime, stoplevels[stopidx], o2, he);
+			/* re-create the dive */
+			delete_single_dive(dive_table.nr - 1);
+			dive = create_dive_from_plan(diveplan);
+			record_dive(dive);
+		}
+		wait_time = time_at_last_depth(dive, stoplevels[stopidx - 1]);
+		if (wait_time)
+			plan_add_segment(diveplan, wait_time, stoplevels[stopidx], o2, he);
+		transitiontime = (stoplevels[stopidx] - stoplevels[stopidx - 1]) / 150;
+		plan_add_segment(diveplan, transitiontime, stoplevels[stopidx - 1], o2, he);
+		/* re-create the dive */
+		delete_single_dive(dive_table.nr - 1);
+		dive = create_dive_from_plan(diveplan);
+		record_dive(dive);
+		stopidx--;
+	}
+	/* now make the dive visible as last dive of the dive list */
+	report_dives(FALSE, FALSE);
+}
+
+void test_planner()
+{
+	struct diveplan diveplan = {};
+	int end_of_last_dive = dive_table.dives[dive_table.nr - 1]->when + dive_table.dives[dive_table.nr -1]->duration.seconds;
+	diveplan.when = end_of_last_dive + 50 * 3600; /* don't take previous dives into account for deco calculation */
+	diveplan.surface_pressure = 1013;
+	plan_add_segment(&diveplan, 120, 36000, 209, 0);
+	plan_add_segment(&diveplan, 1800, 36000, 209, 0);
+	plan_add_segment(&diveplan, 59, 27000, 209, 0);
+	plan_add_segment(&diveplan, 1, 27000, 400, 0);
+
+	plan(&diveplan);
+}
diff --git a/profile.c b/profile.c
index 2725c27d9..070902a06 100644
--- a/profile.c
+++ b/profile.c
@@ -1755,15 +1755,15 @@ static struct plot_info *create_plot_info(struct dive *dive, struct divecomputer
 			int j;
 			int t0 = (entry - 1)->sec;
 			int t1 = entry->sec;
-			float ceiling_pressure = 0;
+			double tissue_tolerance = 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, entry->po2);
-				if (min_pressure > ceiling_pressure)
-					ceiling_pressure = min_pressure;
+				if (min_pressure > tissue_tolerance)
+					tissue_tolerance = min_pressure;
 			}
-			entry->ceiling = deco_allowed_depth(ceiling_pressure, surface_pressure, dive, !prefs.calc_ceiling_3m_incr);
+			entry->ceiling = deco_allowed_depth(tissue_tolerance, surface_pressure, dive, !prefs.calc_ceiling_3m_incr);
 		}
 	}
 #if DECO_CALC_DEBUG & 1