Title: Calculate Earth’s Obliquity and Precession in the Past
Version: 3.10.1
Description: Easily calculate precession and obliquity from an orbital solution (defaults to ZB18a from Zeebe and Lourens (2019) <doi:10.1126/science.aax0612>) and assumed or reconstructed values for tidal dissipation (Td) and dynamical ellipticity (Ed). This is a translation and adaptation of the 'C'-code in the supplementary material to Zeebe and Lourens (2022) <doi:10.1029/2021PA004349>, with further details on the methodology described in Zeebe (2022) <doi:10.3847/1538-3881/ac80f8>. The name of the 'C'-routine is 'snvec', which refers to the key units of computation: spin vector s and orbit normal vector n.
License: GPL (≥ 3)
Encoding: UTF-8
RoxygenNote: 7.3.2
Suggests: astrochron, ggplot2, tidyr, testthat, roxygen2, knitr, withr, curl, rmarkdown
Imports: deSolve, cli (≥ 3.4.0), stats, dplyr, tibble, purrr, readr, tidyselect, rlang (≥ 0.4.11), glue, tools, backports (≥ 1.1.6), stringr
Config/testthat/edition: 3
Depends: R (≥ 3.6.0)
VignetteBuilder: knitr
URL: https://japhir.github.io/snvecR/, https://github.com/japhir/snvecR
BugReports: https://github.com/japhir/snvecR/issues
NeedsCompilation: no
Packaged: 2025-03-05 15:07:09 UTC; japhir
Author: Ilja Kocken ORCID iD [aut, cre, trl, cph], Richard Zeebe ORCID iD [aut]
Maintainer: Ilja Kocken <ikocken@hawaii.edu>
Repository: CRAN
Date/Publication: 2025-03-06 12:10:02 UTC

snvecR: Calculate Earth’s Obliquity and Precession in the Past

Description

Easily calculate precession and obliquity from an orbital solution (defaults to ZB18a from Zeebe and Lourens (2019) doi:10.1126/science.aax0612) and assumed or reconstructed values for tidal dissipation (Td) and dynamical ellipticity (Ed). This is a translation and adaptation of the 'C'-code in the supplementary material to Zeebe and Lourens (2022) doi:10.1029/2021PA004349, with further details on the methodology described in Zeebe (2022) doi:10.3847/1538-3881/ac80f8. The name of the 'C'-routine is 'snvec', which refers to the key units of computation: spin vector s and orbit normal vector n.

Author(s)

Maintainer: Ilja Kocken ikocken@hawaii.edu (ORCID) [translator, copyright holder]

Authors:

See Also

Useful links:

Astronomical Solutions PT-ZB18a(x.xxxx,y.yyyy) for the past 100 Myr

Description

Astronomical Solutions PT-ZB18a(x.xxxx,y.yyyy) for the past 100 Myr

Format

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

A data frame with 249,480 rows and 4 columns:

time

Time in thousands of years (kyr).

epl

Obliqity \epsilon (radians).

phi

Axial Precession phi (radians).

cp

Climatic Precession e sin(\bar{\omega}) (unitless).

Source

All astronomical solutions by Zeebe can be found on http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/Astro.html.

They can be loaded into R quickly, using get_solution().

References

Zeebe, R. E. and Lourens, L. J. (2022). Geologically constrained astronomical solutions for the Cenozoic era. Earth and Planetary Science Letters. doi:10.1016/j.epsl.2022.117595

Astronomical Solutions ZB17 for the past 100 Myr

Description

Astronomical Solutions ZB17 for the past 100 Myr

Format

get_solution("ZB17x")

A data frame with 62,501 rows and 3 columns:

time

Time in thousands of years (kyr).

ecc

Eccentricity e (unitless).

inc

Inclination I (degrees).

Source

All astronomical solutions by Zeebe can be found on http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/Astro.html.

They can be loaded into R quickly, using get_solution().

References

Zeebe, R. E. (2017). Numerical Solutions for the orbital motion of the Solar System over the Past 100 Myr: Limits and new results. The Astronomical Journal. doi:10.3847/1538-3881/aa8cce

Astronomical Solution ZB18a for the Past 100 Myr

Description

Astronomical Solution ZB18a for the Past 100 Myr

Format

get_solution("ZB18a-100")

A data frame with 62,501 rows and 3 columns:

time

Time in thousands of years (kyr).

ecc

Eccentricity e (unitless).

inc

Inclination I (degrees).

Source

All astronomical solutions by Zeebe can be found on http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/Astro.html.

They can be loaded into R quickly, using get_solution().

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.'

Astronomical Solution ZB18a for the Past 300 Myr

Description

Astronomical Solution ZB18a for the Past 300 Myr

Format

get_solution("ZB18a-300")

A data frame with 187,501 rows and 3 columns:

time

Time in thousands of years (kyr).

ecc

Eccentricity e (unitless).

inc

Inclination I (degrees).

Source

All astronomical solutions by Zeebe can be found on http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/Astro.html.

