palinsol (version 0.92)

astro: Compute astronomical parameters in the past or in the future

Description

--

Usage

astro(t,solution,degree=FALSE)
ber78(t,degree=FALSE)
ber90(t,degree=FALSE)
la04(t,degree=FALSE)
precession(t,solution)
obliquity(t,solution,degree=FALSE)

Arguments

t

Time, years after 1950

solution

solution used. One of ber78, ber90 or la04

degree

returns angles in degrees if TRUE

Value

• A vector of 3 (la04) or 4 (ber78 and ber90) astronomical elements ll{ 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 d'Alambert decomposition).

  Angles are returned in radians unless degree=TRUE

Details

Both ber78 and ber90 compute astronomical elements based on a d'Alembert 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.

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 valide for about 50 Myr around the present.

precession and obliquity do as astro, but only return precession (e sin varpi) parameter and obliquity, respectively.

References

Berger, A. L. (1978). Long-term variations of daily insolation and Quaternary climatic changes, J. Atmos. Sci., 35, 2362-2367.

Berger and M.F. Loutre (1991), Insolation values for the climate of the last 10 million years, Quaternary Science Reviews, 10, 297 - 317.

J. Laskar et al. (2004), A long-term numerical solution for the insolation quantities of the Earth, Astron. Astroph., 428, 261-285.

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)