Insol_l1l2 | R Documentation |
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
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, ...)
orbit |
Output from a solution, such as |
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 |
d1 |
lower calendar day (360-day year) of the time-integral |
d2 |
upper calendar day (360-day year) of the time-integral |
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
Time-integrated insolation in kJ/m2 if avg=TRUE
, else
time-average in W/m2
Insol_d1d2()
: Mean and integrated insolation over an interval bounded by calendar days
Michel Crucifix, U. catholique de Louvain, Belgium.
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
## 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)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.