[PATCH] gas-model: add R compressibility script

[Linus Torvalds]

On Wed, Mar 9, 2016 at 2:07 PM, Lubomir I. Ivanov <neolit123 at gmail.com>
wrote:
>
> the source which wasn't consistent was the helium experimental data:
> https://shareok.org/bitstream/handle/11244/2062/6614196.PDF?sequence=1
>
> (in-document page 108).

Ok. So what I would suggest doing for that data is to actually do a
least-square fit of *both* of those temperatures together.

That should do a fit right in between. As it indeed does:


​which seems to be a simple and reasonable approach to get a line that goes
in the middle of the 273K and 323K data.

It goes all the way up to 700 bar, because that's what the source does.

With R, that gets me the virial coefficients of

  p1  4.87320026468e-04
  p2 -8.83632921053e-08
  p3  5.33304543646e-11

and the plot doesn't look that far off what we have now.

Patch against the previous R script attached, if people want to verify
this. Because it *would* be good to have somebody sanity-check my
calculations..

            Linus
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.subsurface-divelog.org/pipermail/subsurface/attachments/20160309/30d7f861/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: helium-nls.png
Type: image/png
Size: 21732 bytes
Desc: not available
URL: <http://lists.subsurface-divelog.org/pipermail/subsurface/attachments/20160309/30d7f861/attachment-0001.png>
-------------- next part --------------
 subsurface-core/compressibility.r | 33 ++++++++++++++++++++++++++++++++-
 1 file changed, 32 insertions(+), 1 deletion(-)

diff --git a/subsurface-core/compressibility.r b/subsurface-core/compressibility.r
index 15aca6ea27c9..05784dd3a790 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,35 @@ 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)
+
+new = data.frame(x2 = seq(min(x2),max(x2),len=200))
+lines(new$x2,predict(hefit2,newdata=new))