CCR patch: Calculate the correct partial gas pressures for CCR dives

This patch adds code to the function fillpressures() in dive.c to
allow calculating o2 pressures, based on the data from the po2
sensors in the system. The following changes were made:
1) add code to perform po2 calculations for CCR with 1, 2 or 3
   oxygen sesnors.
2) Add four fields to the gas_pressures structure in dive.h. This
   allows communication of data between the function that calls
   get_pressures() and the return of partail pressure values to the
   calling function.
3) Delete the fields for setpoint and gas partial pressures from
   the structure plot_info. All partial pressures (from instruments
   as well as calculated) now reside in the pressures structure
   that forms part of plot_info.
4) Perform changes in several parts of profile.c to make use of the
   pressures structure in plot_info.

[Dirk Hohndel: yet again massive whitespace cleanup]

Signed-off-by: willem ferguson <willemferguson@zoology.up.ac.za>
Signed-off-by: Dirk Hohndel <dirk@hohndel.org>

3 files changed, 111 insertions, 27 deletions

diff --git a/dive.c b/dive.c
index 65efd49cf..4d906e5ff 100644
--- a/dive.c
+++ b/dive.c
@@ -1525,23 +1525,105 @@ 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 */
+/* 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 gass pressures applicable to a single point on the dive profile.
+ * The structure "pressures" is used to obtain data and to return calculated gas pressures to the calling software.
+ * Call parameters:	po2 = po2 value applicable to the record in calling function
+ *			amb_pressure = ambient pressure applicable to the record in calling function
+ *			*pressures = structure for communicating o2 sensor values from and gas pressures to the calling function.
+ *			*mix = structure containing cylinder gas mixture information.
+			*dc = pointer to divecomputer structure.
+ * This function called by: calculate_gas_information_new() in profile.c; add_segment() in deco.c.
+ */
 extern void fill_pressures(struct gas_pressures *pressures, const double amb_pressure, const struct gasmix *mix, double po2, const struct divecomputer *dc)
 {
-	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->o2 = po2;
-			pressures->he = (amb_pressure - pressures->o2) * (double)get_he(mix) / (1000 - get_o2(mix));
-			pressures->n2 = amb_pressure - pressures->o2 - pressures->he;
+	double sensor_j, sump, midp, minp, maxp;
+	double diff_limit = 100; // The limit beyond which O2 sensor differences are considered significant (default = 100 mbar)
+	int num_of_diffs;	 // The number of unacceptable differences among the ogygen sensor partial pressure measurements
+	int i, j, np;
+	bool maxdif = false, mindif = false;
+
+	if (dc->dctype == CCR) { // for CCR rebreathers..
+		// Estimate the most reliable PO2, given the different oxygen partial pressure values from the O2 sensors
+		switch (dc->no_o2sensors) {
+		case 2: { // For 2 sensors: take the mean value of the two partial pressure sensors.
+			np = 0;
+			sump = 0;		  // This calculation allows that for some samples (especially at start of dive) with
+			for (j = 0; j < 2; j++) { // an inactive sensor, the sensor(s) with zero values are not used.
+				sensor_j = pressures->sensor[j];
+				if (sensor_j) {
+					np++;
+					sump = sump + sensor_j;
+				}
+			}
+			if (np > 0)			   // if there is at least one sensor value
+				pressures->o2 = sump / np; // then calculate the mean, else
+//			else pressures->o2 = po2; // if there are no valid sensor values, the use the po2 parameter.
+			break;
 		}
-	} else {
+		case 3: { /* For 3 sensors: diff_limit is the critical limit indicating unacceptable difference (default = 100 mbar).
+			   * a) If all three readings are within a range of diff_limit, then take the mean value. This
+			   *     includes the case where reading 1 is within diff_limit of reading 2; and reading 2 is
+			   *     within diff_limit of reading 3, but readings 1 and 3 differ by more than diff_limit.
+			   * b) If one sensor differs by more than diff-limit from the other two, then take the mean
+			   *     of the closer two sensors and disregard the 3rd sensor, considered as an outlier.
+			   * c) If all 3 sensors differ by more than diff_limit then take the mean of the 3 readings. */
+			for (minp = 9999999999, maxp = -1, sump = 0, j = 0; j < 3; j++) {
+				sensor_j = pressures->sensor[j];
+				if (sensor_j < minp)
+					minp = sensor_j;
+				if (sensor_j > maxp)
+					maxp = sensor_j;
+				sump = sump + sensor_j; // Find min, max and mid of p-values
+			}
+			midp = sump - minp - maxp;
+			num_of_diffs = 0;
+			if ((maxp - midp) > diff_limit) {
+				num_of_diffs++;
+				maxdif = true;
+			}
+			if ((midp - minp) > diff_limit) {
+				num_of_diffs++;
+				mindif = true; // Find no of unacceptable differences
+			}
+			switch (num_of_diffs) {
+			case 0:
+				;
+			case 2: {
+				pressures->o2 = sump / 3; // 0 or 2 unacceptable differences: find mean of three values.
+				break;
+			}
+			case 1: {
+				if (maxdif) // 1 unacceptable difference: find mean of two closer values
+					pressures->o2 = (minp + midp) / 2;
+				if (mindif)
+					pressures->o2 = (maxp + midp) / 2;
+				break;
+			}
+			} //switch no_of_diffs
+		}
+		default: { // if # of sensors is 1 or unknown, simply take the value of the first sensor
+			if (pressures->sensor[0])
+				pressures->o2 = pressures->sensor[0];
+			else
+				pressures->o2 = po2; // if no sensor value found, then go to next section, esimating PO2 using depth.
+		}
+		} // switch dc->n0_o2_sensors
+
+		pressures->he = (amb_pressure - pressures->o2) * (double)get_he(mix) / (1000 - get_o2(mix));
+		pressures->n2 = amb_pressure - pressures->o2 - pressures->he;
+
+
+	} // if dc->type == CCR;   Now pressures->o2 should be defined, based on direct measurements.
+
+	if (po2) {						// If we had a CCR dive (above) then pressures->o2 is defined, therefore this option is a fallback,
+		if (po2 >= amb_pressure || get_o2(mix) == 1000) // dependent on the po2 parameter in the call to this
+			pressures->o2 = amb_pressure;		// function and applicable to non-CCR dives.
+		else
+			pressures->o2 = po2;
+		pressures->he = (amb_pressure - pressures->o2) * (double)get_he(mix) / (1000 - get_o2(mix));
+		pressures->n2 = amb_pressure - pressures->o2 - pressures->he;
+	} else if (!pressures->o2) { // Open circuit dives: no gas pressure values available, they need to be calculated
 		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;
diff --git a/dive.h b/dive.h
index 68421fec1..af5049fab 100644
--- a/dive.h
+++ b/dive.h
@@ -135,6 +135,8 @@ static inline int get_he(const struct gasmix *mix)
 
 struct gas_pressures {
 	double o2, n2, he;
+	double sensor[3];
+	double setpoint;
 };
 
 extern void sanitize_gasmix(struct gasmix *mix);
diff --git a/profile.c b/profile.c
index 653c5b9dd..10aa9c898 100644
--- a/profile.c
+++ b/profile.c
@@ -553,10 +553,10 @@ struct plot_data *populate_plot_entries(struct dive *dive, struct divecomputer *
 		entry->in_deco = sample->in_deco;
 		entry->cns = sample->cns;
 		entry->pressures.o2 = sample->po2.mbar / 1000.0;
-		entry->o2setpoint = sample->o2setpoint.mbar / 1000.0;   // for rebreathers
-		entry->o2sensor[0] = sample->o2sensor[0].mbar / 1000.0; // for up to three rebreather O2 sensors
-		entry->o2sensor[1] = sample->o2sensor[1].mbar / 1000.0;
-		entry->o2sensor[2] = sample->o2sensor[2].mbar / 1000.0;
+		entry->pressures.setpoint = sample->o2setpoint.mbar / 1000.0;   // for rebreathers
+		entry->pressures.sensor[0] = sample->o2sensor[0].mbar / 1000.0; // for up to three rebreather O2 sensors
+		entry->pressures.sensor[1] = sample->o2sensor[1].mbar / 1000.0;
+		entry->pressures.sensor[2] = sample->o2sensor[2].mbar / 1000.0;
 
 		/* FIXME! sensor index -> cylinder index translation! */
 		entry->cylinderindex = sample->sensor;
@@ -882,21 +882,21 @@ void fill_o2_values(struct divecomputer *dc, struct plot_info *pi)
 
 		// For 1st iteration, initialise the last_ values
 		if (i == 0) {
-			last_setpoint = pi->entry->o2setpoint;
+			last_setpoint = pi->entry->pressures.setpoint;
 			for (j = 0; j < dc->no_o2sensors; j++)
-				last_sensor[j] = pi->entry->o2sensor[j];
+				last_sensor[j] = pi->entry->pressures.sensor[j];
 		} else {
 			// Now re-insert the missing oxygen pressure values
-			if (entry->o2setpoint)
-				last_setpoint = entry->o2setpoint;
+			if (entry->pressures.setpoint)
+				last_setpoint = entry->pressures.setpoint;
 			else
-				entry->o2setpoint = last_setpoint;
+				entry->pressures.setpoint = last_setpoint;
 
 			for (j = 0; j < dc->no_o2sensors; j++)
-				if (entry->o2sensor[j])
-					last_sensor[j] = entry->o2sensor[j];
+				if (entry->pressures.sensor[j])
+					last_sensor[j] = entry->pressures.sensor[j];
 				else
-					entry->o2sensor[j] = last_sensor[j];
+					entry->pressures.sensor[j] = last_sensor[j];
 		}
 	}
 }
@@ -918,7 +918,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->pressures.o2);
+				entry->o2setpoint, entry->pressures->sensor[0], entry->pressures->sensor[1], entry->pressures->sensor[2], entry->pressures.o2);
 		}
 		fclose(f1);
 	}