Compare R and C implementations
In snvecR: Calculate Earth’s Obliquity and Precession in the Past

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8, fig.height = 4
)
## modern_r <- getRversion() >= "4.1.0"
pth <- withr::local_tempdir(pattern = "snvecR")
withr::local_options(list(snvecR.cachedir = pth))

library(tibble)  # nice dataframes
library(ggplot2) # nice plots
library(snvecR)  # this package

Introduction

The function snvec() uses some of the parameters of a full astronomical solution (AS) such as ZB18a from Zeebe and Lourens (2019) in combination with values for tidal dissipation (T_d) and dynamical ellipticity (E_d) to calculate precession and obliquity (or tilt).

In this vignette we show how we can run snvec() and contrast the result to pre-computed solutions, which were calculated using the c-routine.

Apply `snvec`

For the full 100 Myr available:

dat <- snvec(-1e5, 1, 1, astronomical_solution = "full-ZB18a")

Load PT-solution

pt <- get_solution("PT-ZB18a(1,1)")

Inspect results

We find that despite the different ODE solvers and timesteps, the C- and R-implementations are almost identical up to -60 Myr.

pl <- ggplot(dat, aes(x = time / 1000, y = cp)) +
  labs(x = "Time (Myr)",
       y = "Climatic precession") +
  geom_line(aes(colour = "snvecR ZB18a(1,1)")) +
  geom_line(aes(colour = "snvec  ZB18a(1,1)"),
            data = pt) +
  # add eccentricity
  geom_line(aes(y = ee, colour = "ZB18a eccentricity"),
            linetype = "solid",
            data = get_solution("full-ZB18a")) +
  labs(colour = "")
pl + xlim(-60, -59)

plo <- ggplot(dat, aes(x = time / 1000, y = epl)) +
  labs(x = "Time (Myr)",
       y = "Obliquity (rad)") +
  geom_line(aes(colour = "snvecR ZB18a(1,1)")) +
  geom_line(aes(colour = "snvec  ZB18a(1,1)"), data = pt) +
  labs(colour = "")
plo + xlim(-60, -59)

But note the subtle differences at around -100 Myr. This is not significant, however, because this difference occurs far beyond the horizon of predictability in the orbital solutions (the eccentricity curves).

pl + xlim(-100, -99)

plo + xlim(-100, -99)

References

Zeebe, R. E., & Lourens, L. J. (2019). Solar System chaos and the Paleocene–Eocene boundary age constrained by geology and astronomy. Science, 365(6456), 926–929. doi:10.1126/science.aax0612.

Zeebe, R. E., & Lourens, L. J. (2022). A deep-time dating tool for paleo-applications utilizing obliquity and precession cycles: The role of dynamical ellipticity and tidal dissipation. Paleoceanography and Paleoclimatology, e2021PA004349. doi:10.1029/2021PA004349.