Add chilling hours and clearsky_radiation functions

NAMESPACE

@@ -2,6 +2,8 @@
 
 export("%>%")
 export(cdd)
+export(chillinghours)
+export(clearsky_radiation)
 export(convert_units)
 export(et0pm)
 export(hdd)

R/chillinghours.R

@@ -1,3 +1,20 @@
+#' Compute chilling hours
+#'
+#' Function to compute chilling hours based on the daily temperature data.
+#'
+#' @param t_h numeric. A vector of daily temperature data.
+#' @param threshold numeric. A vector of two values representing the lower and upper threshold temperatures.
+#' @param tunit character. The unit of the temperature data. Default is "celsius".
+#' @param thresholdunit character. The unit of the threshold temperature. Default is "celsius".
+#' @param output character. The type of output. Default is "chillhours".
+#' @param na.rm logical. Should missing values be removed? Default is TRUE.
+#'
+#' @return A numeric value representing the chilling hours or chilling units.
+#'
+#' @examples
+#' chillinghours(c(5, 6, 7, 8, 9, 10))
+#'
+#' @export
 chillinghours = function(t_h, threshold = c(0,7.2),
                          tunit = "celsius",
                          thresholdunit = "celsius",

R/clearsky_radiation.R

@@ -0,0 +1,28 @@
+#' Calculate the clear sky radiation.
+#'
+#' This function computes the clear sky radiation using the formula proposed by Allen et al. (1998).
+#'
+#' @param j numeric. Day of the year.
+#' @param lat numeric. Latitude in radians.
+#' @param elev numeric. Elevation in meters.
+#'
+#' @return A numeric value of the clear sky radiation in W/m2/day.
+#'
+#' @examples
+#' # date='2024-03-09'
+#' # lat=-30.631261
+#' # elev=80
+#'
+#' lat = convert_units(-30.631261, "angle", "deg", "rad")
+#' clearsky_radiation(j=69, lat=lat, elev=80)
+#'
+#' @export
+clearsky_radiation = function(j, lat, elev){
+  dr = 1 + 0.033 * cos(2 * pi * j / 365)
+  delta_s = 0.409 * sin((2 * pi * j / 365) - 1.39)
+  ws = acos(-tan(lat) * tan(delta_s))
+  ra = (24 * 60 / pi) * 0.0820 * dr * (ws * sin(lat) * sin(delta_s) +
+                                         cos(lat) * cos(delta_s) * sin(ws))
+  rso = (0.75 + 0.00002 * elev) * ra
+  return(rso)
+}

R/et0pm.R

@@ -174,13 +174,7 @@ et0pm = function(date,
     }
   }
   j = as.numeric(format(date, format = "%j"))
-  dr = 1 + 0.033 * cos(2 * pi * j / 365)
-  delta_s = 0.409 * sin((2 * pi * j / 365) - 1.39)
-  # lat_r = lat * pi / 180
-  ws = acos(-tan(lat) * tan(delta_s))
-  ra = (24 * 60 / pi) * 0.0820 * dr * (ws * sin(lat) * sin(delta_s) +
-                                         cos(lat) * cos(delta_s) * sin(ws))
-  rso = (0.75 + 0.00002 * elev) * ra
+  rso = clearsky_radiation(j, lat, elev)
   rns = (1 - 0.23) * r_s
   rnl = 0.000000004903 * (((t_max + 273.16) ^ 4 + (t_min + 273.16) ^ 4) /
                             2) * (0.34 - 0.14 * sqrt(ea)) * ((1.35 * r_s / rso) - 0.35)