# Copyright (c) 2012 Michel Crucifix <michel.crucifix@uclouvain.be>
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject
# the following conditions:
# The above copyright notice and this permission notice shall be
# incluudedin all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND INFRINGEMENT
# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# When using this package for actual applications, always
# cite the authors of the original insolation solutions
# Berger, Loutre and/or Laskar, see details in man pages
## The present file is an implementation of
## Andre Berger, Marie-France Loutre, and Qiuzhen Yin,
## Total irradiation during any time interval of the year using elliptic integrals,
## Quaternary Science Reviews, 29, 1968 - 1982 2010
# ------------------------------------------------------------------
# R Code developed for R version 2.15.2 (2012-10-26) -- "Trick or Treat"
# ------------------------------------------------------------------
SIDERAL_YEAR = 365.25636
TROPIC_YEAR = 365.24219876
YEAR = SIDERAL_YEAR
W <- function (phi, eps, ecc, lambda,S0=1365,n=3)
{
# generalisation of W that will cope with negative Lambda or several rotations
nrot <- floor(lambda / (2*pi))
return ( W_(phi,eps,ecc, lambda %% (2*pi), S0, n) +
nrot * W_(phi,eps,ecc, 2*pi, S0, n) )
}
W_ <- function (phi, eps, ecc, lambda,S0=1365,n=3)
{
# require(gsl) ## gnu scientific library ported by Robin Hankin. Thank you Robin !!!
pi2=pi/2
H00 = pi2
seps = sin(eps)
ceps = cos(eps)
sphi = sin(phi)
cphi = cos(phi)
tphi = sphi/cphi
eq34 <- function (lambda)
{
slambda = sin(lambda)
clambda = cos(lambda)
k = seps/cphi
sesl = seps*slambda
tdelta = sesl / sqrt(1-sesl*sesl)
H0 = acos(-tphi*tdelta)
Flk = gsl::ellint_F(lambda,k)
Elk = gsl::ellint_E(lambda,k)
sphi*seps*(H00-clambda*H0) + cphi*Elk +
sphi*tphi*Flk - sphi*tphi*ceps*ceps*
gsl::ellint_P(lambda, k, -seps*seps)
}
#
eq40 <- function(lambda)
{
## the max(-1, min(1,
## is to account for a numerical artefact when lambda = lambda1,2,3,4
slambda = sin(lambda)
clambda = cos(lambda)
k = seps/cphi
k1 = cphi/seps
psi = asin(max(-1,min(1,k*slambda)))
if (clambda < 0) psi = pi-psi
sesl = seps*slambda
tdelta = sesl / sqrt(1-sesl*sesl)
H0 = acos(max(-1,min(1,-tphi*tdelta)))
Fpk = gsl::ellint_F(psi,k1)
Epk = gsl::ellint_E(psi,k1)
Pipk = gsl::ellint_P(psi,k1, -cphi*cphi)
( sphi*seps*(H00-clambda*H0) + seps*Epk +
ceps* ceps/seps * Fpk - sphi*sphi*ceps*ceps/seps*
Pipk)
}
#
eq38 <- function(lambda) { - pi * sphi*seps*cos(lambda) }
T = YEAR * 0.086400 * 1000
xes = sqrt(1-ecc*ecc)
W0 = S0*T/(2*pi*pi*xes)
if (phi >= (pi2-eps) | phi <= -(pi2-eps) )
{
## above polar circle
lambda1 = (asin(cphi/seps) )
lambda2 = pi - lambda1
lambda3 = pi + lambda1
lambda4 = 2*pi - lambda1
WW=0
## THERE IS A BUG HERE. HAPPY HUNTING !!!
## ref. is p. 1974 of Berger et al. 2010
if (lambda > 0) WW = WW + eq40(min(lambda,lambda1))
if (lambda > lambda2) WW = WW + eq40(min(lambda,lambda3)) - eq40(lambda2)
if (lambda > lambda4) WW = WW + eq40(lambda) - eq40(lambda4)
if (phi >= (pi2-eps) ) { ## northern hemisphere
if (lambda > lambda1) WW = WW + eq38(min(lambda,lambda2)) - eq38(lambda1)
} else
{ ## Southern hemisphere
if (lambda > lambda3) WW = WW + eq38(min(lambda,lambda4)) - eq38(lambda3)
}
WW = W0*WW
} else ## outside polar circle
{ WW = W0 * eq34(lambda) }
WW
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.