gas-model: add R compressibility script

It annoyed me that we hand-waved a bit about how the virial factors were
actually computed in the gas-model.c file, so here's an actual R script
that computes them and plots the results.

You can run it with (for example):

    R --vanilla < compressibility.r

and it will generate a Rplots.pdf of the plots, and the coefficients
will be shown on stdout.

The result actually differs in insignificant ways from the values that
Lubomir computed, which is likely just due to tools.  I used R, Lubomir
seems to have used

    http://polynomialregression.drque.net/online.php

but the actual curve is pretty much the same.

NOTE! R is not entirely happy about the non-linear fit of the Helium
curve: the fit is *so* precise that it failes the R relative-offset
convergence criterion.  That is apparently generally a sign of
artificial data.

That is probably because Lubomir generated them from the linear mix of
two polynomial fits, rather than a linear mix of the original data.  But
maybe the original data was artificial?

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

1 files changed, 78 insertions, 0 deletions

diff --git a/subsurface-core/compressibility.r b/subsurface-core/compressibility.r
new file mode 100644
index 000000000..15aca6ea2
--- /dev/null
+++ b/subsurface-core/compressibility.r
@@ -0,0 +1,78 @@
+# Compressibility data gathered by Lubomir I Ivanov:
+#
+#  "Data obtained by finding two books online:
+#
+#   [1]
+#    PERRY’S CHEMICAL ENGINEERS’ HANDBOOK SEVENTH EDITION
+#    pretty serious book, from which the wiki AIR values come from!
+#
+#     http://www.unhas.ac.id/rhiza/arsip/kuliah/Sistem-dan-Tekn-Kendali-Proses/PDF_Collections/REFERENSI/Perrys_Chemical_Engineering_Handbook.pdf
+#     page 2-165
+#
+#     [*](Computed from pressure-volume-temperature tables in Vasserman monographs)
+#     ^ i have no idea idea what this means, but the values might not be exactly
+#     experimental?!
+#
+#   the only thing this book is missing is helium, thus [2]!
+#
+#   [2]
+#    VOLUMETRIC BEHAVIOR OF HELIUM-ARGON MIXTURES AT HIGH PRESSURE AND MODERATE TEMPERATURE.
+#
+#     https://shareok.org/bitstream/handle/11244/2062/6614196.PDF?sequence=1
+#     page 108
+#
+#
+#     the book has some tables with pressure values in atmosphere units. i'm
+#     converting them bars. one of the relevant tables is for 323K and one for 273K
+#     (both almost equal distance from 300K).
+#
+#     this again is a linear mix operation between isotherms, which is probably not
+#     the most accurate solution but it works.
+#
+# all data sets contain Z values at 300k, while the pressures are in bars in
+# the 1 to 500 range
+#
+#
+
+x  = c(1, 5, 10, 20, 40, 60, 80, 100, 200, 300, 400, 500)
+o2 = c(0.9994, 0.9968, 0.9941, 0.9884, 0.9771, 0.9676, 0.9597, 0.9542, 0.9560, 0.9972, 1.0689, 1.1572)
+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)
+
+#
+# Get the O2 virial coefficients
+#
+plot(x,o2)
+o2fit = nls(o2 ~ 1.0 + p1*x + p2 *x^2 + p3*x^3, start=list(p1=0,p2=0,p3=0))
+summary(o2fit)
+
+new = data.frame(x = seq(min(x),max(x),len=200))
+lines(new$x,predict(o2fit,newdata=new))
+
+#
+# Get the N2 virial coefficients
+#
+plot(x,n2)
+n2fit = nls(n2 ~ 1.0 + p1*x + p2 *x^2 + p3*x^3, start=list(p1=0,p2=0,p3=0))
+summary(n2fit)
+
+new = data.frame(x = seq(min(x),max(x),len=200))
+lines(new$x,predict(n2fit,newdata=new))
+
+#
+# Get the He virial coefficients
+#
+# NOTE! This will not confirm convergence, thus the warnOnly.
+# That may be a sign that the data is possibly artificial.
+#
+plot(x,he)
+hefit = nls(he ~ 1.0 + p1*x + p2 *x^2 + p3*x^3,
+	start=list(p1=0,p2=0,p3=0),
+	control=nls.control(warnOnly=TRUE))
+summary(hefit)
+
+new = data.frame(x = seq(min(x),max(x),len=200))
+lines(new$x,predict(hefit,newdata=new))