Use PSCR gas when computing O2 toxicity

Both the calculations for CNS and OTU did not take
into account the pO2 drop when using a PSCR. Furthermore,
there was some unit confusion due to not using internal
units.

Reported-by: arosl
Signed-off-by: Robert C. Helling <helling@atdotde.de>

4 files changed, 33 insertions, 9 deletions

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 68ae46ab1..c8436d4e1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,4 @@
+- core: Take into account PSCR pO2 drop when computing CNS and OTU
 - planner: allow automatic setpoint changes at specified depths in CCR mode
 - core: fix OTU und CNS calculations for CCR dives at shallow depths
 - desktop: make selection of multiple dives reasonable fast
diff --git a/core/divelist.c b/core/divelist.c
index a6debfb20..379e33eea 100644
--- a/core/divelist.c
+++ b/core/divelist.c
@@ -114,9 +114,16 @@ static int calculate_otu(const struct dive *dive)
 				po2f = MIN((int) sample->setpoint.mbar,
 					   depth_to_mbar(sample->depth.mm, dive));
 			} else {						// For OC and rebreather without o2 sensor/setpoint
-				int o2 = active_o2(dive, dc, psample->time);	// 	... calculate po2 from depth and FiO2.
-				po2i = lrint(o2 * depth_to_bar(psample->depth.mm, dive));	// (initial) po2 at start of segment
-				po2f = lrint(o2 * depth_to_bar(sample->depth.mm, dive));	// (final) po2 at end of segment
+				double amb_presure = depth_to_bar(sample->depth.mm, dive);
+				double pamb_pressure = depth_to_bar(psample->depth.mm , dive);
+				if (dc->divemode == PSCR) {
+					po2i = pscr_o2(pamb_pressure, get_gasmix_at_time(dive, dc, psample->time));
+					po2f = pscr_o2(amb_presure, get_gasmix_at_time(dive, dc, sample->time));
+				} else {
+					int o2 = active_o2(dive, dc, psample->time);	// 	... calculate po2 from depth and FiO2.
+					po2i = lrint(o2 * pamb_pressure);	// (initial) po2 at start of segment
+					po2f = lrint(o2 * amb_presure);	// (final) po2 at end of segment
+				}
 			}
 		}
 		if ((po2i > 500) || (po2f > 500)) {			// If PO2 in segment is above 500 mbar then calculate otu
@@ -173,9 +180,16 @@ static double calculate_cns_dive(const struct dive *dive)
 			trueo2 = true;
 		}
 		if (!trueo2) {
-			int o2 = active_o2(dive, dc, psample->time);			// For OC and rebreather without o2 sensor:
-			po2i = lrint(o2 * depth_to_bar(psample->depth.mm, dive));	// (initial) po2 at start of segment
-			po2f = lrint(o2 * depth_to_bar(sample->depth.mm, dive));	// (final) po2 at end of segment
+			double amb_presure = depth_to_bar(sample->depth.mm, dive);
+			double pamb_pressure = depth_to_bar(psample->depth.mm , dive);
+			if (dc->divemode == PSCR) {
+				po2i = pscr_o2(pamb_pressure, get_gasmix_at_time(dive, dc, psample->time));
+				po2f = pscr_o2(amb_presure, get_gasmix_at_time(dive, dc, sample->time));
+			} else {
+				int o2 = active_o2(dive, dc, psample->time);	// 	... calculate po2 from depth and FiO2.
+				po2i = lrint(o2 * pamb_pressure);	// (initial) po2 at start of segment
+				po2f = lrint(o2 * amb_presure);	// (final) po2 at end of segment
+			}
 		}
 		int po2avg = (po2i + po2f) / 2;	// po2i now holds the mean po2 of initial and final po2 values of segment.
 		/* Don't increase CNS when po2 below 500 matm */
diff --git a/core/gas.c b/core/gas.c
index e62bed1ea..8dbfe61cb 100644
--- a/core/gas.c
+++ b/core/gas.c
@@ -96,6 +96,15 @@ fraction_t get_gas_component_fraction(struct gasmix mix, enum gas_component comp
 	}
 }
 
+// O2 pressure in mbar according to the steady state model for the PSCR
+// NB: Ambient pressure comes in bar!
+int pscr_o2(const double amb_pressure, struct gasmix mix) {
+	int o2 = get_o2(mix) * amb_pressure - (int)((1.0 - get_o2(mix) / 1000.0) * prefs.o2consumption / (prefs.bottomsac * prefs.pscr_ratio) * 1000000);
+	if (o2 < 0.0) // He's dead, Jim.
+		o2 = 0.0;
+	return o2;
+}
+
 /* fill_pressures(): Compute partial gas pressures in bar from gasmix and ambient pressures, possibly for OC or CCR, to be
  * extended to PSCT. This function does the calculations of gas pressures applicable to a single point on the dive profile.
  * The structure "pressures" is used to return calculated gas pressures to the calling software.
@@ -123,9 +132,7 @@ void fill_pressures(struct gas_pressures *pressures, const double amb_pressure,
 		}
 	} else {
 		if (divemode == PSCR) { /* The steady state approximation should be good enough */
-			pressures->o2 = get_o2(mix) / 1000.0 * amb_pressure - (1.0 - get_o2(mix) / 1000.0) * prefs.o2consumption / (prefs.bottomsac * prefs.pscr_ratio / 1000.0);
-			if (pressures->o2 < 0) // He's dead, Jim.
-				pressures->o2 = 0;
+			pressures->o2 = pscr_o2(amb_pressure, mix) / 1000.0;
 			if (get_o2(mix) != 1000) {
 				pressures->he = (amb_pressure - pressures->o2) * get_he(mix) / (1000.0 - get_o2(mix));
 				pressures->n2 = (amb_pressure - pressures->o2) * get_n2(mix) / (1000.0 - get_o2(mix));
diff --git a/core/gas.h b/core/gas.h
index cd7c9e389..83eef72e9 100644
--- a/core/gas.h
+++ b/core/gas.h
@@ -59,6 +59,8 @@ static inline int get_n2(struct gasmix mix)
 	return 1000 - get_o2(mix) - get_he(mix);
 }
 
+int pscr_o2(const double amb_pressure, struct gasmix mix);
+
 struct gas_pressures {
 	double o2, n2, he;
 };