ep44_output_txt <- "
MD     TVD      pres    temp   gradient
0       0      100.0    92.7     NA
248.7   248.7   127.9   95.3    0.11212
497.4   497.4   156.9   97.9    0.11652
746.2   746.2   186.9   100.4   0.12078
994.9   994.9   218.0   103.0   0.12486
1243.6  1243.6  250.0   105.5   0.12878
1492.3  1492.3  283.0   108.1   0.13256
1741.0  1741.0  316.8   110.7   0.13622
1989.7  1989.7  351.6   113.2   0.13977
2238.5  2238.5  387.2   115.8   0.14322
2487.2  2487.2  423.7   118.3   0.1466
2735.9  2735.9  461.0   120.9   0.14991
2984.6  2984.6  499.6   123.4   0.15525
3233.3  3233.3  540.3   126.0   0.1637
3482.1  3482.1  583.2   128.5   0.17228
3730.8  3730.8  628.2   131.1   0.18095
3979.5  3979.5  675.3   133.6   0.18966
4228.2  4228.2  724.7   136.2   0.19836
4476.9  4476.9  776.2   138.7   0.20701
4725.6  4725.6  829.8   141.2   0.21555
4974.4  4974.4  885.5   143.7   0.22393
5223.1  5223.1  943.2   146.2   0.23212
5471.8  5471.8  1002.9  148.7   0.24006
5720.5  5720.5  1064.5  151.2   0.24772
5969.2  5969.2  1128.0  153.6   0.25508
6217.9  6217.9  1193.1  156.1   0.26211
6466.7  6466.7  1260.0  158.5   0.26878
6715.4  6715.4  1328.4  160.8   0.27508
6964.1  6964.1  1398.3  163.1   0.28101
7212.8  7212.8  1471.2  165.4   0.29291
7461.5  7461.5  1545.0  167.6   0.29689
7710.3  7710.3  1619.7  169.7   0.30049
7959.0  7959.0  1695.3  171.7   0.30373
8207.7  8207.7  1771.5  173.5   0.30665
8456.4  8456.4  1848.4  175.3   0.30889
8705.1  8705.1  1925.2  176.8   0.30894
8953.8  8953.8  2002.1  178.1   0.30901
9202.6  9202.6  2078.9  179.1   0.30912
9451.3  9451.3  2155.9  179.8   0.30926
9700.0  9700.0  2232.8  180.0   0.30947
"
library(readr)
ep44_output <- readr::read_table2(ep44_output_txt, col_names = TRUE)
ep44_output
# gradient as function of depth (TVD)
library(ggplot2)

ggplot(ep44_output, aes(x=gradient, y=pres)) +
    scale_y_reverse(limits = c(max(ep44_output$pres), 0), 
                    breaks = seq(0, max(ep44_output$pres), 500)) +
    geom_line() +
    geom_point()

Polynomial Regression

# third order polynomial
library(dplyr)
# not NAs allowed
ep44_output_nna <- ep44_output %>% 
    na.omit() %>% 
    print()


x <- ep44_output_nna$pres
y <- ep44_output_nna$gradient

model_3 <- lm(y ~ poly(x, 3, raw = TRUE))
model_3
# function of average pressure
p <- function(x) 9.364e-02 + 1.214e-04 * x +  4.784e-08  * x^2 -2.695e-11 * x^3
# function of depth
# full coefficients
t <- function(x) 1.202e-01 -6.175e-06 *x + 8.329e-09 *x*x -5.871e-13 *x*x*x

# rounded coefficients
u <- function(x) 1.1e-1 - 6.2e-6 *x + 8.3e-9 * x^2 - 6.5e-13 * x^3
model_3$fitted.values
plot(x, y)
lines(x, model_3$fitted.values, col = "red")
lines(x, p(x), col = "blue")

Sixth degree

model <- lm(y ~ poly(x, 6, raw = TRUE))
model
plot(model$fitted.values, x)
lines(y, x, col = "red")
fx <- function(x) {
    1.037e-01 + 3.620e-05 * x - 2.208e-08 * x^2 + 8.760e-12 * x^3 - 
        1.406e-15 * x^4 + 1.025e-19 * x^5 - 2.900e-24 * x^6 }



fx(ep44_output_nna$TVD)
gx <- function(x) {
    1.04e-1 + 3.6e-5 * x - 2.21e-8 * x^2 + 8.7e-12 * x^3 -
        1.41e-15 * x^4 + 1.03e-19 * x^5 - 2.9e-24 * x^6
}

plot(y, x)
lines(model$fitted.values, x, col = "red")
lines(fx(x), x, col = "blue")
# lines(gx(x), x, col = "green")
# function of depth
fD  <- function(x) 1.037e-01 + 3.620e-05 * x - 2.208e-08 * x^2 + 
    8.760e-12 * x^3 - 1.406e-15 * x^4 + 1.025e-19 * x^5 - 2.900e-24 * x^6 

# function of pressure
fP <- function(x) 1.1e-1 - 6.2e-6 *x + 8.3e-9 * x^2 - 6.5e-13 * x^3


f0nzie/rNodal documentation built on May 6, 2019, 10:14 a.m.