
###############################################################.

# Olmo slope correction + DEM horizon mask

# R implementation of the Olmo et al. (1999) correction to get global short-wave on an inclined plane from measured global horizontal irradiance (G_H), plus a small optional ground-reflection term. I also include (1) how we compute the clearness index 
# kt (Olmo needs it) using extraterrestrial horizontal irradiance, (2) solar geometry (sun zenith & azimuth), and (3) brief notes linking the implementation to the papers cited (Olmo; Leuning et al.; McGloin et al.).
# 
# Olmo equation (Eq. 4):  G(ω) = GH * exp(−kt(ω2 −ωH^2))
# 
# where 
# ω is the angular distance between the sun direction and the plane normal (in radians), and 
# ωH is the angular distance between the sun and the horizontal plane (i.e. the solar zenith angle). The code computes 
# ω from plane tilt/azimuth and sun zenith/azimuth (Olmo et al., 1999). See citations after the code. 

# -------- DEM horizon mask -----------.
# integrate a DEM horizon extraction step directly into the function.
# input:
# . a raster DEM file (or a terra/raster object),
# . your station coordinates (lat/lon),
# and it will calculate the horizon profile automatically.
# it:
# . Reproject the DEM to a suitable projected CRS (so distances/azimuths are correct).
# . For each azimuth from 0–359°, cast a radial line from the station point.
# . Sample DEM along the radial line and compute the maximum elevation angle from the station’s perspective (arctan(height diff / horizontal distance)).
# Then that profile is fed directly into the Olmo correction logic.
# the function:
# . Computes sun position internally (no external solar package)
# . Auto-calculates horizon elevation profile from your DEM
# . Masks irradiance when the sun is below that local horizon

# The correction is performed accounting for the sensor’s mounting height above the local ground.

# References
# . Olmo, F.J., Vida, J., Foyo, I., Castro-Diez, Y., Alados-Arboledas, L., 1999. Prediction of global irradiance on inclined surfaces from horizontal global irradiance. Energy 24, 689–704. doi:10.1016/S0360-5442(99)00025-0
# . Leuning, R., van Gorsel, E., Massman, W.J., Isaac, P.R., 2012. Reflections on the surface energy imbalance problem. Agric. For. Meteorol. 156, 65–74. doi:10.1016/j.agrformet.2011.12.002
# . McGloin, R., Šigut, L., Havránková, K., Dušek, J., Pavelka, M., Sedlák, P., 2018. Energy balance closure at a variety of ecosystems in Central Europe with contrasting topographies. Agric. For. Meteorol. 248, 418–431. doi:10.1016/j.agrformet.2017.10.003

###############################################################.




if (!requireNamespace("terra", quietly = TRUE)) install.packages("terra")
library(terra)

#' Olmo slope correction with automatic DEM horizon shadowing
#'
#' @param datetime POSIXct vector in UTC
#' @param Gh numeric vector, global horizontal irradiance (W/m2)
#' @param lat station latitude (decimal degrees, N+)
#' @param lon station longitude (decimal degrees, E+)
#' @param tilt plane tilt (degrees) >> SLOPE
#' @param azimuth_plane plane azimuth (degrees, 0=N, clockwise)  >> ASPECT
#' @param albedo ground albedo
#' @param add_reflection logical
#' @param dem_corr perform the DEM horizon mask: whether to force G values to 0 when the sensor is shadowed 
#' @param dem optional: filename or SpatRaster of DEM
#' @param max_dist maximum search distance for horizon (m)
#' @param az_step azimuth step for horizon calculation (degrees)
#' @param sensor_height sensor height above ground (m)
#' @param night_control whether to force nigh time G values to 0
#' 
#'
#' @return data.frame with tilt-corrected irradiance and sun visibility flag

