Nothing
R/calcDlySolIrr.R
In macroBiome: A Tool for Mapping the Distribution of the Biomes and Bioclimate
Defines functions
calcDlySolIrr
# ~~~~~~~~~~~~
# description:
# ~~~~~~~~~~~~
# This script contains functions to calculate daily irradiance, i.e.:
# calcDtMOptAM(lat, DCL, Z_0, Z_n, t_1, t_80, argCkd = FALSE)
# calcDlySolIrr(year, n, lat, elv = NULL, aprchI_0 = c("Yin1999", "Iqbal1983"), aprchDCL = c("Brock1981", "CBM"),
# argCkd = FALSE)
#
#### DEFINE FUNCTIONS ################################################################################################
# ********************************************************************************************************************
# Name: calcDlySolIrr
# Inputs: - double, year (using astronomical year numbering) (year)
# - double, day of the year (n)
# - double, latitude, deg (lat)
# - double, elevation, m (elv)
# - character, name of the approach for calculating values of the solar constant (aprchI_0)
# - character, name of the approach for calculating values of the solar declination (aprchDCL)
# - logical, indicates whether or not the checking and correction of arguments should be omitted (argCkd)
# Returns: - list object (l.SolIrr)
# $I_0_MJ.m2.hr1 ... solar constant, MJ m-2 hr-1
# $DCL_deg ......... solar declination, deg
# $Z_0_deg ......... solar zenith angle at solar noon, deg
# $Z_n_deg ......... solar zenith angle at sunset (or at solar midnight when the sun does not rise), deg
# $t_1_hr .......... sunrise/sunset hour (solar hour at sunrise and sunset), hr
# $t_80_hr ......... solar hour corresponding to the 80° zenith angle, hr
# $m_a ............. daytime mean optical air mass, unitless
# $tau ............. transmission coefficient, unitless
# $f_m ............. an exponent depending on optical air mass, unitless
# $cosZ_a .......... diurnal average of the cosine of the solar zenith angle, deg
# $R_E_MJ.m2.hr1 ... mean hourly solar irradiance under cloudless-sky conditions, MJ m-2 hr-1
# $DL_hr ........... daylength, hr
# Features: This function calculates daily solar irradiance data.
# Depends: - calcDtMOptAM() ... daytime mean optical air mass
# Ref: - Brock TD (1981) Calculating solar radiation for ecological studies. Ecol Model 14(1–2):1-19.
# DOI: 10.1016/0304-3800(81)90011-9
# - Duffie JA, Beckman WA (1980) Solar Engineering of Thermal Processes. Wiley-Interscience, New York, NY
# - Forsythe WC, Rykiel Jr. EJ, Stahl RS, Wu H, Schoolfield RM (1995) A model comparison for daylength as a
# function of latitude and day of year. Ecol Model 80(1):87–95. DOI: 10.1016/0304-3800(94)00034-F
# - Iqbal M (1983) An Introduction to Solar Radiation. Academic Press, London, UK.
# DOI: 10.1016/B978-0-12-373750-2.X5001-0
# - Linacre E (1992) Climate, Data and Resources. A Reference and Guide. Routledge, New York, NY
# - Yin X (1997a) Calculating daytime mean relative air mass. Agric For Meteorol 87(2-3):85-90.
# DOI: 10.1016/S0168-1923(97)00029-4
# - Yin X (1997b) Optical Air Mass: Daily Integration and its Applications. Meteorol Atmos Phys
# 63(3-4):227-233. DOI: 10.1007/BF01027387
# - Yin X (1999) Bright Sunshine Duration in Relation to Precipitation, Air Temperature and Geographic
# Location. Theor Appl Climatol 64(1–2):61–68. DOI: 10.1007/s007040050111
# ********************************************************************************************************************
calcDlySolIrr <- function(year, n, lat, elv = NULL, aprchI_0 = c("Yin1999", "Iqbal1983"),
aprchDCL = c("Brock1981", "CBM"), argCkd = FALSE) {
# ~~~~ FUNCTION WARNINGS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
aprchI_0 <- strex::match_arg(aprchI_0)
aprchDCL <- strex::match_arg(aprchDCL)
err_han <- errorHandling(year = year, n = n, lat = lat, elv = elv, onlyLgthChecking = argCkd)
list2env(Filter(Negate(is.null), err_han), envir = environment())
cv.arg <- c("year", "n", "lat")
for (i in 1 : length(cv.arg)) {
if (is.null(get(cv.arg[i]))) { stop("Invalid argument: '", cv.arg[i], "' is missing, with no default.") }
}
cv.opv <- c("I_0_MJ.m2.hr1", "DCL_deg", "Z_0_deg", "Z_n_deg", "t_1_hr", "t_80_hr", "m_a", "tau", "f_m", "cosZ_a",
"R_E_MJ.m2.hr1", "DL_hr")
if (missing(elv)) {
incl <- which(stats::complete.cases(year = year, n = n, lat = lat))
elv <- NULL
} else {
incl <- which(stats::complete.cases(year = year, n = n, lat = lat, elv = elv))
cv.arg <- c(cv.arg, c("elv"))
}
# ~~~~ FUNCTION VARIABLES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
l.SolIrr <- stats::setNames(lapply(vector("list", length(cv.opv)), function(x) rep(NA, lgth)), cv.opv)
if (identical(integer(0), incl)) {
return(l.SolIrr)
} else {
for (i in 1 : length(cv.arg)) { assign(cv.arg[i], get(cv.arg[i])[incl]) }
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 01. Calculate the solar constant (I_0), MJ m-2 hr-1
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (aprchI_0 == "Iqbal1983") { # ref: Iqbal (1983)
# 'Average' solar constant, I_SC (in MJ m-2 hr-1)
I_SC <- 4.871 # ref: 'text' in Iqbal (1983)
# Correction factor for the eccentricity of Earth’s orbit (unitless)
# ref: Eq 1.2.3 in Iqbal (1983); Eq 1.4.1 in Duffie and Beckman (1980)
E_0 <- 1. + (0.033 * cos(360. * n / 365. * pi /180.))
# Solar constant, I_0 (in MJ m-2 hr-1)
I_0 <- I_SC * E_0 # ref: Eq 4.2.1 in Iqbal (1983)
} else { # ref: Yin (1999); Brock (1981)
# 'Average' solar constant, I_SC (in MJ m-2 hr-1)
I_SC <- 4.9212 # ref: 'text' in Yin (1997b); Yin (1999)
# Radius vector of the Earth, R_1 (unitless)
# ref: Eq 2 in Brock (1981)
R_1 <- 1. / (1. + (0.033 * cos(360. * (n - 1.) / 365. * pi / 180.)))**0.5
# Solar constant, I_0 (in MJ m-2 hr-1)
I_0 <- I_SC /R_1 ** 2. # ref: 'text' in Brock (1981)
}
l.SolIrr$I_0_MJ.m2.hr1[incl] <- I_0
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 02. Calculate the solar declination (DCL), deg
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (aprchDCL == "CBM") { # ref: Eqs 1 and 2 in Forsythe et al. (1995)
DCL <- rad2deg(asin(0.39795 * cos(0.2163108 + 2. * atan(0.9671396 * tan(0.00860 * (n - 186.))))))
} else { # ref: Eq A7 in Yin (1997a); Eq 1 in Brock (1981); Eq 4 in Forsythe et al. (1995)
DCL <- 23.45 * sin(deg2rad(360. / 365. * (n + 283.)))
}
l.SolIrr$DCL_deg[incl] <- DCL
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 03. Calculate the solar zenith angle at solar noon (Z_0), deg
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ref: Eq 2.10 in Yin (1997b)
Z_0 <- ifelse(abs(lat - DCL) < 90., lat - DCL, 90.)
l.SolIrr$Z_0_deg[incl] <- Z_0
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 04. Calculate the solar zenith angle at sunset (or at solar midnight when the sun does not rise) (Z_n), deg
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ref: Eq 2.9 in Yin (1997b)
Z_n <- ifelse(abs(lat + DCL) < 90., 90.,
rad2deg(acos(sin(deg2rad(DCL)) * sin(deg2rad(lat)) - cos(deg2rad(DCL)) * cos(deg2rad(lat)))))
l.SolIrr$Z_n_deg[incl] <- Z_n
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 05. Calculate the daytime mean optical air mass (and the value of an associated experiential function)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Coefficents for calculating optical air mass
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
c_LE80 <- matrix(c(0.008307, 1.021, -0.01259), nrow = 1, dimnames = list("c_LE80", c("c1", "c2", "c3")))
c_GT80 <- matrix(c(0.03716, 1.538, -1.689), nrow = 1, dimnames = list("c_GT80", c("c1", "c2", "c3")))
coef_cho <- rbind(c_LE80, c_GT80)
# Optical air mass at solar noon (corresponding to cosine zenith angle at its maximum), m_o (unitless)
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
coef_val <- matrix(coef_cho[ifelse(Z_0 <= 80., "c_LE80", "c_GT80"), ], nrow = length(incl), ncol = ncol(coef_cho),
dimnames = list(NULL, colnames(coef_cho)))
# ref: Eq 2.1 in Yin (1997b)
m_o <- as.numeric(coef_val[, "c2"] / (coef_val[, "c1"] + cos(deg2rad(Z_0))) + coef_val[, "c3"])
# Optical air mass corresponding to cosine zenith angle at its medium, m_c (unitless)
# ref: definitions for Eq 3.3 in Yin (1997b)
Z_c <- rad2deg(acos((cos(deg2rad(Z_0)) + cos(deg2rad(Z_n))) / 2.))
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
coef_val <- matrix(coef_cho[ifelse(Z_c <= 80., "c_LE80", "c_GT80"), ], nrow = length(incl), ncol = ncol(coef_cho),
dimnames = list(NULL, colnames(coef_cho)))
# ref: Eq 2.1 in Yin (1997b)
m_c <- as.numeric(coef_val[, "c2"] / (coef_val[, "c1"] + cos(deg2rad(Z_c))) + coef_val[, "c3"])
# Optical air mass corresponding to cosine zenith angle at its bottom-quarter range point, m_l (unitless)
# ref: definitions for Eq 3.3 in Yin (1997b)
Z_l <- rad2deg(acos((cos(deg2rad(Z_0)) + 3. * cos(deg2rad(Z_n))) / 4.))
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
coef_val <- matrix(coef_cho[ifelse(Z_l <= 80., "c_LE80", "c_GT80"), ], nrow = length(incl), ncol = ncol(coef_cho),
dimnames = list(NULL, colnames(coef_cho)))
# ref: Eq 2.1 in Yin (1997b)
m_l <- as.numeric(coef_val[, "c2"] / (coef_val[, "c1"] + cos(deg2rad(Z_l))) + coef_val[, "c3"])
# Sunrise/sunset hour (solar hour at sunrise and sunset), t_1 (hr)
# ref: Eq 2.5 in Yin (1997b)
# t_1 <- ifelse(abs(lat - DCL) >= 90., 0.,
# ifelse(abs(lat + DCL) >= 90., 12.,
# (12. / pi) * acos(-1. * tan(deg2rad(lat)) * tan(deg2rad(DCL)))))
t_1 <- vector(length = length(incl))
t_1[abs(lat - DCL) >= 90. & abs(lat + DCL) < 90.] <- 0.
t_1[abs(lat - DCL) < 90. & abs(lat + DCL) >= 90.] <- 12.
sls <- abs(lat - DCL) < 90. & abs(lat + DCL) < 90.
if (!identical(sls, FALSE)) { t_1[sls] <- (12. / pi) * acos(-1. * tan(deg2rad(lat[sls])) * tan(deg2rad(DCL[sls]))) }
l.SolIrr$t_1_hr[incl] <- t_1
# Solar hour corresponding to the 80° zenith angle, t_80 (hr)
# ref: Eq 2.8 in Yin (1997b)
ratio <- (cos(deg2rad(80.)) - (sin(deg2rad(lat)) * sin(deg2rad(DCL)))) / cos(deg2rad(lat)) * cos(deg2rad(DCL))
t_80 <- vector(length = length(incl))
t_80[ratio >= 1.] <- 0.
t_80[ratio <= -1.] <- 12.
t_80[ratio < 1. & ratio > -1.] <- (1. / 15.) * rad2deg(acos(ratio[ratio < 1. & ratio > -1.]))
l.SolIrr$t_80_hr[incl] <- t_80
# Daytime mean optical air mass, m_a (unitless)
# ref: 2.7 Eq in Yin (1997b)
# print(c("lat", "DCL", "Z_0", "Z_n", "t_1", "t_80"))
# print(c(lat, DCL, Z_0, Z_n, t_1, t_80))
m_a <- calcDtMOptAM(lat, DCL, Z_0, Z_n, t_1, t_80, argCkd = TRUE)
# Checking whether use of elevation-correction for the optical air mass values
if (!is.null(elv)) {
# Station atmospheric pressure, sap (unitless)
# ref: Eq 2.11 in Yin (1997b); Linacre (1992)
sap <- exp(-elv / 8000.)
m_o <- m_o * sap
m_c <- m_c * sap
m_l <- m_l * sap
m_a <- m_a * sap
}
l.SolIrr$m_a[incl] <- m_a
# Transmission coefficient, tau (unitless)
# ref: 'text' Yin (1997b), Linacre (1992)
tau <- 0.75
l.SolIrr$tau[incl] <- tau
# An exponent depending on optical air mass, f_m (unitless)
# ref: Eq 3.3 in Yin (1997b)
f_m <- 0.2323 + 0.01452 * (m_a + m_l) * exp(1.403 * tau) - 0.1528 * m_o + m_c + 0.487 * (m_c - m_l)
l.SolIrr$f_m[incl] <- f_m
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 06. Calculate the diurnal average of the cosine of the solar zenith angle
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Half-length of the illuminated day, Hr (in radians)
# ref: Eq A2 in Yin (1997a)
Hr <- vector(length = length(incl))
Hr[abs(lat - DCL) >= 90. & abs(lat + DCL) < 90.] <- 0.
Hr[abs(lat - DCL) < 90. & abs(lat + DCL) >= 90.] <- pi
sls <- abs(lat - DCL) < 90. & abs(lat + DCL) < 90.
if (!identical(sls, FALSE)) { Hr[sls] <- acos(-1. * tan(deg2rad(lat[sls])) * tan(deg2rad(DCL[sls]))) }
# Diurnal average of the cosine of the solar zenith angle, cosZ_a (unitless)
# ref: Eq A5 in Yin (1997a)
cosZ_a <- vector(length = length(incl))
cosZ_a[abs(lat - DCL) >= 90.] <- 0.
sls <- abs(lat - DCL) < 90.
if (!identical(sls, FALSE)) { cosZ_a[sls] <- sin(deg2rad(lat[sls])) * sin(deg2rad(DCL[sls])) +
(1. / Hr[sls]) * cos(deg2rad(lat[sls])) * cos(deg2rad(DCL[sls])) * sin(Hr[sls]) }
l.SolIrr$cosZ_a[incl] <- cosZ_a
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 07. Calculate the mean hourly solar irradiance under cloudless-sky conditions (R_E), MJ m-2 hr-1
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ref: 'text' in Yin (1999); Eq 3.2 in Yin (1997b)
R_E <- I_0 * cosZ_a * (tau ** f_m)
l.SolIrr$R_E_MJ.m2.hr1[incl] <- R_E
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 08. Calculate the daylength (DL), hr
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DL <- vector(length = length(incl))
DL[abs(lat - DCL) >= 90. & abs(lat + DCL) < 90.] <- 0.
DL[abs(lat - DCL) < 90. & abs(lat + DCL) >= 90.] <- 24.
sls <- abs(lat - DCL) < 90. & abs(lat + DCL) < 90.
if (aprchDCL == "CBM") { # ref: Eqs 3 in Forsythe et al. (1995)
DL[sls] <- 24. - (24. / pi) * acos((sin(deg2rad(0.8333)) + sin(deg2rad(lat[sls])) * sin(deg2rad(DCL[sls]))) /
(cos(deg2rad(lat[sls])) * cos(deg2rad(DCL[sls]))))
} else { # ref: 'text' in Yin (1997a); Eq 4 in Brock (1981); Eqs 5 and 6 in Forsythe et al. (1995)
DL[sls] <- 24. / 180. * rad2deg(acos(-1. * tan(deg2rad(lat[sls])) * tan(deg2rad(DCL[sls]))))
}
l.SolIrr$DL_hr[incl] <- DL
# ~~~~ RETURN VALUES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
return(l.SolIrr)
}
# ********************************************************************************************************************
# Name: calcDtMOptAM
# Inputs: - double, latitude, deg (lat)
# - double, solar declination, deg (DCL)
# - double, solar zenith angle at solar noon, deg (Z_0)
# - double, solar zenith angle at sunset (or at solar midnight when the sun does not rise), deg (Z_n)
# - double, sunrise/sunset hour (solar hour at sunrise and sunset), hr (t_1)
# - double, solar hour corresponding to the 80° zenith angle, hr (t_80)
# - logical, indicates whether or not the checking and correction of arguments should be omitted (argCkd)
# Returns: - double, daytime mean optical air mass, unitless (m_a)
# Features: This function returns the daytime mean optical air mass by using the Equation 2.7 in Yin (1997b).
# Ref: - Yin X (1997b) Optical Air Mass: Daily Integration and its Applications. Meteorol Atmos Phys
# 63(3-4):227-233. DOI: 10.1007/BF01027387
# ********************************************************************************************************************
calcDtMOptAM <- function(lat, DCL, Z_0, Z_n, t_1, t_80, argCkd = FALSE) {
# ~~~~ FUNCTION WARNINGS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
cv.arg <- c("lat", "DCL", "Z_0", "Z_n", "t_1", "t_80")
for (i in 1 : length(cv.arg)) {
if (is.null(get(cv.arg[i]))) { stop("Invalid argument: '", cv.arg[i], "' is missing, with no default.") }
}
if (!argCkd) {
maVarFeatures <- data.frame(lgth = rep(1, 6), lwst_val = c(rep(-90, 2), rep(0, 4)),
hgst_val = c(rep(90, 4), rep(24, 2)),
rng_wrd = c(rep("values from the interval [-90, 90]", 2),
rep("values from the interval [0, 90]", 2),
rep("values from the interval [0, 24]", 2)))
rownames(maVarFeatures) <- cv.arg
# Checking the class and content of objects
for (i in 1 : length(cv.arg)) {
arg <- get(cv.arg[i])
if (!is.list(arg) & length(arg) >= 1L &
length(which(!is.na(match(class(arg), c("numeric", "integer"))))) >= 1L) {
if (!(is.numeric(arg) &
all(arg[!is.na(arg)] >= maVarFeatures[cv.arg[i], "lwst_val"] &
arg[!is.na(arg)] <= maVarFeatures[cv.arg[i], "hgst_val"]))) {
stop("Invalid argument: '", cv.arg[i], "' has to be a vector with ", maVarFeatures[cv.arg[i], "rng_wrd"],
".")
}
} else {
stop("Invalid argument: '", cv.arg[i], "' has to be a vector with ", maVarFeatures[cv.arg[i], "rng_wrd"],
".")
}
}
ls.arg <- stats::setNames(list(lat, DCL, Z_0, Z_n, t_1, t_80), cv.arg)
if (!all(sapply(ls.arg, length) == length(ls.arg[[1]]))) {
stop("Invalid argument: The vectors ",
paste0(paste(cv.arg[1 : (length(cv.arg) - 1)], sep = "", collapse = ", "), " and ",
cv.arg[length(cv.arg)]),
" have not the same length.")
}
}
# ~~~~ FUNCTION VARIABLES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
# Coefficents for calculating optical air mass
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
c_LE80 <- matrix(c(0.008307, 1.021, -0.01259), nrow = 1, dimnames = list("c_LE80", c("c1", "c2", "c3")))
c_GT80 <- matrix(c(0.03716, 1.538, -1.689), nrow = 1, dimnames = list("c_GT80", c("c1", "c2", "c3")))
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 01. Select the equation(s) applied
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
eq <- stats::setNames(lapply(vector("list", 4), function(x) vector()), c("a", "b", "c", "d"))
eq$a <- which(t_1 == 0.)
eq$b <- which(t_1 != 0. & Z_n <= 80.)
eq$c <- which(t_1 != 0. & Z_0 >= 80.)
eq$d <- which(t_1 != 0. & Z_n > 80. & Z_0 < 80.)
applEq <- which(lapply(eq, function(x) { !identical(x, integer(0)) } ) == T)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 02. Set the result vector
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
m_a <- rep(NA, length(lat))
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 03. Run the algorithm for each equation type
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for (i_eq in applEq) {
slctd <- eq[[c("a", "b", "c", "d")[i_eq]]]
sn <- 0
while(sn <= 2) {
if (c("a", "b", "c", "d")[i_eq] == "a") {
m_a[slctd] <- 39.7
break
}
if (c("a", "b", "c", "d")[i_eq] == "b") {
t_1g <- t_1[slctd]
cg <- c_LE80
}
if (c("a", "b", "c", "d")[i_eq] == "c") {
t_1g <- t_1[slctd]
cg <- c_GT80
}
if (c("a", "b", "c", "d")[i_eq] == "d") {
sn <- sn + 1
if (sn == 1) t_1g <- t_80[slctd]
if (sn == 2) t_1g <- t_1[slctd]
if (sn == 3) t_1g <- t_80[slctd]
if (sn == 1) cg <- c_LE80
if (sn == 2) cg <- c_GT80
if (sn == 3) cg <- c_GT80
}
a_ar <- cg[1] + sin(deg2rad(lat[slctd])) * sin(deg2rad(DCL[slctd]))
if (sn < 2) {
b_ar <- cos(deg2rad(lat[slctd])) * cos(deg2rad(DCL[slctd]))
}
mbr2 <- mbr3 <- mbr4 <- rep(NA, length(slctd))
sls <- a_ar > b_ar
if (!identical(sls, integer(0))) {
mbr2[sls] <- cg[2] / sqrt(a_ar[sls] ** (2.) - b_ar[sls] ** (2.))
numer <- b_ar[sls] + a_ar[sls] * cos(deg2rad(15. * t_1g[sls]))
denom <- a_ar[sls] + b_ar[sls] * cos(deg2rad(15. * t_1g[sls]))
mbr3[sls] <- acos(numer / denom)
}
sls <- a_ar < b_ar
if (!identical(sls, integer(0))) {
mbr2[sls] <- cg[2] / sqrt(b_ar[sls] ** (2.) - a_ar[sls] ** (2.))
sqr1 <- sqrt((b_ar[sls] + a_ar[sls]) * (1. + cos(deg2rad(15. * t_1g[sls]))))
sqr2 <- sqrt((b_ar[sls] - a_ar[sls]) * (1. - cos(deg2rad(15. * t_1g[sls]))))
mbr3[sls] <- log((sqr1 + sqr2) / (sqr1 - sqr2))
}
sls <- a_ar == b_ar
if (!identical(sls, integer(0))) {
mbr2[sls] <- cg[2] / a_ar[sls]
mbr3[sls] <- tan(deg2rad(15. * t_1g[sls] / 2.))
}
mbr4 <- cg[3] * t_1g
F_eq6 <- 1. / deg2rad(15.) * mbr2 * mbr3 + mbr4
if (sn == 0) { F0_eq6 <- F_eq6 }
if (sn == 1) { F1_eq6 <- F_eq6 }
if (sn == 2) { F2_eq6 <- F_eq6 }
if (sn == 3) { F3_eq6 <- F_eq6 }
if (sn == 1 | sn == 2) { next }
if (sn == 0) {
m_a[slctd] <- 1. / t_1g * F0_eq6
break
}
if (sn == 3) {
m_a[slctd] <- 1. / t_1[slctd] * (F1_eq6 + F2_eq6 - F3_eq6)
break
}
}
}
# ~~~~ RETURN VALUES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
return(m_a)
}
deg2rad <- function(deg) { deg * (pi / 180.) }
rad2deg <- function(rad) { rad * (180. / pi) }
Try the macroBiome package in your browser
Any scripts or data that you put into this service are public.
macroBiome documentation built on May 29, 2024, 5:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions calcDlySolIrr
# ~~~~~~~~~~~~
# description:
# ~~~~~~~~~~~~
# This script contains functions to calculate daily irradiance, i.e.:
# calcDtMOptAM(lat, DCL, Z_0, Z_n, t_1, t_80, argCkd = FALSE)
# calcDlySolIrr(year, n, lat, elv = NULL, aprchI_0 = c("Yin1999", "Iqbal1983"), aprchDCL = c("Brock1981", "CBM"),
# argCkd = FALSE)
#
#### DEFINE FUNCTIONS ################################################################################################
# ********************************************************************************************************************
# Name: calcDlySolIrr
# Inputs: - double, year (using astronomical year numbering) (year)
# - double, day of the year (n)
# - double, latitude, deg (lat)
# - double, elevation, m (elv)
# - character, name of the approach for calculating values of the solar constant (aprchI_0)
# - character, name of the approach for calculating values of the solar declination (aprchDCL)
# - logical, indicates whether or not the checking and correction of arguments should be omitted (argCkd)
# Returns: - list object (l.SolIrr)
# $I_0_MJ.m2.hr1 ... solar constant, MJ m-2 hr-1
# $DCL_deg ......... solar declination, deg
# $Z_0_deg ......... solar zenith angle at solar noon, deg
# $Z_n_deg ......... solar zenith angle at sunset (or at solar midnight when the sun does not rise), deg
# $t_1_hr .......... sunrise/sunset hour (solar hour at sunrise and sunset), hr
# $t_80_hr ......... solar hour corresponding to the 80° zenith angle, hr
# $m_a ............. daytime mean optical air mass, unitless
# $tau ............. transmission coefficient, unitless
# $f_m ............. an exponent depending on optical air mass, unitless
# $cosZ_a .......... diurnal average of the cosine of the solar zenith angle, deg
# $R_E_MJ.m2.hr1 ... mean hourly solar irradiance under cloudless-sky conditions, MJ m-2 hr-1
# $DL_hr ........... daylength, hr
# Features: This function calculates daily solar irradiance data.
# Depends: - calcDtMOptAM() ... daytime mean optical air mass
# Ref: - Brock TD (1981) Calculating solar radiation for ecological studies. Ecol Model 14(1–2):1-19.
# DOI: 10.1016/0304-3800(81)90011-9
# - Duffie JA, Beckman WA (1980) Solar Engineering of Thermal Processes. Wiley-Interscience, New York, NY
# - Forsythe WC, Rykiel Jr. EJ, Stahl RS, Wu H, Schoolfield RM (1995) A model comparison for daylength as a
# function of latitude and day of year. Ecol Model 80(1):87–95. DOI: 10.1016/0304-3800(94)00034-F
# - Iqbal M (1983) An Introduction to Solar Radiation. Academic Press, London, UK.
# DOI: 10.1016/B978-0-12-373750-2.X5001-0
# - Linacre E (1992) Climate, Data and Resources. A Reference and Guide. Routledge, New York, NY
# - Yin X (1997a) Calculating daytime mean relative air mass. Agric For Meteorol 87(2-3):85-90.
# DOI: 10.1016/S0168-1923(97)00029-4
# - Yin X (1997b) Optical Air Mass: Daily Integration and its Applications. Meteorol Atmos Phys
# 63(3-4):227-233. DOI: 10.1007/BF01027387
# - Yin X (1999) Bright Sunshine Duration in Relation to Precipitation, Air Temperature and Geographic
# Location. Theor Appl Climatol 64(1–2):61–68. DOI: 10.1007/s007040050111
# ********************************************************************************************************************
calcDlySolIrr <- function(year, n, lat, elv = NULL, aprchI_0 = c("Yin1999", "Iqbal1983"),
aprchDCL = c("Brock1981", "CBM"), argCkd = FALSE) {
# ~~~~ FUNCTION WARNINGS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
aprchI_0 <- strex::match_arg(aprchI_0)
aprchDCL <- strex::match_arg(aprchDCL)
err_han <- errorHandling(year = year, n = n, lat = lat, elv = elv, onlyLgthChecking = argCkd)
list2env(Filter(Negate(is.null), err_han), envir = environment())
cv.arg <- c("year", "n", "lat")
for (i in 1 : length(cv.arg)) {
if (is.null(get(cv.arg[i]))) { stop("Invalid argument: '", cv.arg[i], "' is missing, with no default.") }
}
cv.opv <- c("I_0_MJ.m2.hr1", "DCL_deg", "Z_0_deg", "Z_n_deg", "t_1_hr", "t_80_hr", "m_a", "tau", "f_m", "cosZ_a",
"R_E_MJ.m2.hr1", "DL_hr")
if (missing(elv)) {
incl <- which(stats::complete.cases(year = year, n = n, lat = lat))
elv <- NULL
} else {
incl <- which(stats::complete.cases(year = year, n = n, lat = lat, elv = elv))
cv.arg <- c(cv.arg, c("elv"))
}
# ~~~~ FUNCTION VARIABLES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
l.SolIrr <- stats::setNames(lapply(vector("list", length(cv.opv)), function(x) rep(NA, lgth)), cv.opv)
if (identical(integer(0), incl)) {
return(l.SolIrr)
} else {
for (i in 1 : length(cv.arg)) { assign(cv.arg[i], get(cv.arg[i])[incl]) }
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 01. Calculate the solar constant (I_0), MJ m-2 hr-1
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (aprchI_0 == "Iqbal1983") { # ref: Iqbal (1983)
# 'Average' solar constant, I_SC (in MJ m-2 hr-1)
I_SC <- 4.871 # ref: 'text' in Iqbal (1983)
# Correction factor for the eccentricity of Earth’s orbit (unitless)
# ref: Eq 1.2.3 in Iqbal (1983); Eq 1.4.1 in Duffie and Beckman (1980)
E_0 <- 1. + (0.033 * cos(360. * n / 365. * pi /180.))
# Solar constant, I_0 (in MJ m-2 hr-1)
I_0 <- I_SC * E_0 # ref: Eq 4.2.1 in Iqbal (1983)
} else { # ref: Yin (1999); Brock (1981)
# 'Average' solar constant, I_SC (in MJ m-2 hr-1)
I_SC <- 4.9212 # ref: 'text' in Yin (1997b); Yin (1999)
# Radius vector of the Earth, R_1 (unitless)
# ref: Eq 2 in Brock (1981)
R_1 <- 1. / (1. + (0.033 * cos(360. * (n - 1.) / 365. * pi / 180.)))**0.5
# Solar constant, I_0 (in MJ m-2 hr-1)
I_0 <- I_SC /R_1 ** 2. # ref: 'text' in Brock (1981)
}
l.SolIrr$I_0_MJ.m2.hr1[incl] <- I_0
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 02. Calculate the solar declination (DCL), deg
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (aprchDCL == "CBM") { # ref: Eqs 1 and 2 in Forsythe et al. (1995)
DCL <- rad2deg(asin(0.39795 * cos(0.2163108 + 2. * atan(0.9671396 * tan(0.00860 * (n - 186.))))))
} else { # ref: Eq A7 in Yin (1997a); Eq 1 in Brock (1981); Eq 4 in Forsythe et al. (1995)
DCL <- 23.45 * sin(deg2rad(360. / 365. * (n + 283.)))
}
l.SolIrr$DCL_deg[incl] <- DCL
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 03. Calculate the solar zenith angle at solar noon (Z_0), deg
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ref: Eq 2.10 in Yin (1997b)
Z_0 <- ifelse(abs(lat - DCL) < 90., lat - DCL, 90.)
l.SolIrr$Z_0_deg[incl] <- Z_0
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 04. Calculate the solar zenith angle at sunset (or at solar midnight when the sun does not rise) (Z_n), deg
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ref: Eq 2.9 in Yin (1997b)
Z_n <- ifelse(abs(lat + DCL) < 90., 90.,
rad2deg(acos(sin(deg2rad(DCL)) * sin(deg2rad(lat)) - cos(deg2rad(DCL)) * cos(deg2rad(lat)))))
l.SolIrr$Z_n_deg[incl] <- Z_n
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 05. Calculate the daytime mean optical air mass (and the value of an associated experiential function)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Coefficents for calculating optical air mass
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
c_LE80 <- matrix(c(0.008307, 1.021, -0.01259), nrow = 1, dimnames = list("c_LE80", c("c1", "c2", "c3")))
c_GT80 <- matrix(c(0.03716, 1.538, -1.689), nrow = 1, dimnames = list("c_GT80", c("c1", "c2", "c3")))
coef_cho <- rbind(c_LE80, c_GT80)
# Optical air mass at solar noon (corresponding to cosine zenith angle at its maximum), m_o (unitless)
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
coef_val <- matrix(coef_cho[ifelse(Z_0 <= 80., "c_LE80", "c_GT80"), ], nrow = length(incl), ncol = ncol(coef_cho),
dimnames = list(NULL, colnames(coef_cho)))
# ref: Eq 2.1 in Yin (1997b)
m_o <- as.numeric(coef_val[, "c2"] / (coef_val[, "c1"] + cos(deg2rad(Z_0))) + coef_val[, "c3"])
# Optical air mass corresponding to cosine zenith angle at its medium, m_c (unitless)
# ref: definitions for Eq 3.3 in Yin (1997b)
Z_c <- rad2deg(acos((cos(deg2rad(Z_0)) + cos(deg2rad(Z_n))) / 2.))
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
coef_val <- matrix(coef_cho[ifelse(Z_c <= 80., "c_LE80", "c_GT80"), ], nrow = length(incl), ncol = ncol(coef_cho),
dimnames = list(NULL, colnames(coef_cho)))
# ref: Eq 2.1 in Yin (1997b)
m_c <- as.numeric(coef_val[, "c2"] / (coef_val[, "c1"] + cos(deg2rad(Z_c))) + coef_val[, "c3"])
# Optical air mass corresponding to cosine zenith angle at its bottom-quarter range point, m_l (unitless)
# ref: definitions for Eq 3.3 in Yin (1997b)
Z_l <- rad2deg(acos((cos(deg2rad(Z_0)) + 3. * cos(deg2rad(Z_n))) / 4.))
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
coef_val <- matrix(coef_cho[ifelse(Z_l <= 80., "c_LE80", "c_GT80"), ], nrow = length(incl), ncol = ncol(coef_cho),
dimnames = list(NULL, colnames(coef_cho)))
# ref: Eq 2.1 in Yin (1997b)
m_l <- as.numeric(coef_val[, "c2"] / (coef_val[, "c1"] + cos(deg2rad(Z_l))) + coef_val[, "c3"])
# Sunrise/sunset hour (solar hour at sunrise and sunset), t_1 (hr)
# ref: Eq 2.5 in Yin (1997b)
# t_1 <- ifelse(abs(lat - DCL) >= 90., 0.,
# ifelse(abs(lat + DCL) >= 90., 12.,
# (12. / pi) * acos(-1. * tan(deg2rad(lat)) * tan(deg2rad(DCL)))))
t_1 <- vector(length = length(incl))
t_1[abs(lat - DCL) >= 90. & abs(lat + DCL) < 90.] <- 0.
t_1[abs(lat - DCL) < 90. & abs(lat + DCL) >= 90.] <- 12.
sls <- abs(lat - DCL) < 90. & abs(lat + DCL) < 90.
if (!identical(sls, FALSE)) { t_1[sls] <- (12. / pi) * acos(-1. * tan(deg2rad(lat[sls])) * tan(deg2rad(DCL[sls]))) }
l.SolIrr$t_1_hr[incl] <- t_1
# Solar hour corresponding to the 80° zenith angle, t_80 (hr)
# ref: Eq 2.8 in Yin (1997b)
ratio <- (cos(deg2rad(80.)) - (sin(deg2rad(lat)) * sin(deg2rad(DCL)))) / cos(deg2rad(lat)) * cos(deg2rad(DCL))
t_80 <- vector(length = length(incl))
t_80[ratio >= 1.] <- 0.
t_80[ratio <= -1.] <- 12.
t_80[ratio < 1. & ratio > -1.] <- (1. / 15.) * rad2deg(acos(ratio[ratio < 1. & ratio > -1.]))
l.SolIrr$t_80_hr[incl] <- t_80
# Daytime mean optical air mass, m_a (unitless)
# ref: 2.7 Eq in Yin (1997b)
# print(c("lat", "DCL", "Z_0", "Z_n", "t_1", "t_80"))
# print(c(lat, DCL, Z_0, Z_n, t_1, t_80))
m_a <- calcDtMOptAM(lat, DCL, Z_0, Z_n, t_1, t_80, argCkd = TRUE)
# Checking whether use of elevation-correction for the optical air mass values
if (!is.null(elv)) {
# Station atmospheric pressure, sap (unitless)
# ref: Eq 2.11 in Yin (1997b); Linacre (1992)
sap <- exp(-elv / 8000.)
m_o <- m_o * sap
m_c <- m_c * sap
m_l <- m_l * sap
m_a <- m_a * sap
}
l.SolIrr$m_a[incl] <- m_a
# Transmission coefficient, tau (unitless)
# ref: 'text' Yin (1997b), Linacre (1992)
tau <- 0.75
l.SolIrr$tau[incl] <- tau
# An exponent depending on optical air mass, f_m (unitless)
# ref: Eq 3.3 in Yin (1997b)
f_m <- 0.2323 + 0.01452 * (m_a + m_l) * exp(1.403 * tau) - 0.1528 * m_o + m_c + 0.487 * (m_c - m_l)
l.SolIrr$f_m[incl] <- f_m
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 06. Calculate the diurnal average of the cosine of the solar zenith angle
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Half-length of the illuminated day, Hr (in radians)
# ref: Eq A2 in Yin (1997a)
Hr <- vector(length = length(incl))
Hr[abs(lat - DCL) >= 90. & abs(lat + DCL) < 90.] <- 0.
Hr[abs(lat - DCL) < 90. & abs(lat + DCL) >= 90.] <- pi
sls <- abs(lat - DCL) < 90. & abs(lat + DCL) < 90.
if (!identical(sls, FALSE)) { Hr[sls] <- acos(-1. * tan(deg2rad(lat[sls])) * tan(deg2rad(DCL[sls]))) }
# Diurnal average of the cosine of the solar zenith angle, cosZ_a (unitless)
# ref: Eq A5 in Yin (1997a)
cosZ_a <- vector(length = length(incl))
cosZ_a[abs(lat - DCL) >= 90.] <- 0.
sls <- abs(lat - DCL) < 90.
if (!identical(sls, FALSE)) { cosZ_a[sls] <- sin(deg2rad(lat[sls])) * sin(deg2rad(DCL[sls])) +
(1. / Hr[sls]) * cos(deg2rad(lat[sls])) * cos(deg2rad(DCL[sls])) * sin(Hr[sls]) }
l.SolIrr$cosZ_a[incl] <- cosZ_a
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 07. Calculate the mean hourly solar irradiance under cloudless-sky conditions (R_E), MJ m-2 hr-1
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ref: 'text' in Yin (1999); Eq 3.2 in Yin (1997b)
R_E <- I_0 * cosZ_a * (tau ** f_m)
l.SolIrr$R_E_MJ.m2.hr1[incl] <- R_E
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 08. Calculate the daylength (DL), hr
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DL <- vector(length = length(incl))
DL[abs(lat - DCL) >= 90. & abs(lat + DCL) < 90.] <- 0.
DL[abs(lat - DCL) < 90. & abs(lat + DCL) >= 90.] <- 24.
sls <- abs(lat - DCL) < 90. & abs(lat + DCL) < 90.
if (aprchDCL == "CBM") { # ref: Eqs 3 in Forsythe et al. (1995)
DL[sls] <- 24. - (24. / pi) * acos((sin(deg2rad(0.8333)) + sin(deg2rad(lat[sls])) * sin(deg2rad(DCL[sls]))) /
(cos(deg2rad(lat[sls])) * cos(deg2rad(DCL[sls]))))
} else { # ref: 'text' in Yin (1997a); Eq 4 in Brock (1981); Eqs 5 and 6 in Forsythe et al. (1995)
DL[sls] <- 24. / 180. * rad2deg(acos(-1. * tan(deg2rad(lat[sls])) * tan(deg2rad(DCL[sls]))))
}
l.SolIrr$DL_hr[incl] <- DL
# ~~~~ RETURN VALUES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
return(l.SolIrr)
}
# ********************************************************************************************************************
# Name: calcDtMOptAM
# Inputs: - double, latitude, deg (lat)
# - double, solar declination, deg (DCL)
# - double, solar zenith angle at solar noon, deg (Z_0)
# - double, solar zenith angle at sunset (or at solar midnight when the sun does not rise), deg (Z_n)
# - double, sunrise/sunset hour (solar hour at sunrise and sunset), hr (t_1)
# - double, solar hour corresponding to the 80° zenith angle, hr (t_80)
# - logical, indicates whether or not the checking and correction of arguments should be omitted (argCkd)
# Returns: - double, daytime mean optical air mass, unitless (m_a)
# Features: This function returns the daytime mean optical air mass by using the Equation 2.7 in Yin (1997b).
# Ref: - Yin X (1997b) Optical Air Mass: Daily Integration and its Applications. Meteorol Atmos Phys
# 63(3-4):227-233. DOI: 10.1007/BF01027387
# ********************************************************************************************************************
calcDtMOptAM <- function(lat, DCL, Z_0, Z_n, t_1, t_80, argCkd = FALSE) {
# ~~~~ FUNCTION WARNINGS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
cv.arg <- c("lat", "DCL", "Z_0", "Z_n", "t_1", "t_80")
for (i in 1 : length(cv.arg)) {
if (is.null(get(cv.arg[i]))) { stop("Invalid argument: '", cv.arg[i], "' is missing, with no default.") }
}
if (!argCkd) {
maVarFeatures <- data.frame(lgth = rep(1, 6), lwst_val = c(rep(-90, 2), rep(0, 4)),
hgst_val = c(rep(90, 4), rep(24, 2)),
rng_wrd = c(rep("values from the interval [-90, 90]", 2),
rep("values from the interval [0, 90]", 2),
rep("values from the interval [0, 24]", 2)))
rownames(maVarFeatures) <- cv.arg
# Checking the class and content of objects
for (i in 1 : length(cv.arg)) {
arg <- get(cv.arg[i])
if (!is.list(arg) & length(arg) >= 1L &
length(which(!is.na(match(class(arg), c("numeric", "integer"))))) >= 1L) {
if (!(is.numeric(arg) &
all(arg[!is.na(arg)] >= maVarFeatures[cv.arg[i], "lwst_val"] &
arg[!is.na(arg)] <= maVarFeatures[cv.arg[i], "hgst_val"]))) {
stop("Invalid argument: '", cv.arg[i], "' has to be a vector with ", maVarFeatures[cv.arg[i], "rng_wrd"],
".")
}
} else {
stop("Invalid argument: '", cv.arg[i], "' has to be a vector with ", maVarFeatures[cv.arg[i], "rng_wrd"],
".")
}
}
ls.arg <- stats::setNames(list(lat, DCL, Z_0, Z_n, t_1, t_80), cv.arg)
if (!all(sapply(ls.arg, length) == length(ls.arg[[1]]))) {
stop("Invalid argument: The vectors ",
paste0(paste(cv.arg[1 : (length(cv.arg) - 1)], sep = "", collapse = ", "), " and ",
cv.arg[length(cv.arg)]),
" have not the same length.")
}
}
# ~~~~ FUNCTION VARIABLES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
# Coefficents for calculating optical air mass
# ref: Eqs 2.2 and 2.3 in Yin (1997b)
c_LE80 <- matrix(c(0.008307, 1.021, -0.01259), nrow = 1, dimnames = list("c_LE80", c("c1", "c2", "c3")))
c_GT80 <- matrix(c(0.03716, 1.538, -1.689), nrow = 1, dimnames = list("c_GT80", c("c1", "c2", "c3")))
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 01. Select the equation(s) applied
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
eq <- stats::setNames(lapply(vector("list", 4), function(x) vector()), c("a", "b", "c", "d"))
eq$a <- which(t_1 == 0.)
eq$b <- which(t_1 != 0. & Z_n <= 80.)
eq$c <- which(t_1 != 0. & Z_0 >= 80.)
eq$d <- which(t_1 != 0. & Z_n > 80. & Z_0 < 80.)
applEq <- which(lapply(eq, function(x) { !identical(x, integer(0)) } ) == T)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 02. Set the result vector
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
m_a <- rep(NA, length(lat))
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 03. Run the algorithm for each equation type
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for (i_eq in applEq) {
slctd <- eq[[c("a", "b", "c", "d")[i_eq]]]
sn <- 0
while(sn <= 2) {
if (c("a", "b", "c", "d")[i_eq] == "a") {
m_a[slctd] <- 39.7
break
}
if (c("a", "b", "c", "d")[i_eq] == "b") {
t_1g <- t_1[slctd]
cg <- c_LE80
}
if (c("a", "b", "c", "d")[i_eq] == "c") {
t_1g <- t_1[slctd]
cg <- c_GT80
}
if (c("a", "b", "c", "d")[i_eq] == "d") {
sn <- sn + 1
if (sn == 1) t_1g <- t_80[slctd]
if (sn == 2) t_1g <- t_1[slctd]
if (sn == 3) t_1g <- t_80[slctd]
if (sn == 1) cg <- c_LE80
if (sn == 2) cg <- c_GT80
if (sn == 3) cg <- c_GT80
}
a_ar <- cg[1] + sin(deg2rad(lat[slctd])) * sin(deg2rad(DCL[slctd]))
if (sn < 2) {
b_ar <- cos(deg2rad(lat[slctd])) * cos(deg2rad(DCL[slctd]))
}
mbr2 <- mbr3 <- mbr4 <- rep(NA, length(slctd))
sls <- a_ar > b_ar
if (!identical(sls, integer(0))) {
mbr2[sls] <- cg[2] / sqrt(a_ar[sls] ** (2.) - b_ar[sls] ** (2.))
numer <- b_ar[sls] + a_ar[sls] * cos(deg2rad(15. * t_1g[sls]))
denom <- a_ar[sls] + b_ar[sls] * cos(deg2rad(15. * t_1g[sls]))
mbr3[sls] <- acos(numer / denom)
}
sls <- a_ar < b_ar
if (!identical(sls, integer(0))) {
mbr2[sls] <- cg[2] / sqrt(b_ar[sls] ** (2.) - a_ar[sls] ** (2.))
sqr1 <- sqrt((b_ar[sls] + a_ar[sls]) * (1. + cos(deg2rad(15. * t_1g[sls]))))
sqr2 <- sqrt((b_ar[sls] - a_ar[sls]) * (1. - cos(deg2rad(15. * t_1g[sls]))))
mbr3[sls] <- log((sqr1 + sqr2) / (sqr1 - sqr2))
}
sls <- a_ar == b_ar
if (!identical(sls, integer(0))) {
mbr2[sls] <- cg[2] / a_ar[sls]
mbr3[sls] <- tan(deg2rad(15. * t_1g[sls] / 2.))
}
mbr4 <- cg[3] * t_1g
F_eq6 <- 1. / deg2rad(15.) * mbr2 * mbr3 + mbr4
if (sn == 0) { F0_eq6 <- F_eq6 }
if (sn == 1) { F1_eq6 <- F_eq6 }
if (sn == 2) { F2_eq6 <- F_eq6 }
if (sn == 3) { F3_eq6 <- F_eq6 }
if (sn == 1 | sn == 2) { next }
if (sn == 0) {
m_a[slctd] <- 1. / t_1g * F0_eq6
break
}
if (sn == 3) {
m_a[slctd] <- 1. / t_1[slctd] * (F1_eq6 + F2_eq6 - F3_eq6)
break
}
}
}
# ~~~~ RETURN VALUES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
return(m_a)
}
deg2rad <- function(deg) { deg * (pi / 180.) }
rad2deg <- function(rad) { rad * (180. / pi) }
Try the macroBiome package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.