They can be loaded into R quickly, using get_solution().

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. and Lourens, L. J. (2022). Geologically constrained astronomical solutions for the Cenozoic era. Earth and Planetary Science Letters. doi:10.1016/j.epsl.2022.117595

Astronomical Solutions ZB20 for the past 300 Myr

Description

Astronomical Solutions ZB20 for the past 300 Myr

Format

get_solution("ZB20x")

A data frame with 187,501 rows and 3 columns:

time

Time in thousands of years (kyr).

ee

Eccentricity e (unitless).

inc

Inclination I (degrees).

Source

All astronomical solutions by Zeebe can be found on http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/Astro.html.

They can be loaded into R quickly, using get_solution().

References

Zeebe, R. E. and Lourens, L. J. (2022). Geologically constrained astronomical solutions for the Cenozoic era. Earth and Planetary Science Letters. doi:10.1016/j.epsl.2022.117595

Astronomical Solutions ZB23.RXX for the past 3.6 Gyr

Description

Astronomical Solutions ZB23.RXX for the past 3.6 Gyr

Format

get_solution("ZB23.Rxx")

A data frame with 8,750,001 rows and 5 columns:

time

Time in thousands of years (kyr).

ecc

Eccentricity e (unitless).

inc

Inclination I (radians).

obliquity

Obliqity \epsilon (radians).

cp

Climatic Precession e sin(\bar{\omega}) (unitless).

Source

All astronomical solutions by Zeebe can be found on http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/Astro.html.

They can be loaded into R quickly, using get_solution().

References

Zeebe, R. E. and Lourens, L. J. (2022). Geologically constrained astronomical solutions for the Cenozoic era. Earth and Planetary Science Letters. doi:10.1016/j.epsl.2022.117595

Full Astronomical Solution ZB18a for the past 100 Myr

Description

The HNBody output of Zeebe & Lourens (2019) after some pre-processing using prepare_solution(). The wikipedia page on Orbital elements describes what the components relate to in order to uniquely specify an orbital plane. The function asks to download the files to the user's cache directory so that they can be accessed more quickly in the future.

Format

get_solution("full-ZB18a")

A data frame with 250,001 rows and 20 columns:

t

Time t (days).

time

Time in thousands of years (kyr).

aa

Semimajor axis a in astronomical units (au).

ee

Eccentricity e (unitless).

inc

Inclination I (degrees).

lph

Longitude of perihelion \varpi (degrees).

lan

Longitude of the ascending node \Omega (degrees).

arp

Argument of perihelion \omega (degrees).

mna

Mean anomaly M (degrees).

The following columns were computed from the above input with prepare_solution():

lphu

Unwrapped longitude of perihelion \varpi (degrees without jumps).

lanu

Unwrapped longitude of the ascending node \Omega (degrees without jumps).

hh

Variable: e\sin(\varpi).

kk

Variable: e\cos(\varpi).

pp

Variable: 2\sin(0.5I)\sin(\Omega).

qq

Variable: 2\sin(0.5I)\cos(\Omega).

cc

Helper: \cos(I).

dd

Helper: \cos(I)/2.

nnx, nny, nnz

The x, y, and z-components of the Eart's orbit unit normal vector \vec{n}, normal to Earth's instantaneous orbital plane.

Source

All astronomical solutions by Zeebe can be found on http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/Astro.html.

They can be loaded into R quickly, using get_solution().

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. and 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. doi:10.1029/2021PA004349

See Also

prepare_solution()

Get an Astronomical Solution

Description

Download supported astronomical solutions from the web and store it in the user's cache directory. The next use of the function will load the data from the cache rather than downloading it again. This also provides a wrapper for astrochron::getLaskar() if one of their supported solutions is specified, but converts the output to a tibble. Note that we do not cache these solutions locally, however.

Usage

get_solution(
  astronomical_solution = "full-ZB18a",
  quiet = FALSE,
  force = FALSE
)

Arguments

astronomical_solution

Character vector with the name of the desired solution. Defaults to "full-ZB18a".

quiet

Be quiet?

  • If TRUE, hide info messages.

  • If FALSE (the default) print info messages and timing.

force

Force re-downloading the results, even if the solution is saved to the cache.

Value

A tibble with the astronomical solution (and some preprocessed new columns).

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. and 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. doi:10.1029/2021PA004349

See Also

full_ZB18a, ZB17, ZB18a_100, ZB18a_300 ZB20, PT_ZB18a, ZB23

Examples



get_solution("full-ZB18a")    # input for snvec
get_solution("ZB18a-300")     # eccentricity
get_solution("ZB20a")
get_solution("La11")
get_solution("PT-ZB18a(1,1)") # pre-computed precession-tilt
get_solution("ZB23.R01")      # one of the 3.6 Gyr solutions

Prepare Astronomical Solution

Description

Calculates helper columns from an astronomical solution input.

Usage

prepare_solution(data, quiet = FALSE)

Arguments

data

A data frame with the following columns:

  • t Time t (days).

  • ee Eccentricity e (unitless).

  • lph Longitude of perihelion \varpi (degrees).

  • lan Longitude of the ascending node \Omega (degrees).

  • inc Inclination I (degrees).

