Use real gas compressibility in planner

Modify formluas for gas use to take into account the
compressibility correction for real gases. This introduces
also the inverse formula to compute the pressure for a given
amount of gas.

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

3 files changed, 19 insertions, 4 deletions

diff --git a/core/dive.h b/core/dive.h
index 8798b22ff..3a1492ed9 100644
--- a/core/dive.h
+++ b/core/dive.h
@@ -136,6 +136,7 @@ extern int units_to_sac(double volume);
 /* Volume in mliter of a cylinder at pressure 'p' */
 extern int gas_volume(cylinder_t *cyl, pressure_t p);
 extern double gas_compressibility_factor(struct gasmix *gas, double bar);
+extern double isothermal_pressure(struct gasmix *gas, double p1, int volume1, int volume2);
 
 
 static inline int get_o2(const struct gasmix *mix)
diff --git a/core/gas-model.c b/core/gas-model.c
index ad1160f3b..98139337c 100644
--- a/core/gas-model.c
+++ b/core/gas-model.c
@@ -62,3 +62,13 @@ double gas_compressibility_factor(struct gasmix *gas, double bar)
 	 */
 	return Z * 0.001 + 1.0;
 }
+
+/* Compute the new pressure when compressing (expanding) volome v1 at pressure p1 bar to volume v2
+ * taking into account the compressebility (to first order) */
+
+double isothermal_pressure(struct gasmix *gas, double p1, int volume1, int volume2)
+{
+	double p_ideal = p1 * volume1 / volume2 / gas_compressibility_factor(gas, p1);
+
+	return p_ideal * gas_compressibility_factor(gas, p_ideal);
+}
diff --git a/core/planner.c b/core/planner.c
index 118ad751d..7bbb7c0c4 100644
--- a/core/planner.c
+++ b/core/planner.c
@@ -246,7 +246,7 @@ static void update_cylinder_pressure(struct dive *d, int old_depth, int new_dept
 	if (in_deco)
 		cyl->deco_gas_used.mliter += gas_used.mliter;
 	if (cyl->type.size.mliter) {
-		delta_p.mbar = gas_used.mliter * 1000.0 / cyl->type.size.mliter;
+		delta_p.mbar = gas_used.mliter * 1000.0 / cyl->type.size.mliter * gas_compressibility_factor(&cyl->gasmix, cyl->end.mbar / 1000.0);
 		cyl->end.mbar -= delta_p.mbar;
 	}
 }
@@ -830,8 +830,11 @@ static void add_plan_to_notes(struct diveplan *diveplan, struct dive *dive, bool
 		volume = get_volume_units(cyl->gas_used.mliter, NULL, &unit);
 		deco_volume = get_volume_units(cyl->deco_gas_used.mliter, NULL, &unit);
 		if (cyl->type.size.mliter) {
-			deco_pressure = get_pressure_units(1000.0 * cyl->deco_gas_used.mliter / cyl->type.size.mliter, &pressure_unit);
-			pressure = get_pressure_units(1000.0 * cyl->gas_used.mliter / cyl->type.size.mliter, &pressure_unit);
+			int remaining_gas = (double)cyl->end.mbar * cyl->type.size.mliter / 1000.0 / gas_compressibility_factor(&cyl->gasmix, cyl->end.mbar / 1000.0);
+			double deco_pressure_bar = isothermal_pressure(&cyl->gasmix, 1.0, remaining_gas + cyl->deco_gas_used.mliter, cyl->type.size.mliter)
+					- cyl->end.mbar / 1000.0;
+			deco_pressure = get_pressure_units(1000.0 * deco_pressure_bar, &pressure_unit);
+			pressure = get_pressure_units(cyl->start.mbar - cyl->end.mbar, &pressure_unit);
 			/* Warn if the plan uses more gas than is available in a cylinder
 			 * This only works if we have working pressure for the cylinder
 			 * 10bar is a made up number - but it seemed silly to pretend you could breathe cylinder down to 0 */
@@ -840,7 +843,8 @@ static void add_plan_to_notes(struct diveplan *diveplan, struct dive *dive, bool
 					translate("gettextFromC", "Warning:"),
 					translate("gettextFromC", "this is more gas than available in the specified cylinder!"));
 			else
-				if ((float) cyl->end.mbar * cyl->type.size.mliter / 1000.0 < (float) cyl->deco_gas_used.mliter)
+				if ((float) cyl->end.mbar * cyl->type.size.mliter / 1000.0 / gas_compressibility_factor(&cyl->gasmix, cyl->end.mbar / 1000.0)
+				    < (float) cyl->deco_gas_used.mliter)
 					snprintf(warning, sizeof(warning), " &mdash; <span style='color: red;'>%s </span> %s",
 						translate("gettextFromC", "Warning:"),
 						translate("gettextFromC", "not enough reserve for gas sharing on ascent!"));