Helper function for partial pressure calculation

This patch introduces a new structure holding partial pressures (doubles in bar) for
all three gases and a helper function to compute them from gasmix (which holds fractions)
and ambient pressure. Currentlty this works for OC and CCR, to be extended later to PSCR.

Currently the dive_comp_type argument is unused.

Signed-off-by: Robert C. Helling <helling@atdotde.de>
Signed-off-by: Dirk Hohndel <dirk@hohndel.org>

8 files changed, 61 insertions, 60 deletions

diff --git a/deco.c b/deco.c
index ffa8ed4ba..2a0ae940c 100644
--- a/deco.c
+++ b/deco.c
@@ -188,28 +188,16 @@ double add_segment(double pressure, const struct gasmix *gasmix, int period_in_s
 {
 	int ci;
 	int fo2 = get_o2(gasmix), fhe = get_he(gasmix);
-	double pn2 = (pressure - WV_PRESSURE) * (1000 - fo2 - fhe) / 1000.0;
-	double phe = (pressure - WV_PRESSURE) * fhe / 1000.0;
+	struct gas_pressures pressures;
+
+	fill_pressures(&pressures, pressure, gasmix, (double) ccpo2 / 1000.0, dive->dc.dctype);
 
 	if (buehlmann_config.gf_low_at_maxdepth && pressure > gf_low_pressure_this_dive)
 		gf_low_pressure_this_dive = pressure;
 
-	if (ccpo2) { /* CC */
-		double rel_o2_amb, f_dilutent;
-		rel_o2_amb = ccpo2 / pressure / 1000;
-		f_dilutent = (1 - rel_o2_amb) / (1 - fo2 / 1000.0);
-		if (f_dilutent < 0) { /* setpoint is higher than ambient pressure -> pure O2 */
-			pn2 = 0.0;
-			phe = 0.0;
-		} else if (f_dilutent < 1.0) {
-			pn2 *= f_dilutent;
-			phe *= f_dilutent;
-		}
-	}
-
 	for (ci = 0; ci < 16; ci++) {
-		double pn2_oversat = pn2 - tissue_n2_sat[ci];
-		double phe_oversat = phe - tissue_he_sat[ci];
+		double pn2_oversat = pressures.n2 - tissue_n2_sat[ci];
+		double phe_oversat = pressures.he - tissue_he_sat[ci];
 		double n2_f = n2_factor(period_in_seconds, ci);
 		double he_f = he_factor(period_in_seconds, ci);
 		double n2_satmult = pn2_oversat > 0 ? buehlmann_config.satmult : buehlmann_config.desatmult;
diff --git a/dive.c b/dive.c
index 22cf3ad46..85f70e3b1 100644
--- a/dive.c
+++ b/dive.c
@@ -1454,6 +1454,29 @@ int gasmix_distance(const struct gasmix *a, const struct gasmix *b)
 	return delta_he + delta_o2;
 }
 
+/* Compute partial gas pressures in bar from gasmix and ambient pressures, possibly for OC or CCR, to be extended to PSCT */
+extern void fill_pressures(struct gas_pressures *pressures, const double amb_pressure, const struct gasmix *mix, double po2, const enum dive_comp_type type)
+{
+	if (po2) {
+		/* we have an O₂ partial pressure in the sample - so this
+		 * is likely a CC dive... use that instead of the value
+		 * from the cylinder info */
+		if (po2 >= amb_pressure || get_o2(mix) == 1000) {
+			pressures->o2 = amb_pressure;
+			pressures->he = 0;
+			pressures->n2 = 0;
+		} else {
+			pressures->he = (amb_pressure - pressures->o2) * (double)get_he(mix) / (1000 - get_o2(mix));
+			pressures->n2 = amb_pressure - pressures->o2 - pressures->he;
+		}
+	} else {
+		pressures->o2 = get_o2(mix) / 1000.0 * amb_pressure;
+		pressures->he = get_he(mix) / 1000.0 * amb_pressure;
+		pressures->n2 = (1000 - get_o2(mix) - get_he(mix)) / 1000.0 * amb_pressure;
+	}
+
+}
+
 static int find_cylinder_match(cylinder_t *cyl, cylinder_t array[], unsigned int used)
 {
 	int i;
diff --git a/dive.h b/dive.h
index 3c32ee8b2..eb8909118 100644
--- a/dive.h
+++ b/dive.h
@@ -115,6 +115,11 @@ static inline int get_he(const struct gasmix *mix)
 	return mix->he.permille;
 }
 
+struct gas_pressures {
+	double o2, n2, he;
+};
+
+extern void fill_pressures(struct gas_pressures *pressures, const double amb_pressure, const struct gasmix *mix, double po2, const enum dive_comp_type type);
 extern void sanitize_gasmix(struct gasmix *mix);
 extern int gasmix_distance(const struct gasmix *a, const struct gasmix *b);
 extern struct gasmix *get_gasmix_from_event(struct event *ev);
diff --git a/profile.c b/profile.c
index 470057237..8656429d5 100644
--- a/profile.c
+++ b/profile.c
@@ -44,8 +44,8 @@ static void dump_pi(struct plot_info *pi)
 		       entry->pressure[0], entry->pressure[1],
 		       entry->sec / 60, entry->sec % 60,
 		       entry->temperature, entry->depth, entry->stopdepth, entry->stoptime, entry->ndl, entry->smoothed,
-		       entry->po2, entry->phe, entry->pn2,
-		       entry->po2 + entry->phe + entry->pn2);
+		       entry->pressures.o2, entry->pressures.he, entry->pressures.n2,
+		       entry->pressures.o2 + entry->pressures.he + entry->pressures.n2);
 	}
 	printf("   }\n");
 }