rad_topo_correction <- function(datetime, # POSIXct vector in UTC
                                Gh, #  numeric vector, global horizontal irradiance (W/m2)
                                lat, #  station latitude (decimal degrees, N+)
                                lon, #  station longitude (decimal degrees, E+)
                                tilt, #  plane tilt (degrees) >> SLOPE
                                azimuth_plane, #  plane azimuth (degrees, 0=N, clockwise)  >> ASPECT
                                albedo, #  ground albedo
                                add_reflection = TRUE, #  logical
                                dem_corr = FALSE, # the DEM horizon mask: whether to force G values to 0 when the sensor is shadowed 
                                dem = NULL, # optional, filename or SpatRaster of DEM
                                max_dist = 5000,  # maximum search distance for horizon (m)
                                az_step = 1, #  azimuth step for horizon calculation (degrees)
                                sensor_height = 0, #  sensor height above ground (m)
                                night_control = TRUE #  whether to force nigh time Gh values to 0: TRUE: force SW_IN to 0 during night, FALSE leave the SW_IN untouched
) 
{

  stopifnot(length(datetime) == length(Gh))
  dt_utc <- as.POSIXct(datetime, tz = "UTC")
  doy <- as.numeric(format(dt_utc, "%j"))
  
  
  # ---- Solar position ----
  
  # Convert to fractional hour UTC
  y <- as.POSIXlt(dt_utc, tz = "UTC")$year + 1900
  m <- as.POSIXlt(dt_utc, tz = "UTC")$mon + 1
  d <- as.POSIXlt(dt_utc, tz = "UTC")$mday
  H <- as.POSIXlt(dt_utc, tz = "UTC")$hour +
    as.POSIXlt(dt_utc, tz = "UTC")$min / 60 +
    as.POSIXlt(dt_utc, tz = "UTC")$sec / 3600
  
  # Julian Day
  jd <- (1461 * (y + 4800 + (m - 14) %/% 12)) %/% 4 +
    (367 * (m - 2 - 12 * ((m - 14) %/% 12))) %/% 12 -
    (3 * ((y + 4900 + (m - 14) %/% 12) %/% 100)) %/% 4 +
    d - 32075 + H / 24
  n <- jd - 2451545.0
  
  # Mean longitude anomaly
  # distanza angolare dal pericentro che un corpo avrebbe se si movesse in un'orbita circolare, a velocità costante e con lo stesso periodo orbitale di un corpo in orbita ellittica
  L <- (280.46 + 0.9856474 * n) %% 360
  g <- (357.528 + 0.9856003 * n) %% 360
  g_rad <- g * pi/180
  lambda <- (L + 1.915 * sin(g_rad) + 0.020 * sin(2 * g_rad)) %% 360
  lambda_rad <- lambda * pi/180
  
  # Obliquity
  # angolo rispetto all'equatore celeste del cerchio massimo di riferimento delle coordinate eclittiche, oggi pari a 23°27'
  eps <- 23.439 - 0.0000004 * n
  eps_rad <- eps * pi/180
  
  # Ascensione retta
  # distanza angolare fra il meridiano fondamentale e il meridiano passante per l'oggetto scelto
  # in prtica è la longitudine, ma proiettata sulla sfera celeste anziché sulla superficie terrestre
  alpha <- atan2(cos(eps_rad) * sin(lambda_rad), cos(lambda_rad))
  
  # Sun declination (OK)
  # è l'angolo che i raggi del Sole formano con il piano equatoriale terrestre, positivo se il Sole è al di sopra del piano stesso
  # insieme all'ascensione retta, serve per determinare l'altezza di un astro sulla sfera celeste
  # in pratica è la latitudine proiettata sulla sfera celeste, anziché sulla superficie terrestre
  delta <- asin(sin(eps_rad) * sin(lambda_rad))  # (d.r in codice McGloin)
  
  # Equation of time (minutes)
  # lo scostamento, durante il corso dell'anno, tra il tempo indicato da un orologio solare e il suo scorrere regolare indicato da un orologio
  E_time <- L/15 - (alpha * 180/pi)/15
  
  # Hour angle
  # lo spostamento angolare del Sole a est o a ovest del meridiano locale a causa della rotazione della terra sul suo asse a 15° per ora, con la mattina negativa e il pomeriggio positivo
  solar_time <- (H + E_time) %% 24
  H_angle <- (solar_time - 12) * 15 * pi/180   # (sha.r in codice McGloin)
  
  # Solar zenith angle
  # The solar zenith angle is the zenith angle of the sun, i.e., the angle between the sun’s rays and the vertical direction. 
  # It is the complement to the solar altitude or solar elevation, which is the altitude angle or elevation angle between the sun’s rays and a horizontal plane.
  # normally used in combination with the solar azimuth angle to determine the position of the Sun as observed from a given location on the surface of the Earth
  lat_rad <- lat * pi/180
  cos_zen <- sin(lat_rad) * sin(delta) + cos(lat_rad) * cos(delta) * cos(H_angle)    # (E.t2 in codice McGloin)
  cos_zen <- pmin(pmax(cos_zen, -1), 1)
  sun_zenith_rad <- acos(cos_zen)   # (z in codice McGloin)
  sun_zenith_deg <- sun_zenith_rad * 180/pi  
  
  # Solar azimuth angle (0=N, clockwise)
  # The solar azimuth angle is the azimuth (horizontal angle with respect to north) of the Sun's position. This horizontal coordinate defines the Sun's relative direction along the local horizon
  sin_az <- -sin(H_angle) * cos(delta) / sin(sun_zenith_rad)
  cos_az <- (sin(delta) - sin(lat_rad) * cos_zen) / (cos(lat_rad) * sin(sun_zenith_rad))
  sun_azimuth_rad <- atan2(sin_az, cos_az)
  sun_azimuth_deg <- (sun_azimuth_rad * 180/pi + 360) %% 360
  
  # --- extraterrestrial correction factor E0 (earth-sun distance) ---
  E0 <- 1 + 0.033 * cos(360 * doy / 365) # McGloin formulation
  # Solar constant (W/m2). literature uses 1367..1361; we use 1367 (Olmo & many authors)
  I_sc <- 1367
  
  
  # --- Extraterrestrial horizontal irradiance G0h ---
  G0h <- I_sc * E0 * pmax(cos(sun_zenith_rad), 0)
  
  # --- Clearness index ---
  
  # calculated using the ratio of K↓ and estimates of the extra-terrestrial radiation
  kt <- ifelse(G0h > 1e-6,  # avoid division by zero
               Gh / G0h,
               NA)          # using NA here allow to not automatically force to 0 nightime values of SW_IN (Gh). This can be handled later during the correction
  # kt <- Gh / G0h
  kt <- replace(kt, kt > 1, 1) # (k in codice McGloin)

  # Convert kt into daily mean values as as done by McGloin
  df.tmp <- data.frame(datetime_hh=datetime, kt=kt)
  df.tmp.km <- df.tmp %>%
    mutate(days = format(datetime_hh, '%F')) %>%
    group_by(days) %>%
    mutate(km = mean(kt, na.rm=T)) %>%
    # mutate(km = quantile(kt, 0.75, na.rm=T)) %>%
    ungroup()
  kt.dm <- df.tmp.km$'km'

  rm(df.tmp, df.tmp.km)
  
  
  # ---- Tilt geometry ----
  
  # Zenith angle for inclined slope, plane normal zenith (i.e. tilt)
  beta_rad <- tilt * pi/180  # (zs in codice McGloin)
  
  # Azimuth angle (aspect) [rad]
  alpha_plane_rad <- azimuth_plane * pi/180
  
  
  ## --- Angular distance between sun and plane normal (omega) ----
  # angolo ω tra la normale del piano e la direzione del Sole
  # cosω = sinβ sinθz + cosβ cosθz cos(φs − φp)  # Eq. 2 in Olmo 
  # 𝛽= tilt del piano,
  # 𝜑s = azimut solare,
  # 𝜑p = azimut del piano

  # Codice McGloin modified ...........................
  
  # term 1
  B.t1 <- (sin(lat_rad)*cos(H_angle)*(-cos(alpha_plane_rad)*sin(beta_rad)) - sin(H_angle)*(sin(alpha_plane_rad)*sin(beta_rad)) + (cos(lat_rad)*cos(H_angle))*cos(beta_rad))*cos(delta)
  # term 2
  B.t2 <- (cos(lat_rad)*(cos(alpha_plane_rad)*sin(beta_rad)) + sin(lat_rad)*cos(beta_rad))*sin(delta)

  # Angular distance between the solar beam and the normal to the slope
  cos_omega <- B.t1 + B.t2
  
  # clip numeric noise
  cos_omega <- pmin(pmax(cos_omega, -1), 1)

  # Angular distance
  omega_rad <- acos(cos_omega)
  # ...........................
  
  
  # Solar zenith angular distance to horizontal plane is omega_H == solar zenith angle
  omega_H_rad <- sun_zenith_rad
  
  # ---- Olmo slope correction ----
  # Olmo formula: G_omega = G_H * exp(-k_t * (omega^2 - omega_H^2)) ---
  
  # G_tilt_olmo <- Gh * exp(-kt * (omega_rad^2 - omega_H_rad^2)) # using half-hourly kt
  G_tilt_olmo <- Gh * exp(-kt.dm * (omega_rad^2 - omega_H_rad^2)) # using daily kt
  
  # . Nighttime correction:
  # TRUE: force G to 0 during night
  # FALSE: leave the SW_IN untouched
  if (night_control) {
    # G_tilt_olmo[sun_zenith_deg >= 90 | is.na(Gh)] <- 0  # this set to 0 possible daily NAs in SW_IN
    G_tilt_olmo[sun_zenith_deg >= 90] <- 0
  } else {
    G_tilt_olmo[sun_zenith_deg >= 90] <- Gh[sun_zenith_deg >= 90] 
  }
  
  
  ## --- Ground reflection (isotropic) ----
  # classical isotropic ground-reflection on tilted plane: G_ref = rho * G_H * (1 - cos(beta))/2
  if (add_reflection) {
    
    # using eq. 5 in Olmo 1999: G_tilt_olmo_alb <- G_tilt_olmo * Fc
    Fc <- 1 + (albedo/100)*(sin(omega_rad/2))^2
    G_tilt_olmo_alb <- G_tilt_olmo * Fc
    
  } else {
    
    Fc <- rep(NA, length(G_tilt_olmo))
    G_tilt_olmo_alb <- G_tilt_olmo
    
  }
  
  
  # ---- DEM horizon profile ----
  
  horizon_profile <- NULL
  
  if (!is.null(dem)) {
    
    if (inherits(dem, "character")) dem <- terra::rast(dem)
    
    # station point
    pt_ll <- terra::vect(data.frame(lon = lon, lat = lat),
                         geom = c("lon", "lat"), crs = "EPSG:4326")
    
    # reproject DEM to projected CRS for distance calc (UTM)
    utm_crs <- paste0("+proj=utm +zone=", floor((lon + 180)/6) + 1,
                      ifelse(lat >= 0, " +north", " +south"), " +datum=WGS84")
    dem_utm <- terra::project(dem, utm_crs)
    pt_utm <- terra::project(pt_ll, utm_crs)
    st_xy <- terra::crds(pt_utm)
    
    azimuths <- seq(0, 359, by = az_step)
    elev_angles <- numeric(length(azimuths))
    
    for (i in seq_along(azimuths)) {
      az <- azimuths[i] * pi/180
      # create points along the ray
      dists <- seq(10, max_dist, by = 10)
      x <- st_xy[1] + dists * sin(az)
      y <- st_xy[2] + dists * cos(az)
      pts <- terra::vect(cbind(x, y), crs = utm_crs)
      elevs <- terra::extract(dem_utm, pts)[,2]
      elev_station <- terra::extract(dem_utm, pt_utm)[,2]
      # accounting for the sensor’s mounting height above the local ground
      angles <- atan2(elevs - (elev_station + sensor_height), dists) * 180/pi
      
      elev_angles[i] <- max(angles, na.rm = TRUE)
    }
    horizon_profile <- data.frame(azimuth_deg = azimuths,
                                  elev_angle_deg = elev_angles)
  }
  
  
  # ---- Horizon mask ----
  # TRUE points under sun, FALSE point shadowed
  in_sun <- rep(TRUE, length(Gh))
  if (!is.null(horizon_profile)) {
    # nearest horizon elevation for each sun azimuth
    idx <- match(round(sun_azimuth_deg) %% 360, 
                 horizon_profile$azimuth_deg)
    hor_elev <- horizon_profile$elev_angle_deg[idx]
    hor_elev[is.na(hor_elev)] <- 0
    in_sun <- sun_zenith_deg < (90 - hor_elev)
  }
  
  
  
  # . Topography shadowing correction
  G_tilt_olmo_dem <- G_tilt_olmo_alb
  if (dem_corr) {
    G_tilt_olmo_dem[!in_sun] <- 0
  }
  
  
  
  # Output ...................
  return(data.frame(
    datetime = dt_utc,
    Gh = Gh,                            # measured irradiance
    G0h = G0h,                          # Extra terrsestrial irradiance
    kt = kt,                            # clearness index
    sun_zenith_deg = sun_zenith_deg,        # solar zenith angle [°]
    sun_azimuth_deg = sun_azimuth_deg,           # solar azimuth angle [°]
    omega_deg = omega_rad * 180/pi,     # angular distance between the solar beam and the normal to the slope
    G_tilt_olmo = G_tilt_olmo,          # slope corrected irradiance
    G_reflection = Fc,                  # reflection factor
    G_tilt_alb = G_tilt_olmo_alb,       # slope and reflection corrected irradiance
    in_sun = in_sun,                    # insulation flag
    G_tilt_total = G_tilt_olmo_dem      # slope, reflection and topographic shadowing corrected irradiance
  ))
  
}



# # -------------------------.
# # Example usage:
# # -------------------------.
# 
# # Example without DEM
# res <- rad_topo_correction(
#   datetime = as.POSIXct("2023-06-21 08:00:00", tz = "UTC") + 0:5*3600,
#   Gh = c(0, 150, 400, 600, 800, 900),
#   lat = 45, lon = 10,
#   tilt = 30, azimuth_plane = 180
# )

# # Example with DEM
# # (Assume "dem.tif" is a high-res DEM covering your site)
# res_dem <- rad_topo_correction(
#   datetime = as.POSIXct("2023-06-21 08:00:00", tz = "UTC") + 0:5*3600,
#   Gh = c(0, 150, 400, 600, 800, 900),
#   lat = 45, lon = 10,
#   tilt = 30, azimuth_plane = 180,
#   dem = "dem.tif", max_dist = 5000
# )
# The function checks for each time if the sun elevation (90°−zenith) is greater than the horizon elevation at that azimuth.
# If not, in_sun is FALSE and the irradiance is set to zero.




