Merge branch 'for_dirk' of https://github.com/atdotde/subsurface

4 files changed, 148 insertions, 73 deletions

diff --git a/deco.c b/deco.c
index 86acc0351..506d8b58f 100644
--- a/deco.c
+++ b/deco.c
@@ -244,7 +244,7 @@ double tissue_tolerance_calc(const struct dive *dive, double pressure)
 	double lowest_ceiling = 0.0;
 	double tissue_lowest_ceiling[16];
 
-	if (prefs.deco_mode != VPMB || !in_planner()) {
+	if (prefs.deco_mode != VPMB) {
 		for (ci = 0; ci < 16; ci++) {
 			tissue_inertgas_saturation[ci] = tissue_n2_sat[ci] + tissue_he_sat[ci];
 			buehlmann_inertgas_a[ci] = ((buehlmann_N2_a[ci] * tissue_n2_sat[ci]) + (buehlmann_He_a[ci] * tissue_he_sat[ci])) / tissue_inertgas_saturation[ci];
@@ -509,7 +509,7 @@ void add_segment(double pressure, const struct gasmix *gasmix, int period_in_sec
 		tissue_n2_sat[ci] += n2_satmult * pn2_oversat * n2_f;
 		tissue_he_sat[ci] += he_satmult * phe_oversat * he_f;
 	}
-	if(prefs.deco_mode == VPMB && in_planner())
+	if(prefs.deco_mode == VPMB)
 		calc_crushing_pressure(pressure);
 	return;
 }
diff --git a/profile.c b/profile.c
index d39133c21..edba51ebd 100644
--- a/profile.c
+++ b/profile.c
@@ -28,6 +28,9 @@ unsigned int dc_number = 0;
 static struct plot_data *last_pi_entry_new = NULL;
 void populate_pressure_information(struct dive *, struct divecomputer *, struct plot_info *, int);
 
+extern bool in_planner();
+extern pressure_t first_ceiling_pressure;
+
 #ifdef DEBUG_PI
 /* debugging tool - not normally used */
 static void dump_pi(struct plot_info *pi)
@@ -866,6 +869,8 @@ static void calculate_ndl_tts(struct plot_data *entry, struct dive *dive, double
 	int ascent_depth = entry->depth;
 	/* at what time should we give up and say that we got enuff NDL? */
 	int cylinderindex = entry->cylinderindex;
+	/* If iterating through a dive, entry->tts_calc needs to be reset */
+	entry->tts_calc = 0;
 
 	/* If we don't have a ceiling yet, calculate ndl. Don't try to calculate
 	 * a ndl for lower values than 3m it would take forever */
@@ -927,69 +932,133 @@ static void calculate_ndl_tts(struct plot_data *entry, struct dive *dive, double
  */
 void calculate_deco_information(struct dive *dive, struct divecomputer *dc, struct plot_info *pi, bool print_mode)
 {
-	int i;
+	int i, count_iteration = 0;
 	double surface_pressure = (dc->surface_pressure.mbar ? dc->surface_pressure.mbar : get_surface_pressure_in_mbar(dive, true)) / 1000.0;
 	int last_ndl_tts_calc_time = 0;
-	for (i = 1; i < pi->nr; i++) {
-		struct plot_data *entry = pi->entry + i;
-		int j, t0 = (entry - 1)->sec, t1 = entry->sec;
-		int time_stepsize = 20;
-
-		entry->ambpressure = depth_to_bar(entry->depth, dive);
-		entry->gfline = MAX((double)prefs.gflow, (entry->ambpressure - surface_pressure) / (gf_low_pressure_this_dive - surface_pressure) *
-									 (prefs.gflow - prefs.gfhigh) +
-								 prefs.gfhigh) *
-					(100.0 - AMB_PERCENTAGE) / 100.0 + AMB_PERCENTAGE;
-		if (t0 > t1) {
-			fprintf(stderr, "non-monotonous dive stamps %d %d\n", t0, t1);
-			int xchg = t1;
-			t1 = t0;
-			t0 = xchg;
-		}
-		if (t0 != t1 && t1 - t0 < time_stepsize)
-			time_stepsize = t1 - t0;
-		for (j = t0 + time_stepsize; j <= t1; j += time_stepsize) {
-			int depth = interpolate(entry[-1].depth, entry[0].depth, j - t0, t1 - t0);
-			add_segment(depth_to_bar(depth, dive),
-				    &dive->cylinder[entry->cylinderindex].gasmix, time_stepsize, entry->o2pressure.mbar, dive, entry->sac);
-			if ((t1 - j < time_stepsize) && (j < t1))
-				time_stepsize = t1 - j;
-		}
-		if (t0 == t1)
-			entry->ceiling = (entry - 1)->ceiling;
-		else
-			entry->ceiling = deco_allowed_depth(tissue_tolerance_calc(dive, depth_to_bar(entry->depth, dive)), surface_pressure, dive, !prefs.calcceiling3m);
-		for (j = 0; j < 16; j++) {
-			double m_value = buehlmann_inertgas_a[j] + entry->ambpressure / buehlmann_inertgas_b[j];
-			entry->ceilings[j] = deco_allowed_depth(tolerated_by_tissue[j], surface_pressure, dive, 1);
-			entry->percentages[j] = tissue_inertgas_saturation[j] < entry->ambpressure ?
-							tissue_inertgas_saturation[j] / entry->ambpressure * AMB_PERCENTAGE :
-							AMB_PERCENTAGE + (tissue_inertgas_saturation[j] - entry->ambpressure) / (m_value - entry->ambpressure) * (100.0 - AMB_PERCENTAGE);
-		}
+	int first_ceiling = 0;
+	bool first_iteration = true;
+	int final_tts = 0 , time_clear_ceiling = 0, time_deep_ceiling = 0, deco_time = 0, prev_deco_time = 10000000;
+	char *cache_data_initial = NULL;
+	/* For VPM-B outside the planner, cache the initial deco state for CVA iterations */
+	if (prefs.deco_mode == VPMB && !in_planner())
+		cache_deco_state(&cache_data_initial);
+	/* For VPM-B outside the planner, iterate until deco time converges (usually one or two iterations after the initial)
+	 * Set maximum number of iterations to 10 just in case */
+	while ((abs(prev_deco_time - deco_time) >= 30) && (count_iteration < 10)) {
+		for (i = 1; i < pi->nr; i++) {
+			struct plot_data *entry = pi->entry + i;
+			int j, t0 = (entry - 1)->sec, t1 = entry->sec;
+			int time_stepsize = 20;
+
+			entry->ambpressure = depth_to_bar(entry->depth, dive);
+			entry->gfline = MAX((double)prefs.gflow, (entry->ambpressure - surface_pressure) / (gf_low_pressure_this_dive - surface_pressure) *
+										(prefs.gflow - prefs.gfhigh) +
+									prefs.gfhigh) *
+						(100.0 - AMB_PERCENTAGE) / 100.0 + AMB_PERCENTAGE;
+			if (t0 > t1) {
+				fprintf(stderr, "non-monotonous dive stamps %d %d\n", t0, t1);
+				int xchg = t1;
+				t1 = t0;
+				t0 = xchg;
+			}
+			if (t0 != t1 && t1 - t0 < time_stepsize)
+				time_stepsize = t1 - t0;
+			for (j = t0 + time_stepsize; j <= t1; j += time_stepsize) {
+				int depth = interpolate(entry[-1].depth, entry[0].depth, j - t0, t1 - t0);
+				add_segment(depth_to_bar(depth, dive),
+					&dive->cylinder[entry->cylinderindex].gasmix, time_stepsize, entry->o2pressure.mbar, dive, entry->sac);
+				if ((t1 - j < time_stepsize) && (j < t1))
+					time_stepsize = t1 - j;
+			}
+			if (t0 == t1) {
+				entry->ceiling = (entry - 1)->ceiling;
+			} else {
+				/* Keep updating the VPM-B gradients until the start of the ascent phase of the dive. */
+				if (prefs.deco_mode == VPMB && !in_planner() && (entry - 1)->ceiling >= first_ceiling && first_iteration == true) {
+					nuclear_regeneration(t1);
+					vpmb_start_gradient();
+					/* For CVA calculations, start by guessing deco time = dive time remaining */
+					deco_time = pi->maxtime - t1;
+					vpmb_next_gradient(deco_time, surface_pressure / 1000.0);
+				}
+				entry->ceiling = deco_allowed_depth(tissue_tolerance_calc(dive, depth_to_bar(entry->depth, dive)), surface_pressure, dive, !prefs.calcceiling3m);
+				/* If using VPM-B outside the planner, take first_ceiling_pressure as the deepest ceiling */
+				if (prefs.deco_mode == VPMB && !in_planner()) {
+					if  (entry->ceiling >= first_ceiling) {
+								time_deep_ceiling = t1;
+								first_ceiling = entry->ceiling;
+								first_ceiling_pressure.mbar = depth_to_mbar(first_ceiling, dive);
+								if (first_iteration) {
+									nuclear_regeneration(t1);
+									vpmb_start_gradient();
+									/* For CVA calculations, start by guessing deco time = dive time remaining */
+									deco_time = pi->maxtime - t1;
+									vpmb_next_gradient(deco_time, surface_pressure / 1000.0);
+								}
+					}
+					// Use the point where the ceiling clears as the end of deco phase for CVA calculations
+					if (entry->ceiling > 0)
+						time_clear_ceiling = 0;
+					else if (time_clear_ceiling == 0)
+						time_clear_ceiling = t1;
+				}
+			}
+			for (j = 0; j < 16; j++) {
+				double m_value = buehlmann_inertgas_a[j] + entry->ambpressure / buehlmann_inertgas_b[j];
+				entry->ceilings[j] = deco_allowed_depth(tolerated_by_tissue[j], surface_pressure, dive, 1);
+				entry->percentages[j] = tissue_inertgas_saturation[j] < entry->ambpressure ?
+								tissue_inertgas_saturation[j] / entry->ambpressure * AMB_PERCENTAGE :
+								AMB_PERCENTAGE + (tissue_inertgas_saturation[j] - entry->ambpressure) / (m_value - entry->ambpressure) * (100.0 - AMB_PERCENTAGE);
+			}
 
-		/* should we do more calculations?
-		 * We don't for print-mode because this info doesn't show up there */
-		if (prefs.calcndltts && !print_mode) {
-			/* only calculate ndl/tts on every 30 seconds */
-			if ((entry->sec - last_ndl_tts_calc_time) < 30) {
-				struct plot_data *prev_entry = (entry - 1);
-				entry->stoptime_calc = prev_entry->stoptime_calc;
-				entry->stopdepth_calc = prev_entry->stopdepth_calc;
-				entry->tts_calc = prev_entry->tts_calc;
-				entry->ndl_calc = prev_entry->ndl_calc;
-				continue;
+			/* should we do more calculations?
+			* We don't for print-mode because this info doesn't show up there
+			* If the ceiling hasn't cleared by the last data point, we need tts for VPM-B CVA calculation
+			* It is not necessary to do these calculation on the first VPMB iteration, except for the last data point */
+			if ((prefs.calcndltts && !print_mode && (prefs.deco_mode != VPMB || in_planner() || !first_iteration)) ||
+			    (prefs.deco_mode == VPMB && !in_planner() && i == pi->nr - 1)) {
+				/* only calculate ndl/tts on every 30 seconds */
+				if ((entry->sec - last_ndl_tts_calc_time) < 30 && i != pi->nr - 1) {
+					struct plot_data *prev_entry = (entry - 1);
+					entry->stoptime_calc = prev_entry->stoptime_calc;
+					entry->stopdepth_calc = prev_entry->stopdepth_calc;
+					entry->tts_calc = prev_entry->tts_calc;
+					entry->ndl_calc = prev_entry->ndl_calc;
+					continue;
+				}
+				last_ndl_tts_calc_time = entry->sec;
+
+				/* We are going to mess up deco state, so store it for later restore */
+				char *cache_data = NULL;
+				cache_deco_state(&cache_data);
+				calculate_ndl_tts(entry, dive, surface_pressure);
+				if (prefs.deco_mode == VPMB && !in_planner() && i == pi->nr - 1)
+					final_tts = entry->tts_calc;
+				/* Restore "real" deco state for next real time step */
+				restore_deco_state(cache_data);
+				free(cache_data);
 			}
-			last_ndl_tts_calc_time = entry->sec;
-
-			/* We are going to mess up deco state, so store it for later restore */
-			char *cache_data = NULL;
-			cache_deco_state(&cache_data);
-			calculate_ndl_tts(entry, dive, surface_pressure);
-			/* Restore "real" deco state for next real time step */
-			restore_deco_state(cache_data);
-			free(cache_data);
+		}
+		if (prefs.deco_mode == VPMB && !in_planner()) {
+			prev_deco_time = deco_time;
+			// Do we need to update deco_time?
+			if (final_tts > 0)
+				deco_time = pi->maxtime + final_tts - time_deep_ceiling;
+			else if (time_clear_ceiling > 0)
+				deco_time = time_clear_ceiling - time_deep_ceiling;
+			vpmb_next_gradient(deco_time, surface_pressure / 1000.0);
+			final_tts = 0;
+			last_ndl_tts_calc_time = 0;
+			first_ceiling = 0;
+			first_iteration = false;
+			count_iteration ++;
+			restore_deco_state(cache_data_initial);
+		} else {
+			// With Buhlmann, or not in planner, iterating isn't needed.  This makes the while condition false.
+			prev_deco_time = deco_time = 0;
 		}
 	}
+	free(cache_data_initial);
 #if DECO_CALC_DEBUG & 1
 	dump_tissues();
 #endif
diff --git a/qt-ui/profile/profilewidget2.cpp b/qt-ui/profile/profilewidget2.cpp
index 3ccd1bb6d..d47467038 100644
--- a/qt-ui/profile/profilewidget2.cpp
+++ b/qt-ui/profile/profilewidget2.cpp
@@ -87,7 +87,7 @@ ProfileWidget2::ProfileWidget2(QWidget *parent) : QGraphicsView(parent),
 	gasPressureItem(new DiveGasPressureItem()),
 	diveComputerText(new DiveTextItem()),
 	diveCeiling(new DiveCalculatedCeiling()),
-	gradientFactor(new DiveTextItem()),
+	decoModelParameters(new DiveTextItem()),
 	reportedCeiling(new DiveReportedCeiling()),
 	pn2GasItem(new PartialPressureGasItem()),
 	pheGasItem(new PartialPressureGasItem()),
@@ -204,7 +204,7 @@ void ProfileWidget2::addItemsToScene()
 	diveComputerText->setData(SUBSURFACE_OBJ_DATA, SUBSURFACE_OBJ_DC_TEXT);
 	scene()->addItem(diveComputerText);
 	scene()->addItem(diveCeiling);
-	scene()->addItem(gradientFactor);
+	scene()->addItem(decoModelParameters);
 	scene()->addItem(reportedCeiling);
 	scene()->addItem(pn2GasItem);
 	scene()->addItem(pheGasItem);
@@ -288,11 +288,11 @@ void ProfileWidget2::setupItemOnScene()
 	rulerItem->setAxis(timeAxis, profileYAxis);
 	tankItem->setHorizontalAxis(timeAxis);
 
-	// show the gradient factor at the top in the center
-	gradientFactor->setY(0);
-	gradientFactor->setX(50);
-	gradientFactor->setBrush(getColor(PRESSURE_TEXT));
-	gradientFactor->setAlignment(Qt::AlignHCenter | Qt::AlignBottom);
+	// show the deco model parameters at the top in the center
+	decoModelParameters->setY(0);
+	decoModelParameters->setX(50);
+	decoModelParameters->setBrush(getColor(PRESSURE_TEXT));
+	decoModelParameters->setAlignment(Qt::AlignHCenter | Qt::AlignBottom);
 
 	setupItem(reportedCeiling, timeAxis, profileYAxis, dataModel, DivePlotDataModel::CEILING, DivePlotDataModel::TIME, 1);
 	setupItem(diveCeiling, timeAxis, profileYAxis, dataModel, DivePlotDataModel::CEILING, DivePlotDataModel::TIME, 1);
@@ -509,7 +509,10 @@ void ProfileWidget2::plotDive(struct dive *d, bool force)
 
 		// this copies the dive and makes copies of all the relevant additional data
 		copy_dive(d, &displayed_dive);
-		gradientFactor->setText(QString("GF %1/%2").arg(prefs.gflow).arg(prefs.gfhigh));
+		if (prefs.deco_mode == VPMB)
+			decoModelParameters->setText(QString("VPM-B +%1").arg(prefs.conservatism_level));
+		else
+			decoModelParameters->setText(QString("GF %1/%2").arg(prefs.gflow).arg(prefs.gfhigh));
 	} else {
 		DivePlannerPointsModel *plannerModel = DivePlannerPointsModel::instance();
 		plannerModel->createTemporaryPlan();
@@ -518,7 +521,10 @@ void ProfileWidget2::plotDive(struct dive *d, bool force)
 			plannerModel->deleteTemporaryPlan();
 			return;
 		}
-		gradientFactor->setText(QString("GF %1/%2").arg(diveplan.gflow).arg(diveplan.gfhigh));
+		if (prefs.deco_mode == VPMB)
+			decoModelParameters->setText(QString("VPM-B +%1").arg(prefs.conservatism_level));
+		else
+			decoModelParameters->setText(QString("GF %1/%2").arg(prefs.gflow).arg(prefs.gfhigh));
 	}
 
 	// special handling for the first time we display things
@@ -926,7 +932,7 @@ void ProfileWidget2::setEmptyState()
 	toolTipItem->setVisible(false);
 	diveComputerText->setVisible(false);
 	diveCeiling->setVisible(false);
-	gradientFactor->setVisible(false);
+	decoModelParameters->setVisible(false);
 	reportedCeiling->setVisible(false);
 	rulerItem->setVisible(false);
 	tankItem->setVisible(false);
@@ -1053,7 +1059,7 @@ void ProfileWidget2::setProfileState()
 	diveComputerText->setPos(itemPos.dcLabel.on);
 
 	diveCeiling->setVisible(prefs.calcceiling);
-	gradientFactor->setVisible(prefs.calcceiling);
+	decoModelParameters->setVisible(prefs.calcceiling);
 	reportedCeiling->setVisible(prefs.dcceiling);
 
 	if (prefs.calcalltissues) {
@@ -1131,7 +1137,7 @@ void ProfileWidget2::setAddState()
 	/* show the same stuff that the profile shows. */
 	currentState = ADD; /* enable the add state. */
 	diveCeiling->setVisible(true);
-	gradientFactor->setVisible(true);
+	decoModelParameters->setVisible(true);
 	setBackgroundBrush(QColor("#A7DCFF"));
 }
 
@@ -1165,7 +1171,7 @@ void ProfileWidget2::setPlanState()
 	/* show the same stuff that the profile shows. */
 	currentState = PLAN; /* enable the add state. */
 	diveCeiling->setVisible(true);
-	gradientFactor->setVisible(true);
+	decoModelParameters->setVisible(true);
 	setBackgroundBrush(QColor("#D7E3EF"));
 }
 
diff --git a/qt-ui/profile/profilewidget2.h b/qt-ui/profile/profilewidget2.h
index 2d1a7bfb4..23e939869 100644
--- a/qt-ui/profile/profilewidget2.h
+++ b/qt-ui/profile/profilewidget2.h
@@ -171,7 +171,7 @@ private:
 	QList<DiveEventItem *> eventItems;
 	DiveTextItem *diveComputerText;
 	DiveCalculatedCeiling *diveCeiling;
-	DiveTextItem *gradientFactor;
+	DiveTextItem *decoModelParameters;
 	QList<DiveCalculatedTissue *> allTissues;
 	DiveReportedCeiling *reportedCeiling;
 	PartialPressureGasItem *pn2GasItem;