The easiest way to get this is with get_solution().

quiet

Be quiet?

  • If TRUE, hide info messages.

  • If FALSE (the default) print info messages and timing.

Details

New columns include:

Value

A tibble with the new columns added.

See Also

get_solution()

Calculate Earth’s Obliquity and Precession in the Past

Description

snvec() computes climatic precession and obliquity (or tilt) from an astronomical solution (AS) input and input values for dynamical ellipticity (E_{d}) and tidal dissipation (T_{d}). It solves a set of ordinary differential equations.

Usage

snvec(
  tend = -1000,
  ed = 1,
  td = 0,
  astronomical_solution = "full-ZB18a",
  os_ref_frame = "HCI",
  os_omt = NULL,
  os_inct = NULL,
  tres = -0.4,
  atol = 1e-05,
  rtol = 0,
  solver = "vode",
  quiet = FALSE,
  output = "nice"
)

Arguments

tend

Final timestep in thousands of years (kyr). Defaults to -1000 kyr.

ed

Dynamical ellipticity E_{d}, normalized to modern. Defaults to 1.0.

td

Tidal dissipation T_{d}, normalized to modern. Defaults to 0.0.

astronomical_solution

Character vector with the name of the desired solution. Defaults to "full-ZB18a".

os_ref_frame

Character vector with the reference frame of the astronomical solution. Either "HCI" for heliocentric inertial reference frame or "J2000" for ecliptic J2000 reference frame. Defaults to "HCI" for HNBody output.

os_omt

Longitude of ascending node of the solar equator relative to ECLIPJ2000.

os_inct

Inclination of the solar equator relative to ECLIPJ2000.

tres

Output timestep resolution in thousands of years (kyr). Defaults to -0.4. To determine the sign, think of from 0 to tend by timestep tres.

atol

Numerical absolute tolerance passed to deSolve::ode()'s atol. Defaults to 1e-5.

rtol

Numerical relative tolerance passed to deSolve::ode()'s rtol. Defaults to 0.

solver

Character vector specifying the method passed to deSolve::ode()'s method. Defaults to "vode" for stiff problems with a variable timestep.

quiet

Be quiet?

  • If TRUE, hide info messages.

  • If FALSE (the default) print info messages and timing.

output

Character vector with name of desired output. One of:

  • "nice" (the default) A tibble with the columns time, eei, epl, phi, cp.

  • "full" A tibble with all the computed and interpolated columns.

  • "ode" A matrix with the output of the ODE solver.

Details

This is a re-implementation of the C-code in the supplementary information of Zeebe & Lourens (2022). The terms are explained in detail in Zeebe (2022).

Value

snvec() returns different output depending on the outputs argument.

If output = "nice" (the default), returns a tibble with the following columns:

where \bar{\omega} is the longitude of perihelion relative to the moving equinox.

If output = "all" (for developers), additional columns are included, typically interpolated to output timescale.

See the source code for descriptions of all the intermediate computational steps.

If output = "ode", it will return the raw output of the ODE solver, which is an object of class deSolve and matrix, with columns time, sx, sy, and sz. This can be useful for i.e. deSolve::diagnostics().

Reference Frames of Astronomical Solutions

NASA provides their asteroid and planet positions in the ecliptic J2000 reference frame, while long-term astronomical solution integrations are often performed in the heliocentric inertial reference frame (HCI) or in the inertial reference frame. This is to align the reference frame with the spin vector of the Sun, making J2 corrections intuitive to implement.

Obliquity is typically given in the ecliptic reference frame, so snvec converts all outputs to J2000 if the os_ref_frame is equal to "HCI" and does no transformations if it is already in "J2000".

For this, it uses \Omega_{\odot} = 75.5940 and i_{\odot} = 7.155 as in Zeebe (2017). You can overwrite these defaults with os_omt and os_inct if desired.

ODE Solver

Note that the different ODE solver algorithm we use (Soetaert et al., 2010) means that the R routine returns an evenly-spaced time grid, whereas the C-routine has a variable time-step. This means we need to explicitly set the stepsize tres.

References

Zeebe, R.E. (2017). Numerical Solutions for the Orbital Motion of the Solar System over the Past 100 Myr: Limits and New Results. The Astronomical Journal, 154(5), doi:10.3847/1538-3881/aa8cce.

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.

Zeebe, R. E. (2022). Reduced Variations in Earth’s and Mars’ Orbital Inclination and Earth’s Obliquity from 58 to 48 Myr ago due to Solar System Chaos. The Astronomical Journal, 164(3), doi:10.3847/1538-3881/ac80f8.

Wikipedia page on Orbital Elements: https://en.wikipedia.org/wiki/Orbital_elements

Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer (2010). Solving Differential Equations in R: Package deSolve. Journal of Statistical Software, 33(9), 1–25. doi:10.18637/jss.v033.i09.

See Also

Examples



# default call
snvec(tend = -1e3, ed = 1, td = 0, tres = -0.4)