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()
# 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")
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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.