datos <- data.frame(
  
  Fecha = as.Date(c(
    "2026-09-17",
    "2026-09-18",
    "2026-09-19",
    "2026-09-20",
    "2026-09-21"
  )),
  
  Tmax = c(18, 19, 18, 20, 19),
  
  Tmin = c(8, 9, 8, 9, 8),
  
  HR_max = c(92, 90, 91, 88, 90),
  
  HR_min = c(55, 53, 57, 50, 54),
  
  Viento = c(1.5, 1.7, 1.4, 1.8, 1.6),
  
  Radiacion = c(14.2, 15.0, 13.8, 16.1, 14.7),
  
  Latitud = c(4.71, 4.71, 4.71, 4.71, 4.71),
  
  Altitud = c(2640, 2640, 2640, 2640, 2640),
  
  Kc = c(0.85, 0.85, 0.85, 0.85, 0.85),
  
  Eficiencia = c(0.90, 0.90, 0.90, 0.90, 0.90),
  
  Precipitacion = c(0, 0, 2, 0, 0)
)

datos
##        Fecha Tmax Tmin HR_max HR_min Viento Radiacion Latitud Altitud   Kc
## 1 2026-09-17   18    8     92     55    1.5      14.2    4.71    2640 0.85
## 2 2026-09-18   19    9     90     53    1.7      15.0    4.71    2640 0.85
## 3 2026-09-19   18    8     91     57    1.4      13.8    4.71    2640 0.85
## 4 2026-09-20   20    9     88     50    1.8      16.1    4.71    2640 0.85
## 5 2026-09-21   19    8     90     54    1.6      14.7    4.71    2640 0.85
##   Eficiencia Precipitacion
## 1        0.9             0
## 2        0.9             0
## 3        0.9             2
## 4        0.9             0
## 5        0.9             0
datos <- datos %>%
  mutate(
    
    es_Tmax = 0.6108 *
      exp((17.27 * Tmax) /
            (Tmax + 237.3)),
    
    es_Tmin = 0.6108 *
      exp((17.27 * Tmin) /
            (Tmin + 237.3))
    
  )

datos
##        Fecha Tmax Tmin HR_max HR_min Viento Radiacion Latitud Altitud   Kc
## 1 2026-09-17   18    8     92     55    1.5      14.2    4.71    2640 0.85
## 2 2026-09-18   19    9     90     53    1.7      15.0    4.71    2640 0.85
## 3 2026-09-19   18    8     91     57    1.4      13.8    4.71    2640 0.85
## 4 2026-09-20   20    9     88     50    1.8      16.1    4.71    2640 0.85
## 5 2026-09-21   19    8     90     54    1.6      14.7    4.71    2640 0.85
##   Eficiencia Precipitacion  es_Tmax  es_Tmin
## 1        0.9             0 2.063989 1.072769
## 2        0.9             0 2.197393 1.148060
## 3        0.9             2 2.063989 1.072769
## 4        0.9             0 2.338281 1.148060
## 5        0.9             0 2.197393 1.072769
datos <- datos %>%
  mutate(
    Tmedia = (Tmax + Tmin) / 2
  )
datos <- datos %>%
  mutate(
    
    es = (es_Tmax + es_Tmin) / 2,
    
    ea = (
      es_Tmin * HR_max / 100 +
      es_Tmax * HR_min / 100
    ) / 2
    
  )
datos <- datos %>%
  mutate(
    VPD = es - ea
  )
datos <- datos %>%
  mutate(
    
    Delta =
      4098 *
      (
        0.6108 *
          exp(
            (17.27 * Tmedia) /
              (Tmedia + 237.3)
          )
      ) /
      (Tmedia + 237.3)^2
    
  )
datos <- datos %>%
  mutate(
    
    P =
      101.3 *
      (
        (293 - 0.0065 * Altitud) /
          293
      )^5.26
    
  )
datos <- datos %>%
  mutate(
    Gamma = 0.000665 * P
  )
calcular_Ra <- function(fecha, latitud){

  J <- yday(fecha)

  dr <- 1 +
    0.033 *
    cos(
      (2*pi/365) * J
    )

  delta_solar <-
    0.409 *
    sin(
      (2*pi/365) * J - 1.39
    )

  phi <- latitud * pi / 180

  ws <-
    acos(
      -tan(phi) *
        tan(delta_solar)
    )

  Ra <-
    ((24 * 60) / pi) *
    0.0820 *
    dr *
    (
      ws *
        sin(phi) *
        sin(delta_solar)
      +
        cos(phi) *
        cos(delta_solar) *
        sin(ws)
    )

  return(Ra)
}

datos$Ra <- mapply(
  calcular_Ra,
  datos$Fecha,
  datos$Latitud
)
datos <- datos %>%
  mutate(
    
    Rns =
      (1 - 0.23) *
      Radiacion
    
  )
