

## Vegetation emissivity function ---------
# Questa funzione gestisce diversi casi:
# . input vettoriali di Lmis (W/m²) e di temperature (puoi passare T_surface oppure T_air);
# . opzionale L_down (radiation down) e parametri tau, L_path;
# . ritorna vettori di emissività e incertezza ε stimata (se fornite incertezza LW_OUT e T)

# Incertezza dovuta all'uso di Tair anzichè Tsurf:
# Una possibilità per migliorare la stima con T_air, è stimare T_surface come: 
# Ts ≈ Tair + ΔT
# dove ΔT può essere:
# . costante (es. +1 K di giorno, 0 K di notte);
# . funzione empirica di irradianza shortwave incidente e vento (se hai i dati);
# . o una fascia di incertezza: Tair ± δT.
# In questa funzione usiamo un parametro opzionale deltaT_surface che applica una correzione a T_air quando T_surface non è fornita. 




# Funzione per stimare emissività


surf_emissivity <- function(
    L_up_meas,           # W/m2   (misura upwelling)
    T_surface = NULL,    # K o NULL
    T_air = NULL,        # K (usato se T_surface NULL)
    deltaT_surface = 0,  # correzione da applicare a T_air se T_surface mancante
    simplified = FALSE,  # se TRUE (sensore vicino a surface) usa la funzione semplice ε = L↑ / σTs^4, se FALSE (sensore misura a distanza) 
    L_down = 0,          # W/m2   (downwelling incident sulla superficie), scalare o vettore
    tau = 1,             # trasmittanza atmosfera (0-1), scalar o vettore
    L_path = 0,          # W/m2   emissione lungo la linea di vista (scalar o vect)
    sigma_L = NA,        # errore su L_up_meas (W/m2), scalar o vect, opzionale
    sigma_T = NA         # errore su T (K), scalar o vect, opzionale
) {

  stefan <- 5.670374419e-8
  
  # riproduci scalari a lunghezza n
  n <- length(L_up_meas)
  vecrep <- function(x) if(length(x)==1) rep(x, n) else x
  L_up_meas <- vecrep(L_up_meas)
  L_down <- vecrep(L_down)
  tau <- vecrep(tau)
  L_path <- vecrep(L_path)
  deltaT_surface <- vecrep(deltaT_surface)
  # gestione sigma (incertezza sensori LW e AT)
  sigma_L <- if(all(is.na(sigma_L))) rep(NA, n) else vecrep(sigma_L)
  sigma_T <- if(all(is.na(sigma_T))) rep(NA, n) else vecrep(sigma_T)
  
  # scegliere T da usare
  if (!is.null(T_surface)) {
    T_use <- vecrep(T_surface)
    T_label <- "T_surface (fornita)"
  } else if (!is.null(T_air)) {
    T_use <- vecrep(T_air) + deltaT_surface
    # T_label <- paste0("T_air + deltaT (", unique(deltaT_surface), " K)")
    if (any(deltaT_surface == 0)) {
      warning("Stai usando T_air come proxy per T_surface. Considera incertezza maggiore.")
    }
  } else {
    stop("Devi fornire almeno T_surface o T_air (K).")
  }

  
  # Calcolo epsilon (emissivity) ........
  # simplified: ε = L↑ / σ * Ts^4
  # rigourous:  ε = L↑ - τ * L↓ - Lpath / τ(σ * Ts^4 - L↓)
  
  epsilon <- numeric(n)
  epsilon_unc <- rep(NA, n)
  T_label <- rep(NA, n)
  
  i <- 1
  for (i in seq_len(n)) {
    
    if (simplified == TRUE) { 
      # denominatore (NON tenendo conto L_down)
      denom <- stefan * T_use[i]^4
    } else {
      # denominatore (tenendo conto L_down)
      denom <- tau[i] * (stefan * T_use[i]^4 - L_down[i])
    }
    
    if (abs(denom) < .Machine$double.eps) {
      epsilon[i] <- NA
      next
    }
    
    if (simplified == TRUE) {
      # formula generale (senza atmosfera)
      epsilon[i] <- L_up_meas[i] / denom
    } else {
      # formula generale (con atmosfera)
      epsilon[i] <- (L_up_meas[i] - tau[i] * L_down[i] - L_path[i]) / denom
    }
    
    # limita [0,1] (ma non forzare se vuoi analizzare outliers)
    epsilon[i] <- pmin(pmax(epsilon[i], 0), 1)
    
    # incertezza approssimativa (se forniti sigma_L e sigma_T)
    if (!is.na(sigma_L[i]) && !is.na(sigma_T[i]) && !is.na(epsilon[i])) {
      # derivazione semplificata: eps = (L - tau Ld - Lp) / (tau sigma T^4) = N / D
      # d(eps)/eps ≈ dN/N - dD/D ; dominant terms -> approx:
      relL <- ifelse(L_up_meas[i] - tau[i]*L_down[i] - L_path[i] != 0,
                     sigma_L[i] / abs(L_up_meas[i] - tau[i]*L_down[i] - L_path[i]), 0)
      relT <- (4 * sigma_T[i] / T_use[i])
      epsilon_unc[i] <- abs(epsilon[i]) * sqrt(relL^2 + relT^2)
    } else {
      epsilon_unc[i] <- NA
    }
    
    T_label[i] <- paste0("T_air + deltaT (", deltaT_surface[i], " K)")
  }
  
  data.frame(
    L_up_meas = L_up_meas,
    T_used = T_use,
    L_down = L_down,
    tau = tau,
    L_path = L_path,
    emissivita = epsilon,
    emissivita_sigma = epsilon_unc,
    T_label = T_label,
    stringsAsFactors = FALSE
  )
  
}





