[PATCH 1/2] gas model: add proper He compressibility data and do a least-squares fit

[Linus Torvalds]

From: Linus Torvalds <torvalds at linux-foundation.org>
Date: Wed, 9 Mar 2016 16:41:45 -0800
Subject: [PATCH 1/2] gas model: add proper He compressibility data and do a least-squares fit

Lubomir pointed to exactly where he got his data from, so I added that
raw Helium data to the R script, and let the least-squares fit just take
care of the interpolation between 273K and 323K.

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

So doing the least-square fitting to *also* look for a middle ground 
between 273K and 323K may be questionable, but it works, and has the 
advantage of directly using the experimental data that we have rather than 
doing any other intermediate munging of the experimental data.

The He compressibility curve itself doesn't change noticeably. Yes, you 
can see it is two different curves, but just barely. 

 subsurface-core/compressibility.r | 39 ++++++++++++++++++++++++++++++++++++++-
 1 file changed, 38 insertions(+), 1 deletion(-)

diff --git a/subsurface-core/compressibility.r b/subsurface-core/compressibility.r
index 15aca6ea27c9..66310f3aa8d4 100644
--- a/subsurface-core/compressibility.r
+++ b/subsurface-core/compressibility.r
@@ -39,7 +39,6 @@ o2 = c(0.9994, 0.9968, 0.9941, 0.9884, 0.9771, 0.9676, 0.9597, 0.9542, 0.9560, 0
 n2 = c(0.9998, 0.9990, 0.9983, 0.9971, 0.9964, 0.9973, 1.0000, 1.0052, 1.0559, 1.1422, 1.2480, 1.3629)
 he = c(1.0005, 1.0024, 1.0048, 1.0096, 1.0191, 1.0286, 1.0381, 1.0476, 1.0943, 1.1402, 1.1854, 1.2297)
 
-
 options(digits=15)
 
 #
@@ -76,3 +75,41 @@ summary(hefit)
 
 new = data.frame(x = seq(min(x),max(x),len=200))
 lines(new$x,predict(hefit,newdata=new))
+
+#
+# Raw data from VOLUMETRIC BEHAVIOR OF HELIUM-ARGON MIXTURES [..]
+# T=323.15K (50 C)
+p323atm = c(674.837, 393.223, 237.310, 146.294, 91.4027, 57.5799, 36.4620, 23.1654, 14.7478, 9.4017, 5.9987, 3.8300,
+            540.204, 319.943, 195.008, 120.951, 75.8599, 47.9005, 30.3791, 19.3193, 12.3080, 7.8495, 5.0100, 3.1992)
+
+Hez323 = c(1.28067, 1.16782, 1.10289, 1.06407, 1.04028, 1.02548, 1.01617, 1.01029, 1.00656, 1.00418, 1.00267, 1.00171,
+           1.22738, 1.13754, 1.08493, 1.05312, 1.03349, 1.02122, 1.01349, 1.00859, 1.00548, 1.00349, 1.00223, 1.00143)
+
+
+# T=273.15 (0 C)
+p273atm = c(683.599, 391.213, 233.607, 143.091, 89.0521, 55.9640, 35.3851, 22.4593, 14.2908, 9.1072, 5.8095, 3.7083,
+            534.047, 312.144, 188.741, 116.508, 72.8529, 45.9194, 29.0883, 18.4851, 11.7702, 7.5040, 4.7881, 3.0570)
+
+Hez273 = c(1.33969, 1.19985, 1.12121, 1.07494, 1.04689, 1.02957, 1.01874, 1.01191, 1.00758, 1.00484, 1.00309, 1.00197,
+           1.26914, 1.16070, 1.09837, 1.06118, 1.03843, 1.02429, 1.01541, 1.00980, 1.00625, 1.00398, 1.00254, 1.00162)
+
+p323 = p323atm * 1.01325
+p273 = p273atm * 1.01325
+
+x2=append(p323,p273)
+he2=append(Hez323,Hez273)
+
+plot(x2,he2)
+
+hefit2 = nls(he2 ~ 1.0 + p1*x2 + p2*x2^2 + p3*x2^3,
+        start=list(p1=0,p2=0,p3=0))
+summary(hefit2)
+
+he3 = function(bar)
+{
+  1.0 +0.00047961098687979363 * bar -0.00000004077670019935 * bar^2 +0.00000000000077707035 * bar^3
+}
+
+new = data.frame(x2 = seq(min(x2),max(x2),len=200))
+lines(new$x2,predict(hefit2,newdata=new))
+curve(he3, min(x2),max(x2),add=TRUE)
-- 
2.8.0.rc2