[PATCH 2/2] gas model: replace Redlich-Kwong with least-square quintic

[Linus Torvalds]

From: Linus Torvalds <torvalds at linux-foundation.org>
Date: Wed, 2 Mar 2016 14:11:53 -0800
Subject: [PATCH 2/2] gas model: replace Redlich-Kwong with least-square quintic

This goes back to just doing air compressibility, but using the
least-squares quintic polynomial equation that Lubomir generated based
on the Wikipedia table for air at 300K in the 1-500 bar range.

We might be able to do similar things for mixed gases..

Signed-off-by: Linus Torvalds <torvalds at linux-foundation.org>
---

I "verified" Lubomir's quintic by checking that yes, the AL80 now is 
reported as "80(77)cuft" again. Other than that I just took his values on 
faith.

Lubomir, what tool did you use for the least-square fitting?

Because it would be interesting to see if we can do this for the table at

  http://www.baue.org/library/zfactor_table.php

and see how reasonable it is to approximate gas mixes by doing the 
O2/N2/He calculations and then just "mixing" the compressibility (looking 
at the O2 and N2 tables and then the air table, I think a stupid and 
straightforward linear mixing is probably _entirely_ unphysical, but not 
necessarily too far off).

For example, looking at 3000 psi, and then just doing

  0.947*0.21 + 1.054*0.79

to do a linear mixing of the O2/N2 values gives you 1.0315, which is quite 
close to the actual air compressibility value.

So I think it might be possible to take those baue.org tables, create 
polynomials for them, and then have a simple linear mix for the gases, and 
get reaonable approximations for all relevant scuba gases..

Am I crazy? Maybe.

 subsurface-core/gas-model.c | 67 +++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 65 insertions(+), 2 deletions(-)

diff --git a/subsurface-core/gas-model.c b/subsurface-core/gas-model.c
index 9cffd5e470e0..c9a6a323955e 100644
--- a/subsurface-core/gas-model.c
+++ b/subsurface-core/gas-model.c
@@ -13,7 +13,7 @@
  *
  */
 
-double redlich_kwong_equation(double t_red, double p_red, double z_init)
+static double redlich_kwong_equation(double t_red, double p_red, double z_init)
 {
 	return (1.0/(1.0 - 0.08664*p_red/(t_red * z_init)) -
 		 0.42748/(sqrt(t_red * t_red * t_red) * ((t_red*z_init/p_red + 0.08664))));
@@ -25,7 +25,7 @@ double redlich_kwong_equation(double t_red, double p_red, double z_init)
  */
 #define STANDARD_TEMPERATURE 293.0
 
-double gas_compressibility_factor(struct gasmix *gas, double bar)
+static double redlich_kwong_compressibility_factor(struct gasmix *gas, double bar)
 {
 	/* Critical points according to https://en.wikipedia.org/wiki/Critical_point_(thermodynamics) */
 
@@ -46,3 +46,66 @@ double gas_compressibility_factor(struct gasmix *gas, double bar)
 				       redlich_kwong_equation(STANDARD_TEMPERATURE/tc, bar/pc,
 							      redlich_kwong_equation(STANDARD_TEMPERATURE/tc, bar/pc,1.0))));
 }
+
+/*
+ * This is a quintic formula by Lubomir I. Ivanov that has
+ * been optimized for the least-square error to the air
+ * compressibility factor table (at 300K) taken from Wikipedia:
+ *
+ * bar  z_factor
+ * ---  ------
+ *   1: 0.9999
+ *   5: 0.9987
+ *  10: 0.9974
+ *  20: 0.9950
+ *  40: 0.9917
+ *  60: 0.9901
+ *  80: 0.9903
+ * 100: 0.9930
+ * 150: 1.0074
+ * 200: 1.0326
+ * 250: 1.0669
+ * 300: 1.1089
+ * 400: 1.2073
+ * 500: 1.3163
+ */
+static double air_compressibility_factor(double bar)
+{
+	double	x0 = 1.0,
+		x1 = bar,
+		x2 = x1*x1,
+		x3 = x2*x1,
+		x4 = x2*x2,
+		x5 = x2*x3;
+
+	return  + x0 * 1.0002556612420115
+		- x1 * 0.0003115084635183305
+		+ x2 * 0.00000227808965401253
+		+ x3 * 1.91596422989e-9
+		- x4 * 8.78421542e-12
+		+ x5 * 6.77746e-15;
+}
+
+/*
+ * We end up using specialized functions for known gases, because
+ * we have special tables for them.
+ *
+ * For now, let's do just air.
+ *
+ * We have other tables for other gases, see for example:
+ *
+ *   http://ww.baue.org/library/zfactor_table.php
+ *
+ * and then we have the Redlich-Kwong function, but that seems
+ * to be almost too generic, and not specific enough to the very
+ * particular pressure and temperature ranges we care about..
+ */
+double gas_compressibility_factor(struct gasmix *gas, double bar)
+{
+#if 1
+	return air_compressibility_factor(bar);
+#else
+	/* Fall back on generic function */
+	return redlich_kwong_compressibility_factor(gas, bar);
+#endif
+}
-- 
2.7.2.334.g7c0da37