palinsol (version 1.0)

Insol_l1l2: Time-integrated insolation

Description

Computes time-integrated incoming solar radiation (Insol) either between given true solar longitudes (Insol_l1l2) or days of year (Insol_d1d2) for a given orbit and latitude

Usage

Insol_l1l2(
  orbit,
  l1 = 0,
  l2 = 2 * pi,
  lat = 65 * pi/180,
  avg = FALSE,
  ell = TRUE,
  ...
)
Insol_d1d2(orbit, d1, d2, lat = 65 * pi/180, avg = FALSE, ...)

Value

Time-integrated insolation in kJ/m2 if avg=TRUE, else time-average in W/m2

Arguments

orbit

Output from a solution, such as ber78, ber90 or la04

l1

lower true solar longitude bound of the time-integral

l2

upper true solar longitude bound of the time-integral

lat

latitude

avg

performs a time-average.

ell

uses elliptic integrals for the calculation (much faster)

...

other arguments to be passed to Insol

d1

lower calendar day (360-day year) of the time-integral

d2

upper calendar day (360-day year) of the time-integral

Functions

• Insol_d1d2(): Mean and integrated insolation over an interval bounded by calendar days

Author

Michel Crucifix, U. catholique de Louvain, Belgium.

Details

All angles input measured in radians.

Note that in contrast to Berger (2010) we consider the tropic year as the reference, rather than the sideral year, which partly explains some of the small differences with the original publication

References

Berger, A., Loutre, M.F. and Yin Q. (2010), Total irradiation during any time interval of the year using elliptic integrals, Quaternary Science Reviews, 29, 1968 - 1982, doi:10.1016/j.quascirev.2010.05.007

Examples

## reproduces Table 1a of Berger et al. 2010:
lat <- seq(85, -85, -10) * pi/180. ## angles in radians. 
orbit=c(eps=  23.446 * pi/180., ecc= 0.016724, varpi= (102.04 - 180)* pi/180. )
T <-  sapply(lat, function(x) c(lat = x * 180/pi, 
        m1 =  Insol_l1l2(orbit, 0, 70 * pi/180, lat=x, ell= TRUE, S0=1368) / 1e3,
        m2 =  Insol_l1l2(orbit, 0, 70 * pi/180, lat=x, ell=FALSE, S0=1368) / 1e3) ) 
data.frame(t(T))
 ## reproduces Table 1b of Berger et al. 2010:
lat <- c(85, 55, 0, -55, -85) * pi/180. ## angles in radians. 
T <-  sapply(lat, function(x) c(lat = x * 180/pi, 
         m1 =  Insol_l1l2(orbit, 30 * pi/180. , 75 * pi/180, 
               lat=x, ell= TRUE, S0=1368) / 1e3,
         m2 =  Insol_l1l2(orbit, 30 * pi/180. , 75 * pi/180, 
               lat=x, ell=FALSE, S0=1368) / 1e3) ) 
 ## reproduces Table 2a of Berger et al. 2010:
lat <- seq(85, -85, -10) * pi/180. ## angles in radians. 

## 21 march in a 360-d year. By definition : day 80 = 21 march at 12u
d1 = 79.5 
d2 = 79.5 + (10 + 30 + 30 ) * 360/365.2425 ## 30th May in a 360-d year

T <-  sapply(lat, function(x) c(lat = x * 180/pi, 
        m1 =  Insol_d1d2(orbit, d1,d2, lat=x, ell= TRUE, S0=1368) / 1e3,
        m2 =  Insol_d1d2(orbit, d1,d2, lat=x, ell= FALSE, S0=1368) / 1e3))
                          
## I did not quite get the same results as on the table 
## on this one; probably a matter of calendar
## note : the authors in fact used S0=1368 (pers. comm.) 
## 1366 in the paper is a misprint

data.frame(t(T))

## annual mean insolation at 65N North, as a function of the longitude of the perihelion
## (expected to be invariant)

varpis <- seq(0,360,10)*pi/180.  
sapply(varpis, function(varpi)
   {   orbit=c(eps=  23.446 * pi/180., ecc= 0.016724, varpi= varpi )
       amean <- Insol_l1l2 (orbit, lat=65*pi/180., avg=TRUE)
       return(amean)
   })