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>

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 */
+		/* milliliters at 1 atm: "true size" */
-		volume = volume_of_air / atm;			/* milliliters at 1 atm: "true size" */
+		volume = volume_of_air / surface_volume_multiplier(type->workingpressure);
 		type->size.mliter = volume + 0.5;
 	}
 

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)
+/* Volume in mliter of a cylinder at pressure 'p' */
-{
+extern int gas_volume(cylinder_t *cyl, pressure_t p);
-	return pressure.mbar / (double) SURFACE_PRESSURE;
-}
 
 static inline int mbar_to_PSI(int mbar)
 {

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 += gas_volume(cyl, start) - gas_volume(cyl, end);
-		airuse += kilo_atm * size;
 	}
-	return airuse;
+	return airuse / 1000.0;
 }
 
 /* this only uses the first divecomputer to calculate the SAC rate */

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;
+	cylinder_t *cyl;
-	double depth;
+	int duration = entry2->sec - entry1->sec;
-	long delta_pressure = GET_PRESSURE(entry1) - GET_PRESSURE(entry2);
+	int depth, airuse;
-	long mliter;
+	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 /
+	cyl = dive->cylinder + index;
-		(delta_time / 60.0) /
+
-		depth_to_mbar(depth, dive);
+	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 */
 

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;
+		start = cyl->start.mbar ? cyl->start : cyl->sample_start;
-		end = cyl->end.mbar ? : cyl->sample_end.mbar;
+		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)
+		if (start.mbar && end.mbar)
-			gas_used += cyl->type.size.mliter / 1000.0 * (start - end);
+			gas_used += gas_volume(cyl, start) - gas_volume(cyl, end);
 	}
 	set_label(single_w.o2he, buf);
 	if (gas_used) {