
# ..................................................................................
# Solar zenith angle, sun elevation angle, daytime/nighttime splitting
# ..................................................................................

# lo zenith angle θz è l’angolo tra la verticale locale e la direzione del Sole.
# la solar elevation (altezza del Sole) è h=90°−θz
# se h>0° -> Sole sopra l’orizzonte (giorno, geometricamente).
# se h≤0° -> Sole sotto l’orizzonte (notte, geometricamente).
# Per osservazioni pratiche si usano spesso soglie “apparente” che tengono conto di rifrazione atmosferica e raggio solare visibile:
# alba/tramonto apparente: h=−0.833° (includendo rifrazione + semidiametro solare).
# Classi di crepuscolo (utile per certe analisi):
# Civil twilight: −6° ≤ h ≤ 0°
# Nautical twilight:  −12° ≤ h ≤ 6°
# Astronomical twilight:  −8° ≤ h ≤ 12°
# Night (astronomica):  h < 18°

library(dplyr)

# gradi in radianti
deg2rad <- function(deg) {return(deg * pi / 180)}

# radianti in gradi
rad2deg <- function(rad) {return(rad * 180 / pi)}


# . Declinazione solare (in radianti) 
# J = DoY
solar_declination <- function(J) {
  delta <- (0.409 * sin(2 * pi / 365 * J - 1.39))
  return(delta)
}


# Solar hour angle in radians for given time (decimal hours) and longitude
solar_hour_angle <- function(solar_time_decimal) {
  # solar time decimal hours from 0 to 24 with solar noon at 12
  hour_angle <- (deg2rad(15 * (solar_time_decimal - 12)))
  return(hour_angle)
}


# . Solar zenith angle ----------------------------------------------------------

solar_zenith <- function(latitude,           # latitudine (deg)
                         J,                  # DOY
                         solar_time_decimal, # HoD as decimal 0-24
                         sza.out)            # angle in 'degrees' or 'radians'
{
  
  phi <- deg2rad(latitude)
  delta <- solar_declination(J)
  h <- solar_hour_angle(solar_time_decimal)
  
  # solar elevation angle:
  cos_theta_s <- sin(phi) * sin(delta) + cos(phi) * cos(delta) * cos(h)  # Eq. A2 in Leuning et al., 2012
  
  # Limitare cos_theta_s tra -1 e 1 per sicurezza numerica
  # cos_theta_s <- max(min(cos_theta_s, 1), -1) # original function (return only one value in case of a vec)
  if (any(cos_theta_s < -1)) { cos_theta_s[cos_theta_s < -1] <- -1}
  if (any(cos_theta_s > 1)) { cos_theta_s[cos_theta_s > 1] <- 1}
  
  # zenith angle
  theta_s <- acos(cos_theta_s)
  
  if (sza.out == 'radians') {
    return(theta_s) # in radianti
  } else {
    return(rad2deg(theta_s)) # in gradi
  }
  
}

# . daytime/nighttime/twilight -------------------------------------------------
# theta_z: vettore angoli zenith in radianti
# datetime: vettore POSIXct (opzionale, utile per aggregazioni per giorno)
# method: "geometric" (h>0) o "apparent" (usa -0.833 deg per alba/tramonto apparente)
#          per definire “notte” in modo conservativo per analisi radiometriche, usa method = "apparent" (incluso l’effetto della rifrazione)
# returns: data.frame con elevation_deg e flags

classify_daynight <- function(theta_z, datetime = NULL, method = c("geometric", "apparent")) {
  
  method <- match.arg(method)
  # conversioni
  zen_deg <- theta_z * 180 / pi
  elev_deg <- 90 - zen_deg
  
  # correzione apparente (rifrazione + semidiametro)
  apparent_offset <- ifelse(method == "apparent", -0.833, 0)  # degrees
  
  # flags
  is_day_geom  <- elev_deg > 0
  is_day_appar <- elev_deg > apparent_offset
  
  # twilight classes (based on elevation)
  twilight_class <- rep(NA_character_, length(elev_deg))
  twilight_class[elev_deg > 0] <- "day"
  twilight_class[elev_deg <= 0 & elev_deg > -6] <- "civil_twilight"
  twilight_class[elev_deg <= -6 & elev_deg > -12] <- "nautical_twilight"
  twilight_class[elev_deg <= -12 & elev_deg > -18] <- "astronomical_twilight"
  twilight_class[elev_deg <= -18] <- "night"
  
  if (!is.null(datetime)) {
    out <- tibble(datetime = as.POSIXct(datetime), theta_z = theta_z,
                  zen_deg = zen_deg, elev_deg = elev_deg,
                  is_day_geom = is_day_geom, is_day_appar = is_day_appar,
                  twilight = twilight_class)
  } else {
    out <- tibble(theta_z = theta_z, zen_deg = zen_deg, elev_deg = elev_deg,
                  is_day_geom = is_day_geom, is_day_appar = is_day_appar,
                  twilight = twilight_class)
  }
  return(out)
}
