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.
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"
)snvec() returns different output depending on the outputs argument.
If output = "nice" (the default), returns a
tibble with the following columns:
time Time in thousands of years (kyr).
epl Calculated Obliquity \(\epsilon\) (radians).
phi Calculated Precession \(\phi\) (radians) from ECLIPJ2000.
lpx Calculated Longitude of Perihelion with respect to the moving node \(\bar{\omega}\).
cp Calculated Climatic precession (-) as \(e\sin\bar{\omega}\).
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.
sx, sy, sz The \(x\), \(y\), and \(z\)-components of Earth's
spin axis unit vector \(\vec{s}\) in the heliocentric inertial
reference frame.
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().
Final timestep in thousands of years (kyr).
Defaults to -1000 kyr.
Dynamical ellipticity \(E_{d}\), normalized to modern.
Defaults to 1.0.
Tidal dissipation \(T_{d}\), normalized to modern. Defaults
to 0.0.
Character vector with the name of the desired
solution. Defaults to "full-ZB18a".
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.
Longitude of ascending node of the solar equator relative to ECLIPJ2000.
Inclination of the solar equator relative to ECLIPJ2000.
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.
Numerical absolute tolerance passed to deSolve::ode()'s
atol. Defaults to 1e-5.
Numerical relative tolerance passed to deSolve::ode()'s
rtol. Defaults to 0.
Character vector specifying the method passed to
deSolve::ode()'s method. Defaults to "vode" for stiff problems with
a variable timestep.
Be quiet?
If TRUE, hide info messages.
If FALSE (the default) print info messages and timing.
Character vector with name of desired output. One of:
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.
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.
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).
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), tools:::Rd_expr_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. tools:::Rd_expr_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. tools:::Rd_expr_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), tools:::Rd_expr_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. tools:::Rd_expr_doi("10.18637/jss.v033.i09").
deSolve::ode() from Soetaert et al., (2010) for the ODE solver that we
use.
get_solution() Get astronomical solutions.
# \donttest{
# \dontshow{
# set the cachedir to a temporary directory
pth <- withr::local_tempdir(pattern = "snvecR")
withr::local_options(snvecR.cachedir = pth)
# }
# default call
snvec(tend = -1e3, ed = 1, td = 0, tres = -0.4)
# }
Run the code above in your browser using DataLab