datos <- datos %>%
  mutate(
    
    Rso =
      (
        0.75 +
          2e-5 * Altitud
      ) * Ra,
    
    Rnl =
      4.903e-9 *
      (
        (
          Tmax + 273.16
        )^4 +
          (
            Tmin + 273.16
          )^4
      ) / 2 *
      (
        0.34 -
          0.14 *
          sqrt(ea)
      ) *
      (
        1.35 *
          pmin(
            Radiacion / Rso,
            1
          ) -
          0.35
      )
    
  )
datos <- datos %>%
  mutate(
    Rn = Rns - Rnl
  )
datos <- datos %>%
  mutate(
    
    ETo =
      (
        0.408 *
          Delta *
          Rn +
          Gamma *
          (
            900 /
              (Tmedia + 273)
          ) *
          Viento *
          (es - ea)
      ) /
      (
        Delta +
          Gamma *
          (
            1 +
              0.34 *
              Viento
          )
      )
    
  )
datos <- datos %>%
  mutate(
    ETc = ETo * Kc
  )
datos <- datos %>%
  mutate(
    
    Necesidad_neta =
      pmax(
        ETc - Precipitacion,
        0
      )
    
  )
datos <- datos %>%
  mutate(
    
    Necesidad_bruta =
      Necesidad_neta /
      Eficiencia
    
  )
resultados <- datos %>%
  select(
    Fecha,
    Tmax,
    Tmin,
    HR_max,
    HR_min,
    Viento,
    Radiacion,
    Tmedia,
    es,
    ea,
    VPD,
    Delta,
    P,
    Gamma,
    Ra,
    Rn,
    ETo,
    Kc,
    ETc,
    Precipitacion,
    Necesidad_neta,
    Necesidad_bruta
  )

resultados
##        Fecha Tmax Tmin HR_max HR_min Viento Radiacion Tmedia       es       ea
## 1 2026-09-17   18    8     92     55    1.5      14.2   13.0 1.568379 1.061071
## 2 2026-09-18   19    9     90     53    1.7      15.0   14.0 1.672727 1.098936
## 3 2026-09-19   18    8     91     57    1.4      13.8   13.0 1.568379 1.076347
## 4 2026-09-20   20    9     88     50    1.8      16.1   14.5 1.743171 1.089717
## 5 2026-09-21   19    8     90     54    1.6      14.7   13.5 1.635081 1.076042
##         VPD      Delta        P      Gamma       Ra       Rn      ETo   Kc
## 1 0.5073083 0.09797057 73.74675 0.04904159 37.26951 9.059272 2.787726 0.85
## 2 0.5737905 0.10373567 73.74675 0.04904159 37.26148 9.439806 3.034089 0.85
## 3 0.4920323 0.09797057 73.74675 0.04904159 37.25164 8.875034 2.706451 0.85
## 4 0.6534539 0.10672477 73.74675 0.04904159 37.23999 9.937963 3.301275 0.85
## 5 0.5590389 0.10081806 73.74675 0.04904159 37.22654 9.289453 2.945020 0.85
##        ETc Precipitacion Necesidad_neta Necesidad_bruta
## 1 2.369567             0       2.369567        2.632852
## 2 2.578976             0       2.578976        2.865529
## 3 2.300483             2       0.300483        0.333870
## 4 2.806084             0       2.806084        3.117871
## 5 2.503267             0       2.503267        2.781407
resumen <- data.frame(

  Indicador = c(
    "ETo promedio",
    "ETc promedio",
    "ETo acumulada",
    "ETc acumulada",
    "Necesidad neta acumulada",
    "Necesidad bruta acumulada"
  ),

  Valor = c(
    mean(datos$ETo, na.rm = TRUE),
    mean(datos$ETc, na.rm = TRUE),
    sum(datos$ETo, na.rm = TRUE),
    sum(datos$ETc, na.rm = TRUE),
    sum(datos$Necesidad_neta, na.rm = TRUE),
    sum(datos$Necesidad_bruta, na.rm = TRUE)
  ),

  Unidad = c(
    "mm/día",
    "mm/día",
    "mm",
    "mm",
    "mm",
    "mm"
  )

)

resumen
##                   Indicador     Valor Unidad
## 1              ETo promedio  2.954912 mm/día
## 2              ETc promedio  2.511675 mm/día
## 3             ETo acumulada 14.774561     mm
## 4             ETc acumulada 12.558377     mm
## 5  Necesidad neta acumulada 10.558377     mm
## 6 Necesidad bruta acumulada 11.731530     mm
ggplot(
  datos,
  aes(x = Fecha)
) +

  geom_line(
    aes(
      y = ETo,
      linetype = "ETo"
    ),
    linewidth = 1
  ) +

  geom_line(
    aes(
      y = ETc,
      linetype = "ETc"
    ),
    linewidth = 1
  ) +

  labs(
    title = "Evapotranspiración de referencia y del cultivo",
    x = "Fecha",
    y = "Evapotranspiración (mm/día)",
    linetype = "Variable"
  ) +

  theme_minimal()