# ==============================================================================.
# FUNZIONI GEOMETRIE -----------------------------------------------------------
# APPROCCIO: Il traliccio blocca (disturba) un SETTORE del cono FoV del radiometro
# ==============================================================================.



library(dplyr)
library(tibble)


# ------------------------------------------------------------------------------.
# Funzione per calcolare l'angolo solido di un settore rettangolare di un cono con apertura half_angle
# VERSIONE VETTORIZZATA
# ------------------------------------------------------------------------------.

omega_sector <- function(theta_w,         # theta_w: apertura angolare orizzontale (azimutale) del settore [rad]
                         theta_h_min,     # theta_h_min, theta_h_max: angoli zenitali min e max del settore [rad]
                         theta_h_max, 
                         half_angle_fov   # half_angle_fov: semi-apertura del cono FoV [rad]
) {
  
  # Per un settore di un cono, l'angolo solido è:
  # Ω = ∫∫ sin(θ) dθ dφ = Δφ × [cos(θ_min) - cos(θ_max)]
  # dove:
  # φ è l'angolo azimutale
  # Δφ = θ_w è l'apertura azimutale (larghezza angolare del traliccio)
  # θ_min = θ_top` è l'angolo zenitale verso la cima del traliccio
  # θ_max = half_ang` è il limite del FoV
  
  # Limitiamo l'integrazione al range del traliccio: usa pmax/pmin per vettorizzazione corretta
  theta_min <- pmax(0, theta_h_min)
  theta_max <- pmin(half_angle_fov, theta_h_max)
  # theta_min <- max(0, theta_h_min)
  # theta_max <- min(half_angle_fov, theta_h_max)
  
  # Integrale: Ω = Δφ * ∫[θ_min, θ_max] sin(θ) dθ
  #           = Δφ * [-cos(θ)]_θ_min^θ_max
  #           = Δφ * (cos(θ_min) - cos(θ_max))
  
  # Usa ifelse per condizione vettorizzata
  # delta_phi == theta_w
  omega <- ifelse(
    theta_max <= theta_min,
    0,
    theta_w * (cos(theta_min) - cos(theta_max))
  )
  
  return(omega)
}


# ------------------------------------------------------------------------------.
# Funzione principale per calcolo disturbo e stats
# ------------------------------------------------------------------------------.

