Discontinuity and stall on tank pressure interpolated lines

Interpolated lines for tank pressures were presenting discontinuities
(sudden drops) and stalls (horizontal lines) with certain dive profiles.

The main reason seems to be that the discrete interpolation of tank
pressure was adding small pressure increments that could be rounded down
or up repeatedly generating cumulative rounding errors that would mean
either a delay on pressure drop that would be drawn as a sudden drop or as
a premature pressure drop that would result in a flat line.

This patch changes the way the discrete interpolation is done, so that we
don't have cumulative rounding errors distorting tank pressure lines.

To calculate accumulated pressure_time values the get_pr_interpolate_data
function was created. The fact that get_pr_interpolate_data transverses
the beginning of the plot_info entry list for each entry that needs
interpolated tank pressure isn't optimal at all. There might be a way to
properly track the data necessary to interpolate tank pressures from
inside the main pi->entry loop in fill_missing_tank_pressures.
Unfortunately I didn't manage to do it inside fill_missing_tank_pressures
so we have get_pr_interpolate_data.

The SURFACE_THRESHOLD test from pressure_time function was also removed as
no matter how shallow the diver is, if he is using the cylinder to breathe
the cylinders tank pressure should be affected.

Signed-off-by: Rodrigo Severo <rodrigo@fabricadeideias.com>
Signed-off-by: Dirk Hohndel <dirk@hohndel.org>

1 files changed, 77 insertions, 13 deletions

diff --git a/profile.c b/profile.c
index ea387bf6d..39ad2a1e1 100644
--- a/profile.c
+++ b/profile.c
@@ -417,6 +417,22 @@ static void dump_pr_track(pr_track_t **track_pr)
 }
 #endif
 
+typedef struct pr_interpolate_struct pr_interpolate_t;
+struct pr_interpolate_struct {
+	int start;
+	int end;
+	int pressure_time;
+        int acc_pressure_time;
+};
+
+#ifdef DEBUG_PR_INTERPOLATE
+static void dump_pr_interpolate(pr_interpolate_t interpolate_pr)
+{
+    printf("INTERPOLATE: start %d - end %d - pt %d - acc_pt %d\n",
+            interpolate_pr.start, interpolate_pr.end, interpolate_pr.pressure_time, interpolate_pr.acc_pressure_time);
+}
+#endif
+
 /*
  * This looks at the pressures for one cylinder, and
  * calculates any missing beginning/end pressures for
@@ -473,7 +489,7 @@ static void fill_missing_segment_pressures(pr_track_t *list)
 			if (pt_sum)
 				pressure -= (start-end)*(double)pt/pt_sum;
 			list->end = pressure;
-			if (list == tmp)
+                        if (list == tmp)
 				break;
 			list = list->next;
 			list->start = pressure;
@@ -500,12 +516,57 @@ static inline int pressure_time(struct dive *dive, struct divecomputer *dc, stru
 	int time = b->sec - a->sec;
 	int depth = (a->depth + b->depth)/2;
 
-	if (depth <= SURFACE_THRESHOLD)
-		return 0;
-
 	return depth_to_mbar(depth, dive) * time;
 }
 
+static struct pr_interpolate_struct get_pr_interpolate_data(pr_track_t *segment, struct plot_info *pi, int cur)
+{
+	struct pr_interpolate_struct interpolate;
+	int i;
+	struct plot_data *entry, *cur_entry;
+
+	interpolate.start = segment->start;
+	interpolate.end = segment->end;
+	interpolate.acc_pressure_time = 0;
+	interpolate.pressure_time = 0;
+	cur_entry = pi->entry + cur;
+
+	for (i = 0; i < pi->nr; i++) {
+		entry = pi->entry + i;
+		if (entry->sec < segment->t_start)
+			continue;
+		if (entry->sec >= segment->t_end) {
+			interpolate.pressure_time += entry->pressure_time;
+			break;
+		}
+		if (entry->sec == segment->t_start) {
+			interpolate.acc_pressure_time = 0;
+			interpolate.pressure_time = 0;
+			if (SENSOR_PRESSURE(entry))
+				interpolate.start = SENSOR_PRESSURE(entry);
+		} else if (i < cur) {
+			if (SENSOR_PRESSURE(entry)) {
+				interpolate.start = SENSOR_PRESSURE(entry);
+				interpolate.acc_pressure_time = 0;
+				interpolate.pressure_time = 0;
+			} else {
+				interpolate.acc_pressure_time += entry->pressure_time;
+				interpolate.pressure_time += entry->pressure_time;
+			}
+		} else if (i == cur) {
+			interpolate.acc_pressure_time += entry->pressure_time;
+			interpolate.pressure_time += entry->pressure_time;
+		} else {
+			interpolate.pressure_time += entry->pressure_time;
+			if (SENSOR_PRESSURE(entry)) {
+				interpolate.end = SENSOR_PRESSURE(entry);
+				break;
+			}
+		}
+	}
+	return interpolate;
+}
+
 static void fill_missing_tank_pressures(struct dive *dive, struct plot_info *pi, pr_track_t **track_pr)
 {
 	int cyl, i;
@@ -526,9 +587,9 @@ static void fill_missing_tank_pressures(struct dive *dive, struct plot_info *pi,
 	/* The first two are "fillers", but in case we don't have a sample
 	 * at time 0 we need to process the second of them here */
 	for (i = 1; i < pi->nr; i++) {
-		double magic, cur_pt;
+		double magic;
 		pr_track_t *segment;
-		int pressure;
+                pr_interpolate_t interpolate;
 
 		entry = pi->entry + i;
 		cyl = entry->cylinderindex;
@@ -538,7 +599,7 @@ static void fill_missing_tank_pressures(struct dive *dive, struct plot_info *pi,
 			continue;
 		}
 
-		/* Find the right pressure segment for this entry.. */
+                /* Find the right pressure segment for this entry.. */
 		segment = track_pr[cyl];
 		while (segment && segment->t_end < entry->sec)
 			segment = segment->next;
@@ -549,14 +610,16 @@ static void fill_missing_tank_pressures(struct dive *dive, struct plot_info *pi,
 			continue;
 		}
 
+                interpolate = get_pr_interpolate_data(segment, pi, i);
+#ifdef DEBUG_PR_INTERPOLATE
+                dump_pr_interpolate(interpolate);
+#endif
 		/* Overall pressure change over total pressure-time for this segment*/
-		magic = (segment->end - segment->start) / (double) segment->pressure_time;
+		magic = (interpolate.end - interpolate.start) / (double) interpolate.pressure_time;
 
 		/* Use that overall pressure change to update the current pressure */
-		cur_pt = pressure_time(dive, &dive->dc, entry-1, entry);
-		pressure = cur_pr[cyl] + cur_pt * magic + 0.5;
-		INTERPOLATED_PRESSURE(entry) = pressure;
-		cur_pr[cyl] = pressure;
+		cur_pr[cyl] = interpolate.start + magic * interpolate.acc_pressure_time + 0.5;
+		INTERPOLATED_PRESSURE(entry) = cur_pr[cyl];
 	}
 }
 
@@ -906,7 +969,8 @@ static void populate_pressure_information(struct dive *dive, struct divecomputer
 
 		/* discrete integration of pressure over time to get the SAC rate equivalent */
 		if (current) {
-			current->pressure_time += pressure_time(dive, dc, entry-1, entry);
+                        entry->pressure_time = pressure_time(dive, dc, entry-1, entry);
+			current->pressure_time += entry->pressure_time;
 			current->t_end = entry->sec;
 		}