Both ber78 and ber90 compute astronomical elements based on a spectral
decomposition (sum of sines and cosines) of obliquity and planetary precession parameters.
ber78 uses the Berger (1978) algorithm and is accurate for +/- 1e6 years about the
present. ber90 uses the Berger and Loutre (1991) algorithm and is accurate
for +/- 3e6 years about the present (but with a tiny accuracy over the last 50 kyr,
usually negligible for any palaeo application, see example below).
la04 interpolates tables provided by Laskar (2004), obtained by a simplectic
numerical integration of the planetary system, in which the Moon is considered as a planet.
This solution is valid for about 50 Myr around the present.
precession and obliquity do as astro, but only return
precession (e sin varpi) parameter and obliquity, respectively.
Value
A vector of 3 (la04) or 4 (ber78 and ber90) astronomical elements
eps
obliquity,
ecc
eccentricity and
varpi true solar longitude of the perihelion.
ber78 and ber90 also return
epsp, the Hilbert transform of obliquity (sines changed
in cosines in the spectral decomposition).
Angles are returned in radians unless degree=TRUE
Author(s)
Michel Crucifix, U. catholique de Louvain, Belgium.
References
Berger, A. L. (1978). Long-term variations of daily
insolation and Quaternary climatic changes, J. Atmos. Sci., 35,
2362-2367, doi:10.1175/1520-0469(1978)035<2362:LTVODI>2.0.CO;2
Berger and M.F. Loutre (1991), Insolation values for the climate of the last 10 million years, Quaternary Science Reviews, 10, 297 - 317, doi:10.1016/0277-3791(91)90033-Q
J. Laskar et al. (2004), A long-term numerical solution for the insolation quantities of the Earth, Astron. Astroph., 428, 261-285, doi:10.1051/0004-6361:20041335
Examples
## compare the obliquity over the last 2 Myr with the three solutions
times <- seq(-2e6,0,1e3)
Obl <- function(t) {c(time=t,ber78=ber78(t)['eps'],
ber90=ber90(t)['eps'], la04=la04(t)['eps'])}
Obls <- data.frame(t(sapply(times,Obl)))
## may take about 10 seconds to run
with(Obls, {
plot(times/1e3, ber78.eps, type='l', xlab='time (kyr)',
ylab='Obliquity (radians)')
lines(times/1e3, ber90.eps, type='l', col='red')
lines(times/1e3, la04.eps, type='l', col='green')
})
legend('topright', c('ber78','ber90','la04'), col=c('black','red','green'), lty=1)
## same but with a zoom over the last 300 000 years:
T <- which (times > -3e5)
with(Obls, {
plot(times[T]/1e3, ber78.eps[T], type='l', xlab='time (kyr)',
ylab='Obliquity (radians)')
lines(times[T]/1e3, ber90.eps[T], type='l', col='red')
lines(times[T]/1e3, la04.eps[T], type='l', col='green')
})
legend('topright', c('ber78','ber90','la04'), col=c('black','red','green'), lty=1)