@@ -550,7 +550,7 @@ struct plot_data *populate_plot_entries(struct dive *dive, struct divecomputer *
 		pi->has_ndl |= sample->ndl.seconds;
 		entry->in_deco = sample->in_deco;
 		entry->cns = sample->cns;
-		entry->po2 = sample->po2.mbar / 1000.0;
+		entry->pressures.o2 = sample->po2.mbar / 1000.0;
 		/* FIXME! sensor index -> cylinder index translation! */
 		entry->cylinderindex = sample->sensor;
 		SENSOR_PRESSURE(entry) = sample->cylinderpressure.mbar;
@@ -703,7 +703,7 @@ static void calculate_ndl_tts(double tissue_tolerance, struct plot_data *entry,
 		while (entry->ndl_calc < max_ndl && deco_allowed_depth(tissue_tolerance, surface_pressure, dive, 1) <= 0) {
 			entry->ndl_calc += time_stepsize;
 			tissue_tolerance = add_segment(depth_to_mbar(entry->depth, dive) / 1000.0,
-						       &dive->cylinder[cylinderindex].gasmix, time_stepsize, entry->po2 * 1000, dive);
+						       &dive->cylinder[cylinderindex].gasmix, time_stepsize, entry->pressures.o2 * 1000, dive);
 		}
 		/* we don't need to calculate anything else */
 		return;
@@ -715,7 +715,7 @@ static void calculate_ndl_tts(double tissue_tolerance, struct plot_data *entry,
 	/* Add segments for movement to stopdepth */
 	for (; ascent_depth > next_stop; ascent_depth -= ascent_mm_per_step, entry->tts_calc += ascent_s_per_step) {
 		tissue_tolerance = add_segment(depth_to_mbar(ascent_depth, dive) / 1000.0,
-					       &dive->cylinder[cylinderindex].gasmix, ascent_s_per_step, entry->po2 * 1000, dive);
+					       &dive->cylinder[cylinderindex].gasmix, ascent_s_per_step, entry->pressures.o2 * 1000, dive);
 		next_stop = ROUND_UP(deco_allowed_depth(tissue_tolerance, surface_pressure, dive, 1), deco_stepsize);
 	}
 	ascent_depth = next_stop;
@@ -733,13 +733,13 @@ static void calculate_ndl_tts(double tissue_tolerance, struct plot_data *entry,
 
 		entry->tts_calc += time_stepsize;
 		tissue_tolerance = add_segment(depth_to_mbar(ascent_depth, dive) / 1000.0,
-					       &dive->cylinder[cylinderindex].gasmix, time_stepsize, entry->po2 * 1000, dive);
+					       &dive->cylinder[cylinderindex].gasmix, time_stepsize, entry->pressures.o2 * 1000, dive);
 
 		if (deco_allowed_depth(tissue_tolerance, surface_pressure, dive, 1) <= next_stop) {
 			/* move to the next stop and add the travel between stops */
 			for (; ascent_depth > next_stop; ascent_depth -= ascent_mm_per_deco_step, entry->tts_calc += ascent_s_per_deco_step)
 				add_segment(depth_to_mbar(ascent_depth, dive) / 1000.0,
-					    &dive->cylinder[cylinderindex].gasmix, ascent_s_per_deco_step, entry->po2 * 1000, dive);
+					    &dive->cylinder[cylinderindex].gasmix, ascent_s_per_deco_step, entry->pressures.o2 * 1000, dive);
 			ascent_depth = next_stop;
 			next_stop -= deco_stepsize;
 		}