disturbo_radiometro <- function(
    fov_deg,        # Field of view del radiometro (gradi)
    h_radiometro,   # Altezza sensore dal suolo (m)
    w_traliccio,    # Larghezza traliccio (m)
    d_braccio,      # Distanza orizzontale radiometro-traliccio (m)
    h_c,            # Vettore altezza copertura vegetale (m)
    albedo,         # Albedo del traliccio (0-1)
    irradianza      # Vettore irradianza incidente (W/m²)
) {
  
  stopifnot(length(h_c) == length(irradianza))
  stopifnot(fov_deg > 0 && fov_deg <= 180)
  stopifnot(h_radiometro > 0)
  stopifnot(d_braccio > 0)
  stopifnot(w_traliccio > 0)
  
  # Conversioni angolari
  fov_rad <- fov_deg * pi/180
  half_ang <- fov_rad/2
  
  # Angolo solido TOTALE del FoV (cono circolare)
  omega_fov <- 2*pi*(1 - cos(half_ang))
  
  # Angolo azimutale sotteso dalla larghezza del traliccio (visto dall'alto)
  theta_w <- 2 * atan(w_traliccio / (2 * d_braccio))
  
  # Tibble con parametri output disturbo
  tibble(
    
    h_c = h_c,
    irradianza = irradianza
    
  ) %>%
    
    mutate(
      
      # Altezza utile del traliccio (sopra la vegetazione)
      h_utile = pmax(0, h_radiometro - h_c),
      
      # Angoli zenitali dei punti estremi del traliccio visibile 
      # theta = 0 è la verticale verso il basso
      # theta aumenta allontanandosi dalla verticale
      
      # 1. Angolo zenitale verso il top del traliccio (alla quota h_c, cima vegetazione)
      theta_top = atan(d_braccio / pmax(h_utile, 0.01)),  # evita divisione per zero
      
      # 2. Angolo zenitale verso il bottom del traliccio visibile
      # Il traliccio si estende verso il basso fino al limite del FoV
      # Se il traliccio è molto alto (h_utile grande), il bottom potrebbe essere fuori dal FoV
      # L'altezza massima visibile nel FoV è: h_max_visible = d_braccio / tan(half_ang)
      # h_max_visible = d_braccio / tan(half_ang)
      # Se h_utile > h_max_visible, il bottom del traliccio è al limite del FoV (theta = half_ang)
      # Altrimenti, il bottom del traliccio è "all'infinito" verso il basso → theta = half_ang
      # In pratica, per tralicci verticali, theta_bottom è sempre al limite del FoV
      theta_bottom = if_else(
        h_utile > 0,
        atan(d_braccio / 0.01),  # praticamente infinito → va fino al limite del FoV
        0
      ),
      theta_bottom = pmin(theta_bottom, half_ang),  # limitato dal FoV (in pratica sempre)
      
      # Verifica se il traliccio è visibile nel FoV
      # È visibile se theta_top < half_ang (il top del traliccio è dentro il cono)
      visibile = (h_utile > 0) & (theta_top < half_ang),
      
      # Calcolo dell'angolo solido occupato dal traliccio
      # Il traliccio blocca un settore del cono da theta_top a theta_bottom con apertura azimutale theta_w
      angolo_solido_traliccio = if_else(
        visibile,
        omega_sector(theta_w, theta_top, theta_bottom, half_ang),
        0
      ),
      
      # Frazione del FoV occupata dal traliccio
      frazione_disturbo = angolo_solido_traliccio / omega_fov,
      
      # Area proiettata del traliccio vista dal radiometro
      # Approssimazione: w_traliccio * h_visibile
      h_visibile = if_else(
        visibile,
        h_utile - d_braccio * tan(theta_top),  # altezza del traliccio visibile
        0
      ),
      h_visibile = pmax(0, h_visibile),
      
      area_proiettata_m2 = w_traliccio * h_visibile,
      
      # Potenza parassita riflessa dal traliccio
      # Usando il modello Lambertiano: L = (E * albedo) / π
      # Potenza per unità di area del sensore: P/A = L * Ω
      potenza_parassita_Wm2 = (irradianza * albedo / pi) * angolo_solido_traliccio,
      
      
      # Output info geometriche
      angolo_solido_FoV = omega_fov,
      theta_w_deg = theta_w * 180/pi,
      theta_top_deg = theta_top * 180/pi,
      theta_bottom_deg = theta_bottom * 180/pi
      
    ) %>%
    
    select(
      angolo_solido_FoV,
      angolo_solido_traliccio,
      frazione_disturbo,
      area_proiettata_m2,
      potenza_parassita_Wm2,
      h_utile,
      h_visibile,
      theta_w_deg,
      theta_top_deg,
      theta_bottom_deg,
      visibile,
      h_c,
      irradianza
    )
}





# # # ESEMPIO CON I DATI  -----------------
# # cat("=== ANALISI DISTURBO RADIOMETRO (VERSIONE CORRETTA) ===\n\n")
# # 
# fov_deg <- 170
# h_radiometro <- 39.5
# w_traliccio <- 1.93
# d_braccio <- 2.04
# h_c_val <- c(1,10,38)
# albedo <- 0.2
# irradianza <- rep(265.302, length(h_c_val))
# 
# # Calcolo
# risultati <- disturbo_radiometro(
#   fov_deg = fov_deg,
#   h_radiometro = h_radiometro,
#   w_traliccio = w_traliccio,
#   d_braccio = d_braccio,
#   # h_c = rep(h_c_val, length(irradianza)),
#   h_c = h_c_val,
#   albedo = albedo,
#   irradianza = irradianza
# )
# 
# # cat("RISULTATI:\n")
# # print(risultati, n = 6)
# #
# # cat("\n\nSTATISTICHE:\n")
# cat(sprintf("Frazione di disturbo: %.4f%% (%.1f ppm)\n",
#             risultati$frazione_disturbo[1]*100,
#             risultati$frazione_disturbo[1]*1e6))
# # cat(sprintf("Rapporto: 1 parte su %.0f\n", 1/risultati$frazione_disturbo[1]))
