gas model: use virial cubic polynomial form

The "virial" form of the Z compression factor is of the form

   Z = 1.0 + A*p + B*p^2 + C*p^3 + ..

and it's considered the "right" polynomial form to use.  It happens to
also make for one constant less per gas (since the 1.0 can be added
later), and can be used to simplify the expression and avoid a few
floating point operations.

However, in order for that kind of expression simplification to make
sense, we need to make sure that we don't calculate the powers of the
pressure multiple times either, and that means we have to inline all the
actual calculations.

Our compiler options still mean that the generated code isn't optimal,
but that's a separate issue. And it is a lot better than it used to be.

Being clever about this does potentially make the code a tiny bit less
legible, but maybe that's not too bad for something that we'd expect to
not ever touch once we get it right.

Signed-off-by: Linus Torvalds <torvalds@linux-foundation.org>
Signed-off-by: Dirk Hohndel <dirk@hohndel.org>

1 files changed, 42 insertions, 56 deletions

diff --git a/subsurface-core/gas-model.c b/subsurface-core/gas-model.c
index 420233612..81765e003 100644
--- a/subsurface-core/gas-model.c
+++ b/subsurface-core/gas-model.c
@@ -4,72 +4,58 @@
 #include <stdlib.h>
 #include "dive.h"
 
-/*
- * Generic cubic polynomial
- */
-static double cubic(double bar, const double coefficient[])
-{
-	double	x0 = 1.0,
-		x1 = bar,
-		x2 = x1*x1,
-		x3 = x2*x1;
-
-	return	x0 * coefficient[0] +
-		x1 * coefficient[1] +
-		x2 * coefficient[2] +
-		x3 * coefficient[3];
-}
+/* "Virial minus one" - the virial cubic form without the initial 1.0 */
+#define virial_m1(C, x1, x2, x3) (C[0]*x1+C[1]*x2+C[2]*x3)
 
 /*
- * Cubic least-square coefficients for O2/N2/He based on data from
+ * Cubic virial least-square coefficients for O2/N2/He based on data from
  *
  *    PERRY’S CHEMICAL ENGINEERS’ HANDBOOK SEVENTH EDITION
  *
  * with the lookup and curve fitting by Lubomir.
  *
+ * The "virial" form of the compression factor polynomial is
+ *
+ *      Z = 1.0 + C[0]*P + C[1]*P^2 + C[2]*P^3 ...
+ *
+ * and these tables do not contain the initial 1.0 term.
+ *
  * NOTE! Helium coefficients are a linear mix operation between the
  * 323K and one for 273K isotherms, to make everything be at 300K.
  */
-static const double o2_coefficients[4] = {
-	+1.00117935180448264158,
-	-0.00074149079841471884,
-	+0.00000291901111247509,
-	-0.00000000162092217461
-};
-
-static const double n2_coefficients[4] = {
-	+1.00030344355797817778,
-	-0.00022528077251905598,
-	+0.00000295430303276288,
-	-0.00000000210649996337
-};
-
-static const double he_coefficients[4] = {
-	+1.00000137322788319261,
-	+0.000488393024887620665,
-	-0.000000054210166760015,
-	+0.000000000010908069275
-};
-
-static double o2_compressibility_factor(double bar) { return cubic(bar, o2_coefficients); }
-static double n2_compressibility_factor(double bar) { return cubic(bar, n2_coefficients); }
-static double he_compressibility_factor(double bar) { return cubic(bar, he_coefficients); }
-
-/*
- * We end up using a simple linear mix of the gas-specific functions.
- */
 double gas_compressibility_factor(struct gasmix *gas, double bar)
 {
-	double o2, n2, he;	// Z factors
-	double of, nf, hf;	// gas fractions ("partial pressures")
-
-	o2 = o2_compressibility_factor(bar);
-	n2 = n2_compressibility_factor(bar);
-	he = he_compressibility_factor(bar);
-
-	of = get_o2(gas) / 1000.0;
-	hf = get_he(gas) / 1000.0;
-	nf = 1.0 - of - nf;
-
-	return o2*of + n2*nf + he*hf;
+	static const double o2_coefficients[3] = {
+		-0.00071809207370164567,
+		+0.00000281852572807643,
+		-0.00000000150290620491
+	};
+	static const double n2_coefficients[3] = {
+		-0.00021926035329221337,
+		+0.00000292844845531647,
+		-0.00000000207613482075
+	};
+	static const double he_coefficients[3] = {
+		+0.00047961098687979363,
+		-0.00000004077670019935,
+		+0.00000000000077707035
+	};
+	double o2, he;
+	double x1, x2, x3;
+	double Z;
+
+	o2 = get_o2(gas) / 1000.0;
+	he = get_he(gas) / 1000.0;
+
+	x1 = bar; x2 = x1*x1; x3 = x2*x1;
+
+	Z = virial_m1(o2_coefficients, x1, x2, x3) * o2 +
+	    virial_m1(he_coefficients, x1, x2, x3) * he +
+	    virial_m1(n2_coefficients, x1, x2, x3) * (1.0 - o2 - he);
+
+	/*
+	 * We add the 1.0 at the very end - the linear mixing of the
+	 * three 1.0 terms is still 1.0 regardless of the gas mix.
+	 */
+	return Z + 1.0;
 }