@@ -764,7 +764,7 @@ void calculate_deco_information(struct dive *dive, struct divecomputer *dc, stru
 		for (j = t0 + time_stepsize; j <= t1; j += time_stepsize) {
 			int depth = interpolate(entry[-1].depth, entry[0].depth, j - t0, t1 - t0);
 			double min_pressure = add_segment(depth_to_mbar(depth, dive) / 1000.0,
-							  &dive->cylinder[entry->cylinderindex].gasmix, time_stepsize, entry->po2 * 1000, dive);
+							  &dive->cylinder[entry->cylinderindex].gasmix, time_stepsize, entry->pressures.o2 * 1000, dive);
 			tissue_tolerance = min_pressure;
 			if (j - t0 < time_stepsize)
 				time_stepsize = j - t0;
@@ -808,6 +808,7 @@ static void calculate_gas_information_new(struct dive *dive, struct plot_info *p
 {
 	int i;
 	double amb_pressure;
+	struct gas_pressures pressures;
 
 	for (i = 1; i < pi->nr; i++) {
 		int fo2, fhe;
@@ -818,23 +819,7 @@ static void calculate_gas_information_new(struct dive *dive, struct plot_info *p
 		fo2 = get_o2(&dive->cylinder[cylinderindex].gasmix);
 		fhe = get_he(&dive->cylinder[cylinderindex].gasmix);
 
-		if (entry->po2) {
-			/* we have an O₂ partial pressure in the sample - so this
-			 * is likely a CC dive... use that instead of the value
-			 * from the cylinder info */
-			if (entry->po2 >= amb_pressure || fo2 == 1000) {
-				entry->po2 = amb_pressure;
-				entry->phe = 0;
-				entry->pn2 = 0;
-			} else {
-				entry->phe = (amb_pressure - entry->po2) * (double)fhe / (1000 - fo2);
-				entry->pn2 = amb_pressure - entry->po2 - entry->phe;
-			}
-		} else {
-			entry->po2 = fo2 / 1000.0 * amb_pressure;
-			entry->phe = fhe / 1000.0 * amb_pressure;
-			entry->pn2 = (1000 - fo2 - fhe) / 1000.0 * amb_pressure;
-		}
+		fill_pressures(&entry->pressures, amb_pressure, &dive->cylinder[cylinderindex].gasmix, entry->pressures.o2, dive->dc.dctype);
 
 		/* Calculate MOD, EAD, END and EADD based on partial pressures calculated before
 		 * so there is no difference in calculating between OC and CC
@@ -845,9 +830,9 @@ static void calculate_gas_information_new(struct dive *dive, struct plot_info *p
 		entry->end = (entry->depth + 10000) * (1000 - fhe) / 1000.0 - 10000;
 		entry->ead = (entry->depth + 10000) * (1000 - fo2 - fhe) / (double)N2_IN_AIR - 10000;
 		entry->eadd = (entry->depth + 10000) *
-				      (entry->po2 / amb_pressure * O2_DENSITY +
-				       entry->pn2 / amb_pressure * N2_DENSITY +
-				       entry->phe / amb_pressure * HE_DENSITY) /
+				      (entry->pressures.o2 / amb_pressure * O2_DENSITY +
+				       entry->pressures.n2 / amb_pressure * N2_DENSITY +
+				       entry->pressures.he / amb_pressure * HE_DENSITY) /
 				      (O2_IN_AIR * O2_DENSITY + N2_IN_AIR * N2_DENSITY) * 1000 - 10000;
 		if (entry->mod < 0)
 			entry->mod = 0;
@@ -877,7 +862,7 @@ static void debug_print_profiledata(struct plot_info *pi)
 			entry = pi->entry + i;
 			fprintf(f1, "%d gas=%8d %8d ; dil=%8d %8d ; o2_sp= %f %f %f %f PO2= %f\n", i, SENSOR_PRESSURE(entry),
 				INTERPOLATED_PRESSURE(entry), DILUENT_PRESSURE(entry), INTERPOLATED_DILUENT_PRESSURE(entry),
-				entry->o2setpoint, entry->o2sensor[0], entry->o2sensor[1], entry->o2sensor[2], entry->po2);
+				entry->o2setpoint, entry->o2sensor[0], entry->o2sensor[1], entry->o2sensor[2], entry->pressures.o2);
 		}
 		fclose(f1);
 	}
@@ -972,11 +957,11 @@ static void plot_string(struct plot_info *pi, struct plot_data *entry, struct me
 	if (entry->cns)
 		put_format(b, translate("gettextFromC", "CNS: %u%%\n"), entry->cns);
 	if (prefs.pp_graphs.po2)
-		put_format(b, translate("gettextFromC", "pO%s: %.2fbar\n"), UTF8_SUBSCRIPT_2, entry->po2);
+		put_format(b, translate("gettextFromC", "pO%s: %.2fbar\n"), UTF8_SUBSCRIPT_2, entry->pressures.o2);
 	if (prefs.pp_graphs.pn2)
-		put_format(b, translate("gettextFromC", "pN%s: %.2fbar\n"), UTF8_SUBSCRIPT_2, entry->pn2);
+		put_format(b, translate("gettextFromC", "pN%s: %.2fbar\n"), UTF8_SUBSCRIPT_2, entry->pressures.n2);
 	if (prefs.pp_graphs.phe)
-		put_format(b, translate("gettextFromC", "pHe: %.2fbar\n"), entry->phe);
+		put_format(b, translate("gettextFromC", "pHe: %.2fbar\n"), entry->pressures.he);
 	if (prefs.mod) {
 		mod = (int)get_depth_units(entry->mod, NULL, &depth_unit);
 		put_format(b, translate("gettextFromC", "MOD: %d%s\n"), mod, depth_unit);
diff --git a/profile.h b/profile.h
index 0a7abf598..384a17bf3 100644
--- a/profile.h
+++ b/profile.h
@@ -38,7 +38,7 @@ struct plot_data {
 	int cns;
 	int smoothed;
 	int sac;
-	double po2, pn2, phe;
+	struct gas_pressures pressures;
 	double o2setpoint, o2sensor[3]; //for rebreathers with up to 3 PO2 sensors
 	double mod, ead, end, eadd;
 	velocity_t velocity;
diff --git a/qt-ui/profile/diveplotdatamodel.cpp b/qt-ui/profile/diveplotdatamodel.cpp
index 4172d0bae..993c09045 100644
--- a/qt-ui/profile/diveplotdatamodel.cpp
+++ b/qt-ui/profile/diveplotdatamodel.cpp
@@ -47,11 +47,11 @@ QVariant DivePlotDataModel::data(const QModelIndex &index, int role) const
 		case SAC:
 			return item.sac;
 		case PN2:
-			return item.pn2;
+			return item.pressures.n2;
 		case PHE:
-			return item.phe;
+			return item.pressures.he;
 		case PO2:
-			return item.po2;
+			return item.pressures.o2;
 		case HEARTBEAT:
 			return item.heartbeat;
 		}
@@ -156,15 +156,15 @@ unsigned int DivePlotDataModel::dcShown() const
 	{                                                             \
 		double ret = -1;                                      \
 		for (int i = 0, count = rowCount(); i < count; i++) { \
-			if (pInfo.entry[i].GAS > ret)                 \
-				ret = pInfo.entry[i].GAS;             \
+			if (pInfo.entry[i].pressures.GAS > ret)                 \
+				ret = pInfo.entry[i].pressures.GAS;             \
 		}                                                     \
 		return ret;                                           \
 	}
 
-MAX_PPGAS_FUNC(phe, pheMax);
-MAX_PPGAS_FUNC(pn2, pn2Max);
-MAX_PPGAS_FUNC(po2, po2Max);
+MAX_PPGAS_FUNC(he, pheMax);
+MAX_PPGAS_FUNC(n2, pn2Max);
+MAX_PPGAS_FUNC(o2, po2Max);
 
 void DivePlotDataModel::emitDataChanged()
 {
diff --git a/qt-ui/profile/diveprofileitem.cpp b/qt-ui/profile/diveprofileitem.cpp
index c5321d770..9ff30ff32 100644
--- a/qt-ui/profile/diveprofileitem.cpp
+++ b/qt-ui/profile/diveprofileitem.cpp
@@ -5,8 +5,8 @@
 #include "divetextitem.h"
 #include "profilewidget2.h"
 #include "animationfunctions.h"
-#include "profile.h"
 #include "dive.h"
+#include "profile.h"
 #include "preferences.h"
 #include "helpers.h"
 #include "diveplanner.h"
diff --git a/qt-ui/profile/divetooltipitem.cpp b/qt-ui/profile/divetooltipitem.cpp
index 40c556cc6..9ee44e4b6 100644
--- a/qt-ui/profile/divetooltipitem.cpp
+++ b/qt-ui/profile/divetooltipitem.cpp
@@ -1,8 +1,8 @@
 #include "divetooltipitem.h"
 #include "divecartesianaxis.h"
 #include "profilewidget2.h"
-#include "profile.h"
 #include "dive.h"
+#include "profile.h"
 #include "membuffer.h"
 #include <QPropertyAnimation>
 #include <QGraphicsSceneMouseEvent>