Take incompressibility of gas into account at higher pressures

This creates a helper function called "gas_volume()" that takes the
cylinder and a particular pressure, and returns the estimated volume of
the gas at surface pressure, including proper approximation of the
incompressibility of gas.

It very much is an approximation, but it's closer to reality than
assuming a pure ideal gas.  See for example compressibility at

    http://en.wikipedia.org/wiki/Compressibility_factor

Suggested-by: Jukka Lind <jukka.lind@iki.fi>
Signed-off-by: Linus Torvalds <torvalds@linux-foundation.org>
Signed-off-by: Dirk Hohndel <dirk@hohndel.org>

5 files changed, 58 insertions, 36 deletions

diff --git a/dive.c b/dive.c
index 798de357c..2fe7cf32a 100644
--- a/dive.c
+++ b/dive.c
@@ -286,6 +286,29 @@ static void update_min_max_temperatures(struct dive *dive, temperature_t tempera
 	}
 }
 
+/*
+ * At high pressures air becomes less compressible, and
+ * does not follow the ideal gas law any more.
+ *
+ * This tries to correct for that, becoming the same
+ * as to_ATM() at lower pressures.
+ *
+ * THIS IS A ROUGH APPROXIMATION! The real numbers will
+ * depend on the exact gas mix and temperature.
+ */
+double surface_volume_multiplier(pressure_t pressure)
+{
+	double bar = pressure.mbar / 1000.0;
+
+	if (bar > 200)
+		bar = 0.00038*bar*bar + 0.51629*bar + 81.542;
+	return bar_to_atm(bar);
+}
+
+int gas_volume(cylinder_t *cyl, pressure_t p)
+{
+	return cyl->type.size.mliter * surface_volume_multiplier(p);
+}
 
 /*
  * If the cylinder tank pressures are within half a bar
@@ -338,7 +361,7 @@ static void match_standard_cylinder(cylinder_type_t *type)
 		return;
 
 	cuft = ml_to_cuft(type->size.mliter);
-	cuft *= to_ATM(type->workingpressure);
+	cuft *= surface_volume_multiplier(type->workingpressure);
 	psi = to_PSI(type->workingpressure);
 
 	switch (psi) {
@@ -381,7 +404,7 @@ static void match_standard_cylinder(cylinder_type_t *type)
  */
 static void sanitize_cylinder_type(cylinder_type_t *type)
 {
-	double volume_of_air, atm, volume;
+	double volume_of_air, volume;
 
 	/* If we have no working pressure, it had *better* be just a physical size! */
 	if (!type->workingpressure.mbar)
@@ -394,8 +417,8 @@ static void sanitize_cylinder_type(cylinder_type_t *type)
 	if (xml_parsing_units.volume == CUFT) {
 		/* confusing - we don't really start from ml but millicuft !*/
 		volume_of_air = cuft_to_l(type->size.mliter);
-		atm = to_ATM(type->workingpressure);		/* working pressure in atm */
-		volume = volume_of_air / atm;			/* milliliters at 1 atm: "true size" */
+		/* milliliters at 1 atm: "true size" */
+		volume = volume_of_air / surface_volume_multiplier(type->workingpressure);
 		type->size.mliter = volume + 0.5;
 	}
 
diff --git a/dive.h b/dive.h
index 201363153..f6677ab2c 100644
--- a/dive.h
+++ b/dive.h
@@ -203,10 +203,8 @@ static inline double bar_to_atm(double bar)
 	return bar / SURFACE_PRESSURE * 1000;
 }
 
-static inline double to_ATM(pressure_t pressure)
-{
-	return pressure.mbar / (double) SURFACE_PRESSURE;
-}
+/* Volume in mliter of a cylinder at pressure 'p' */
+extern int gas_volume(cylinder_t *cyl, pressure_t p);
 
 static inline int mbar_to_PSI(int mbar)
 {
diff --git a/divelist.c b/divelist.c
index 952b03bb2..1c020afe4 100644
--- a/divelist.c
+++ b/divelist.c
@@ -697,26 +697,19 @@ static int calculate_otu(struct dive *dive)
  */
 static double calculate_airuse(struct dive *dive)
 {
-	double airuse = 0;
+	int airuse = 0;
 	int i;
 
 	for (i = 0; i < MAX_CYLINDERS; i++) {
 		pressure_t start, end;
 		cylinder_t *cyl = dive->cylinder + i;
-		int size = cyl->type.size.mliter;
-		double kilo_atm;
-
-		if (!size)
-			continue;
 
 		start = cyl->start.mbar ? cyl->start : cyl->sample_start;
 		end = cyl->end.mbar ? cyl->end : cyl->sample_end;
-		kilo_atm = (to_ATM(start) - to_ATM(end)) / 1000.0;
 
-		/* Liters of air at 1 atm == milliliters at 1k atm*/
-		airuse += kilo_atm * size;
+		airuse += gas_volume(cyl, start) - gas_volume(cyl, end);
 	}
-	return airuse;
+	return airuse / 1000.0;
 }
 
 /* this only uses the first divecomputer to calculate the SAC rate */
diff --git a/profile.c b/profile.c
index 8b5feae23..d3501d0da 100644
--- a/profile.c
+++ b/profile.c
@@ -1011,32 +1011,40 @@ static void set_sac_color(struct graphics_context *gc, int sac, int avg_sac)
 	}
 }
 
