add more d18Osw equations and ensure that they are reversible

R/d18Osw_calc.R

@@ -1,27 +1,50 @@
 #' Calculate the d18Osw from the d18Occ and temperature
 #'
-#' According to Kim & O'neil 1997 as modified by Bemis et al., 1998
+#' Defaults to Kim & O'neil 1997 as modified by Bemis et al., 1998.
 #'
 #' @param temperature  The formation temperature (in °C).
 #' @inheritParams temp_calc
 #' @return The oxygen isotope composition of the sea water in VSMOW.
 #' @author Ilja J. Kocken
+#' @references
+#' Kim, S.-T., & O’Neil, J. R. (1997). Equilibrium and nonequilibrium oxygen isotope effects in synthetic carbonates. Geochimica et Cosmochimica Acta, 61(16), 3461–3475. https://doi.org/10.1016/S0016-7037(97)00169-5
+#'
+#' Bemis, B. E., Spero, H. J., Bijma, J., & Lea, D. W. (1998). Reevaluation of the oxygen isotopic composition of planktonic foraminifera: Experimental results and revised paleotemperature equations. Paleoceanography, 13(2), 150–160. https://doi.org/10.1029/98PA00070
 #' @export
 d18Osw_calc <- function(d18Occ, temperature, equation = NULL) {
   equation <- equation_supported(equation)
 
-  if (equation == "KimONeil1997") {
-    # T (°C) = a + b(δc - δw) + c(δc - δw)^2
-    # with a = 16.1, b = -4.64, and c = 0.09
-    # δw = (sqrt(-4ac + b^2 + 4cT) + b + 2cδc) / 2c
-    d18Osw <- (sqrt(-4 * 16.1 * 0.09 + 4.64^2 + 4 * 0.09 * temperature) + -4.64 + 2 * 0.09 * d18Occ) /
-      (2 * 0.09) +
-      # note the 0.27, which is from conversion from VPDB to VSMOW
-      0.27
+  if (equation == "Shackleton1974") {
+    d18Osw <- general_bemis_d18Osw(d18Occ, temperature, a = 16.9, b = -4.0, c = NA, d = -0.20)
+  } else if (equation == "ErezLuz1983") {
+    d18Osw <- general_bemis_d18Osw(d18Occ, temperature, a = 17.0, b = -4.52, c = 0.03, d = -0.22)
+  } else if (equation == "KimONeil1997") {
+    d18Osw <- general_bemis_d18Osw(d18Occ, temperature, a = 16.1, b = -4.64, c = 0.09, d = -0.27)
+  } else if (equation == "Bemis1998") {
+    # there are a lot of variants in table 1, which one should I choose?
+    # TODO
+    ## d18Osw <- general_bemis_d18Osw(d18Occ, temperature, a = 16.1, b = -4.64, c = 0.09, d = -0.27)
   } else if (equation == "Marchitto2014") {
     # Marchitto et al., 2014 equation 9
     # (δcp - δws + 0.27) = -0.245±0.005t + 0.0011±0.0002t² + 3.58±0.02
     d18Osw <- 0.245 * temperature - 0.0011 * temperature^2 - 3.58 + d18Occ + 0.27
   }
   d18Osw
 }
+
+#' General Bemis et al., 1998 Table 1, inverted
+#'
+#' \eqn{\delta_{w} = \frac{\sqrt{-4ac + b^{2} + 4cT} + b + 2cT}{2c} - d}{δw = (sqrt(-4ac + b^2 + 4cT) + b +2cδc) / 2c - d}
+#'
+#' @param temperature Input temperature.
+#' @inheritParams general_bemis_temp
+#' @returns The calculated \eqn{\delta^{18}O}{δ18O} value in VSMOW.
+#' @references Bemis, B. E., Spero, H. J., Bijma, J., & Lea, D. W. (1998). Reevaluation of the oxygen isotopic composition of planktonic foraminifera: Experimental results and revised paleotemperature equations. Paleoceanography, 13(2), 150–160. https://doi.org/10.1029/98PA00070
+general_bemis_d18Osw <- function(d18Occ, temperature, a, b, c, d) {
+  if (is.na(c)) {
+    d18Osw <- (a + b * (d18Occ - d) - temperature) / b
+  } else {
+    d18Osw <- (sqrt(-4 * a * c + b^2 + 4 * c * temperature) + b + 2 * c * d18Occ) / (2 * c) - d
+  }
+  d18Osw
+}

R/equation_supported.R

@@ -1,7 +1,10 @@
-#' Parse equation name into our internal simplified code.
+#' Check if the provided equation is supported
 #'
+#' Also parses equation name into our internal simplified code.
 #' See details for which equations are currently supported.
 #'
+#' - TODO: add all supported equations and references somewhere.
+#'
 #' - Kim & O'Neil 1997 as recalculated by Bemis et al., 1998:
 #'
 #'   \eqn{T = a + b(\delta_{c} - \delta_{w}) + c(\delta_{c} - \delta_{w})^2}{T = a + b(delta_c - delta_w) + c(delta_c - delta_w)^2},
@@ -16,11 +19,11 @@
 #' @returns Character vector of cleaned up equation reference. Gives an error if it's currently not supported.
 equation_supported <- function(equation) {
   supported_equations <- c(
+    "Shackleton1974",
+    "ErezLuz1983",
     "KimONeil1997", # NOTE: as recalculated in Bemis 1998!
+    ## "Bemis1998",
     "Marchitto2014"
-    # "Bemis1998" TODO
-    # "Erez1983" TODO
-    # "Shackleton1974" TODO
   )
 
   # default
@@ -31,10 +34,17 @@ equation_supported <- function(equation) {
       cli::cli_abort("{.var {equation}} must be a {.type {character}} vector.")
     }
 
-    if (equation %in% c("Marchitto et al., 2014",
-                        "Marchitto2014",
-                        "Marchitto")) {
-      equation <- "Marchitto2014"
+    # allow aliases of how you type the equation reference
+    if (equation %in% c("Shackleton1974",
+                        "Shackleton 1974",
+                        "Shackleton, 1974")) {
+      equation <- "Shackleton1974"
+    }
+    if (equation %in% c("Erez & Luz, 1983",
+                        "Erez & Luz 1983",
+                        "Erez&Luz1983",
+                        "ErezLuz1983")) {
+      equation <- "ErezLuz1983"
     }
     if (equation %in% c("Kim & O'Neil, 1997",
                         "Kim & O'Neil 1997",
@@ -43,6 +53,16 @@ equation_supported <- function(equation) {
                         "KimONeil1997")) {
       equation <- "KimONeil1997"
     }
+    if (equation %in% c("Bemis 1998",
+                        "Bemis 1998",
+                        "Bemis et al., 1998")) {
+      equation <- "Bemis1998"
+    }
+    if (equation %in% c("Marchitto et al., 2014",
+                        "Marchitto2014",
+                        "Marchitto")) {
+      equation <- "Marchitto2014"
+    }
   }
 
   if (!equation %in% supported_equations) {

R/temp_calc.R

@@ -1,21 +1,29 @@
 #' Calculate the temperature from d18Occ and the d18Osw
 #'
-#' This is the relationship from Kim & O'neil 1997,
-#' as updated by Bemis et al., 1998
+#' Defaults to the relationship from Kim & O'neil 1997,
+#' as updated by Bemis et al., 1998.
 #'
 #' @param d18Occ The oxygen isotope composition of the carbonate in VPDB.
 #' @param d18Osw The oxygen isotope composition of the sea water in VSMOW
 #' @param equation Character vector with the equation to use. Defaults to Kim & O'Neil, 1997.
 #' @return The temperature in degrees Celsius.
 #' @author Ilja J. Kocken
 #' @export
+#' @references
+#' Kim, S.-T., & O’Neil, J. R. (1997). Equilibrium and nonequilibrium oxygen isotope effects in synthetic carbonates. Geochimica et Cosmochimica Acta, 61(16), 3461–3475. https://doi.org/10.1016/S0016-7037(97)00169-5
+#'
+#' Bemis, B. E., Spero, H. J., Bijma, J., & Lea, D. W. (1998). Reevaluation of the oxygen isotopic composition of planktonic foraminifera: Experimental results and revised paleotemperature equations. Paleoceanography, 13(2), 150–160. https://doi.org/10.1029/98PA00070
 temp_calc <- function(d18Occ, d18Osw, equation = NULL) {
   equation <- equation_supported(equation)
-  if (equation == "KimONeil1997") {
-    # T (°C) = a + b(δc - δw) + c(δc - δw)^2
-    # with a = 16.1, b = -4.64, and c = 0.09
-    d18Osw <- d18Osw - 0.27 # convert from VSMOW to VPDB
-    temp <- 16.1 - 4.64 * (d18Occ - d18Osw) + 0.09 * (d18Occ - d18Osw)^2
+  if (equation == "Shackleton1974") {
+    temp <- general_bemis_temp(d18Occ, d18Osw, a = 16.9, b = -4.0, c = NA, d = -0.20)
+  } else if (equation == "ErezLuz1983") {
+    temp <- general_bemis_temp(d18Occ, d18Osw, a = 17.0, b = -4.52, c = 0.03, d = -0.22)
+  } else if (equation == "KimONeil1997") {
+    temp <- general_bemis_temp(d18Occ, d18Osw, a = 16.1, b = -4.64, c = 0.09, d = -0.27)
+  } else if (equation == "Bemis1998") {
+    # there are a lot of variants in table 1
+    ## d18Osw <- general_bemis_temp(d18Occ, temperature, a = 17.0, b = -4.52, c = 0, d = -0.20)
   } else if (equation == "Marchitto2014") {
     # Equation 9
     # (δcp - δws + 0.27) = -0.245±0.005t + 0.0011±0.0002t² + 3.58±0.02
@@ -25,3 +33,25 @@ temp_calc <- function(d18Occ, d18Osw, equation = NULL) {
   }
   temp
 }
+
+#' General Bemis et al. 1998 table 1
+#'
+#' \eqn{T (°C) = a + b(\delta_{c} - \delta_{w}) + c(\delta_{c}- \delta_{w})^{2}}{T (°C) = a + b(δc - δw) + c(δc - δw)^2}
+#'
+#' @param d18Occ \eqn{\delta^{18}O}{δ18O} of the calcite in VPDB.
+#' @param d18Osw \eqn{\delta^{18}O}{δ18O} of the sea water in VSMOW.
+#' @param a Intercept.
+#' @param b Slope.
+#' @param c Slope of quadratic term.
+#' @param d The conversion correction from VSMOW to VPDB.
+#' @returns The calculated temperature in °C.
+#'
+#' @references Bemis, B. E., Spero, H. J., Bijma, J., & Lea, D. W. (1998). Reevaluation of the oxygen isotopic composition of planktonic foraminifera: Experimental results and revised paleotemperature equations. Paleoceanography, 13(2), 150–160. https://doi.org/10.1029/98PA00070
+general_bemis_temp <- function(d18Occ, d18Osw, a, b, c, d) {
+  if (is.na(c)) {
+    temperature <- a + b * (d18Occ - (d18Osw + d))
+  } else {
+    temperature <- a + b * (d18Occ - (d18Osw + d)) + c * (d18Occ - (d18Osw + d))^2
+  }
+  temperature
+}

tests/testthat/test-d18Osw_temp_calculations.R

@@ -1,21 +1,53 @@
+test_that("All equations work", {
+  expect_error(d18Osw_calc(5, 3, equation = "hoi"))
+  expect_error(temp_calc(5, 3, equation = "hoi"))
+  temp_calc(d18Occ = 2, d18Osw = 0, equation = "Shackleton1974")
+  temp_calc(d18Occ = 2, d18Osw = 0, equation = "ErezLuz1983")
+  temp_calc(d18Occ = 2, d18Osw = 0, equation = "KimONeil1997")
+  temp_calc(d18Occ = 2, d18Osw = 0, equation = "Marchitto2014")
+  d18Osw_calc(d18Occ = 2, temperature = 5, equation = "Shackleton1974")
+  d18Osw_calc(d18Occ = 2, temperature = 5, equation = "ErezLuz1983")
+  d18Osw_calc(d18Occ = 2, temperature = 5, equation = "KimONeil1997")
+  d18Osw_calc(d18Occ = 2, temperature = 5, equation = "Marchitto2014")
+})
+
 test_that("They're reversible", {
+  # test Shackleton, 1974
+  expect_equal(d18Osw_calc(d18Occ = 3,
+                           temp = temp_calc(d18Occ = 3,
+                                            d18Osw = 0,
+                                            equation = "Shackleton1974"),
+                           equation = "Shackleton1974"),
+               expected = 0)
+  expect_equal(temp_calc(d18Occ = 3,
+                         d18Osw = d18Osw_calc(d18Occ = 3, temp = 5,
+                                              equation = "Shackleton1974"),
+                         equation = "Shackleton1974"),
+               expected = 5)
+
+  # test Erez & Luz, 1983
+  expect_equal(d18Osw_calc(d18Occ = 3,
+                           temp = temp_calc(d18Occ = 3,
+                                            d18Osw = 0,
+                                            equation = "ErezLuz1983"),
+                           equation = "ErezLuz1983"),
+               expected = 0)
+  expect_equal(temp_calc(d18Occ = 3,
+                         d18Osw = d18Osw_calc(d18Occ = 3, temp = 5,
+                                              equation = "ErezLuz1983"),
+                         equation = "ErezLuz1983"),
+               expected = 5)
+
   # test Kim & O'Neil, 1997
-  temp_calc(d18Occ = 2, d18Osw = 0)
-  d18Osw_calc(d18Occ = 2, temperature = 5)
-  # is it reversible?
   expect_equal(d18Osw_calc(d18Occ = 3,
                            temp = temp_calc(d18Occ = 3,
                                             d18Osw = 0)),
                expected = 0)
   expect_equal(temp_calc(d18Occ = 3,
                          d18Osw = d18Osw_calc(d18Occ = 3, temp = 5)),
                expected = 5)
-  expect_error(d18Osw_calc(5, 3, equation = "hoi"))
-  expect_error(temp_calc(5, 3, equation = "hoi"))
 
   # test Marchitto et al., 2014
-  temp_calc(d18Occ = 2, d18Osw = 0, equation = "Marchitto2014")
-  d18Osw_calc(d18Occ = 2, temperature = 5, equation = "Marchitto2014")
   expect_equal(d18Osw_calc(d18Occ = 3,
                            temp = temp_calc(d18Occ = 3,
                                             d18Osw = 0,