R/solar.R

In FabianMitze/FieldClimRevised: Calculates Fluxes for Weather Stations

#' Eccentricity
#'
#' Calculates the eccentricity.
#'
#' @param ... Additional parameters passed to later functions.
#' @return Eccentricity at the date.
#' @export
#'
sol_eccentricity <- function (...) {
  UseMethod("sol_eccentricity")
}

#' @rdname sol_eccentricity
#' @method sol_eccentricity POSIXt
#' @param datetime POSIXt object (POSIXct, POSIXlt).
#' See [base::as.POSIXlt] and [base::strptime] for conversion.
#' @export
sol_eccentricity.POSIXt <- function(datetime, ...) {

  if(!inherits(datetime, "POSIXt")){
    stop("datetime has to be of class POSIXt.")
  }

  # day of year
  doy <- as.numeric(strftime(datetime, format = "%j"))

  x <- 2.0* pi * (doy - 1) / 365.0
  exz <- 1.00011 + 0.034221 * cos(x) + 0.00128 * sin(x) + 0.000719 * cos(2 * x) + 0.000719 * sin(2. * x)
  return(exz)
}

#' @rdname sol_eccentricity
#' @method sol_eccentricity weather_station
#' @param weather_station Object of class weather_station.
#' @export
sol_eccentricity.weather_station <- function(weather_station, ...) {
  check_availability(weather_station, "datetime")

  datetime <- weather_station$measurements$datetime

  return(sol_eccentricity(datetime))
}



#' Solar azimuth and elevation angles
#'
#' Calculates solar azimuth and solar elevation angle.
#'
#' @param ... Additional parameters passed to later functions.
#' @return data.frame with two columns: sol_azimuth and sol_elevation.
#' @export
#'
sol_angles <- function(...) {
  UseMethod("sol_angles")
}

#' @rdname sol_angles
#' @method sol_angles POSIXt
#' @param datetime POSIXt object (POSIXct, POSIXlt).
#' See [base::as.POSIXlt] and [base::strptime] for conversion.
#' @param lat Latitude in decimal degrees.
#' @param lon Longitude in decimal degrees.
#' @export
sol_angles.POSIXt <- function(datetime, lat, lon, ...){

  if(!inherits(datetime, "POSIXt")){
    stop("datetime has to be of class POSIXt.")
  }
  # day of year
  doy <- as.numeric(strftime(datetime, format = "%j"))

  # decimal hour
  lt <- as.POSIXlt(datetime)
  ut <- lt$hour + lt$min / 60 + lt$sec / 3600

  f <- pi / 180   # angle to radian

  # conversion latitude to degrees
  gbr <- lat*f
  glr <- lon*f

  # medium sun time
  t <- ut + lon / 15. # in hours

  # angle in radian
  m <- 356.6 + 0.9856 * doy # in degrees
  m <- m*f                  # in radian

  # time formula
  zt <- 0.1644 * sin(2.* glr) - 0.1277 * sin(m) # in hours

  # hour-angle of the sun
  h <- (15.* f) * (t + zt - 12.)  # in radian

  # geocentric apparent ecliptic longitude of the sun (in radian)
  del <- 279.3 * f + 0.9856 * f * doy + 1.92 * f * sin(356.6 * f + 0.9856 * f * doy)

  # Sine declination of the sun (in radian)
  sde <- sin(23.44 * f) * sin(del)

  # sun height
  shh <- sin(gbr) * sde + cos(gbr) * cos(asin(sde)) * cos(h)
  sh  <- asin(shh) / f # (in degrees)

  # Sun azimuth
  saz <- (sde * cos(gbr) - cos(asin(sde)) * sin(gbr) * cos(h)) / cos((sh * f))

  saz_2 <- rep(NA, length(saz))
  for(i in 1:length(saz_2)){
    if(t[i] <= 12){saz_2[i] <- acos(saz[i])}
    if(t[i]  > 12){saz_2[i] <- 360 * f - acos(saz[i])}
  }
  saz_deg <- saz_2/f
  results <- data.frame(sol_azimuth = saz_deg,
                        sol_elevation = sh)
  return(results)
}

#' @rdname sol_angles
#' @method sol_angles weather_station
#' @param weather_station Object of class weather_station.
#' @export
sol_angles.weather_station <- function(weather_station, ...) {
  check_availability(weather_station, "datetime", "latitude", "longitude")

  datetime <- weather_station$measurements$datetime
  lat <- weather_station$location_properties$latitude
  lon <- weather_station$location_properties$longitude

  return(sol_angles(datetime, lat, lon))
}


#' Solar elevation angle
#'
#' Calculates solar elevation angle for the given date and time.
#'
#' @param ... Additional parameters passed to later functions.
#' @return Solar elevation angle in degrees.
#' @export
#'
sol_elevation <- function (...) {
  UseMethod("sol_elevation")
}

#' @rdname sol_elevation
#' @method sol_elevation POSIXt
#' @param datetime POSIXt object (POSIXct, POSIXlt).
#' See [base::as.POSIXlt] and [base::strptime] for conversion.
#' @param lat Latitude in decimal degrees.
#' @param lon Longitude in decimal degrees.
#' @export
sol_elevation.POSIXt <- function(datetime, lat, lon, ...) {
  angles <- sol_angles(datetime, lat, lon)
  return(angles$sol_elevation)
}

#' @rdname sol_elevation
#' @method sol_elevation weather_station
#' @param weather_station Object of class weather_station.
#' @export
sol_elevation.weather_station <- function(weather_station, ...) {
  check_availability(weather_station, "datetime", "latitude", "longitude")

  datetime <- weather_station$measurements$datetime
  lat <- weather_station$location_properties$latitude
  lon <- weather_station$location_properties$longitude

  return(sol_elevation(datetime, lat, lon))
}



#' Solar azimuth angle
#'
#' Calculates solar azimuth angle for the given date and time.
#'
#' @param ... Additional parameters passed to later functions.
#' @return Solar azimuth angle in degrees.
#' @export
#'
sol_azimuth <- function (...) {
  UseMethod("sol_azimuth")
}

#' @rdname sol_azimuth
#' @method sol_azimuth POSIXt
#' @param datetime POSIXt object (POSIXct, POSIXlt).
#' See [base::as.POSIXlt] and [base::strptime] for conversion.
#' @param lat Latitude in decimal degrees.
#' @param lon Longitude in decimal degrees.
#' @export
sol_azimuth.POSIXt <- function(datetime, lat, lon, ...) {
  angles <- sol_angles(datetime, lat, lon)
  return(angles$sol_azimuth)
}

#' @rdname sol_azimuth
#' @method sol_azimuth weather_station
#' @param weather_station Object of class weather_station.
#' @export
sol_azimuth.weather_station <- function(weather_station, ...) {
  check_availability(weather_station, "datetime", "latitude", "longitude")

  datetime <- weather_station$measurements$datetime
  lat <- weather_station$location_properties$latitude
  lon <- weather_station$location_properties$longitude

  return(sol_azimuth(datetime, lat, lon))
}