diff --git a/R/PEdaily_Oudin.R b/R/PEdaily_Oudin.R index c677966407f1255b33aa3d128c3d8d1d0f415034..b71ac99926c6d2725c2ef712575d2cee14a69cb9 100644 --- a/R/PEdaily_Oudin.R +++ b/R/PEdaily_Oudin.R @@ -1,4 +1,4 @@ -PEdaily_Oudin <- function(JD, Temp, LatRad, Lat, LatUnit = c("rad", "deg"), run_fortran = F) { +PEdaily_Oudin <- function(JD, Temp, LatRad, Lat, LatUnit = c("rad", "deg")) { ## ---------- deprecated function @@ -28,9 +28,9 @@ PEdaily_Oudin <- function(JD, Temp, LatRad, Lat, LatUnit = c("rad", "deg"), run_ if (!any(LatUnit %in% c("rad", "deg"))) { stop("'LatUnit' must be one of \"rad\" or \"deg\"") } - # if (!inherits(Lat, "numeric") | length(Lat) != 1) { - # stop("'Lat' must be a 'numeric' of length one") - # } + if (!inherits(Lat, "numeric") | length(Lat) != 1) { + stop("'Lat' must be a 'numeric' of length one") + } if (LatUnit[1L] == "rad" & ((Lat >= pi/2) | (Lat <= -pi/2))) { stop("'Lat' must be comprised between -pi/2 and +pi/2 degrees") } @@ -49,93 +49,67 @@ PEdaily_Oudin <- function(JD, Temp, LatRad, Lat, LatUnit = c("rad", "deg"), run_ ## ---------- Oudin's formula - - if (run_fortran) - { - if (LatUnit[1L] == "rad") { - Lat = Lat * 180 / pi - } + + PE_Oudin_D <- rep(NA, length(Temp)) + COSFI <- cos(FI) + AFI <- abs(FI / 42) + + for (k in seq_along(Temp)) { - LInputs = as.integer(length(Temp)) + TETA <- 0.4093 * sin(JD[k] / 58.1 - 1.405) + COSTETA <- cos(TETA) + COSGZ <- max(0.001, cos(FI - TETA)) + GZ <- acos(COSGZ) + COSOM <- 1 - COSGZ / COSFI / COSTETA - if (length(Lat) == 1) { - Lat = rep(Lat, LInputs) + if (COSOM < -1) { + COSOM <- -1 + } + if (COSOM > 1) { + COSOM <- 1 } - RESULTS <- .Fortran("frun_etp_oudin", PACKAGE="airGR", - ##inputs - LInputs = LInputs, - InputsLAT = as.double(Lat), - InputsTT = as.double(Temp), - InputsJJ = as.double(JD), - ##outputs - PE_Oudin_D = rep(as.double(-999.999), LInputs) - ) - PE_Oudin_D = RESULTS$PE_Oudin_D - } else { + COSOM2 <- COSOM * COSOM + + if (COSOM2 >= 1) { + SINOM <- 0 + } else { + SINOM <- sqrt(1 - COSOM2) + } - PE_Oudin_D <- rep(NA, length(Temp)) - COSFI <- cos(FI) - AFI <- abs(FI / 42) + OM <- acos(COSOM) + COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1) - for (k in seq_along(Temp)) { - - TETA <- 0.4093 * sin(JD[k] / 58.1 - 1.405) - COSTETA <- cos(TETA) - COSGZ <- max(0.001, cos(FI - TETA)) - GZ <- acos(COSGZ) - COSOM <- 1 - COSGZ / COSFI / COSTETA - - if (COSOM < -1) { - COSOM <- -1 - } - if (COSOM > 1) { - COSOM <- 1 - } - - COSOM2 <- COSOM * COSOM - - if (COSOM2 >= 1) { - SINOM <- 0 - } else { - SINOM <- sqrt(1 - COSOM2) - } - - OM <- acos(COSOM) - COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1) - - if (COSPZ < 0.001) { - COSPZ <- 0.001 - } - - ETA <- 1 + cos(JD[k] / 58.1) / 30 - GE <- 446 * OM * COSPZ * ETA - - if (is.na(Temp[k])) { - PE_Oudin_D[k] <- NA - } else { - if (Temp[k] >= -5.0) { - PE_Oudin_D[k] <- GE * (Temp[k] + 5) / 100 / 28.5 - } else { - PE_Oudin_D[k] <- 0 - } - } - + if (COSPZ < 0.001) { + COSPZ <- 0.001 } - if (any(is.na(Temp))) { - if (any(is.na(PE_Oudin_D))) { - warning("'Temp' time series, and therefore the returned 'PE' time series, contain missing value(s)") - } else { - warning("'Temp' time series contains missing value(s)") - } + ETA <- 1 + cos(JD[k] / 58.1) / 30 + GE <- 446 * OM * COSPZ * ETA + + if (is.na(Temp[k])) { + PE_Oudin_D[k] <- NA + } else { + if (Temp[k] >= -5.0) { + PE_Oudin_D[k] <- GE * (Temp[k] + 5) / 100 / 28.5 + } else { + PE_Oudin_D[k] <- 0 } - if (!any(is.na(Temp)) & any(is.na(PE_Oudin_D))) { - warning("returned 'PE' time series contains missing value(s)") } } + if (any(is.na(Temp))) { + if (any(is.na(PE_Oudin_D))) { + warning("'Temp' time series, and therefore the returned 'PE' time series, contain missing value(s)") + } else { + warning("'Temp' time series contains missing value(s)") + } + } + if (!any(is.na(Temp)) & any(is.na(PE_Oudin_D))) { + warning("returned 'PE' time series contains missing value(s)") + } + return(PE_Oudin_D) }