[PATCH] Take incompressibility of gas into account at higher pressures

Linus Torvalds torvalds at linux-foundation.org
Mon Feb 25 15:29:39 PST 2013


From: Linus Torvalds <torvalds at linux-foundation.org>
Date: Mon, 25 Feb 2013 15:23:16 -0800
Subject: [PATCH] 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 at iki.fi>
Signed-off-by: Linus Torvalds <torvalds at linux-foundation.org>
---

Ok, this uses the suggested correction for pressures above 200 bar. 

The patch is slightly larger than required, I tried to make sure all the 
uses make sense and matched each other. I found another place 
(get_local_sac()) where we didn't convert bar to atm.

 dive.c       | 31 +++++++++++++++++++++++++++----
 dive.h       |  6 ++----
 divelist.c   | 13 +++----------
 profile.c    | 34 +++++++++++++++++++++-------------
 statistics.c | 10 +++++-----
 5 files changed, 58 insertions(+), 36 deletions(-)

diff --git a/dive.c b/dive.c
index 798de357cf99..2fe7cf32a2d2 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 201363153e65..f6677ab2c04a 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 952b03bb222c..1c020afe4379 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 8b5feae2370f..d3501d0da0b7 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 9cbfcc72f61e..b6b618fbdd0c 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) {
-- 
1.8.2.rc0.16.g20a599e



More information about the subsurface mailing list