-static double get_local_sac(struct plot_data *entry1, struct plot_data *entry2, struct dive *dive)
+/* Get local sac-rate (in ml/min) between entry1 and entry2 */
+static int get_local_sac(struct plot_data *entry1, struct plot_data *entry2, struct dive *dive)
 {
 	int index = entry1->cylinderindex;
-	int delta_time = entry2->sec - entry1->sec;
-	double depth;
-	long delta_pressure = GET_PRESSURE(entry1) - GET_PRESSURE(entry2);
-	long mliter;
+	cylinder_t *cyl;
+	int duration = entry2->sec - entry1->sec;
+	int depth, airuse;
+	pressure_t a, b;
+	double atm;
 
 	if (entry2->cylinderindex != index)
 		return 0;
-	if (delta_pressure <= 0)
+	if (duration <= 0)
 		return 0;
-	if (delta_time <= 0)
+	a.mbar = GET_PRESSURE(entry1);
+	b.mbar = GET_PRESSURE(entry2);
+	if (!a.mbar || !b.mbar)
 		return 0;
-	depth = (entry1->depth + entry2->depth) / 2.0;
 
-	mliter = dive->cylinder[index].type.size.mliter;
+	/* Mean pressure in ATM */
+	depth = (entry1->depth + entry2->depth) / 2;
+	atm = (double) depth_to_mbar(depth, dive) / SURFACE_PRESSURE;
 
-	return delta_pressure * mliter /
-		(delta_time / 60.0) /
-		depth_to_mbar(depth, dive);
+	cyl = dive->cylinder + index;
+
+	airuse = gas_volume(cyl, a) - gas_volume(cyl, b);
+
+	/* milliliters per minute */
+	return airuse / atm * 60 / duration;
 }
 
 /* calculate the current SAC in ml/min and convert to int */
 #define GET_LOCAL_SAC(_entry1, _entry2, _dive) \
-	(int) get_local_sac(_entry1, _entry2, _dive)
+	get_local_sac(_entry1, _entry2, _dive)
 
 #define SAC_WINDOW 45	/* sliding window in seconds for current SAC calculation */
 
diff --git a/statistics.c b/statistics.c
index 9cbfcc72f..b6b618fbd 100644
--- a/statistics.c
+++ b/statistics.c
@@ -604,10 +604,10 @@ static void show_single_dive_stats(struct dive *dive)
 	/* for the O2/He readings just create a list of them */
 	for (idx = 0; idx < MAX_CYLINDERS; idx++) {
 		cylinder_t *cyl = &dive->cylinder[idx];
-		unsigned int start, end;
+		pressure_t start, end;
 
-		start = cyl->start.mbar ? : cyl->sample_start.mbar;
-		end = cyl->end.mbar ? : cyl->sample_end.mbar;
+		start = cyl->start.mbar ? cyl->start : cyl->sample_start;
+		end = cyl->end.mbar ?cyl->sample_end : cyl->sample_end;
 		if (!cylinder_none(cyl)) {
 			/* 0% O2 strangely means air, so 21% - I don't like that at all */
 			int o2 = cyl->gasmix.o2.permille ? : O2_IN_AIR;
@@ -621,8 +621,8 @@ static void show_single_dive_stats(struct dive *dive)
 		}
 		/* and if we have size, start and end pressure, we can
 		 * calculate the total gas used */
-		if (cyl->type.size.mliter && start && end)
-			gas_used += cyl->type.size.mliter / 1000.0 * (start - end);
+		if (start.mbar && end.mbar)
+			gas_used += gas_volume(cyl, start) - gas_volume(cyl, end);
 	}
 	set_label(single_w.o2he, buf);
 	if (gas_used) {