# Carregar os dados

#data <- read.xlsx('Data/rrs.xlsx', detectDates = T)
#estacoes = read.xlsx('Data/stations.xlsx', detectDates = T)

QWIP = function(data) {
  
  
  
rrsT <- data[,-1]

# Extrair Rrs de 400 a 700 nm (colunas 52 a 352 no MATLAB, R usa índice baseado em 1)
Rrs_vis <- rrsT[, 1:301]
wave <- 400:700

# Extrair bandas específicas
Rrs_492 <- rrsT$Rrs_492
Rrs_665 <- rrsT$Rrs_665

# Calcular AVW
m <- nrow(Rrs_vis)
AVW <- numeric(m)
for (i in 1:m) {
  AVW[i] <- sum(Rrs_vis[i, ]) / sum(Rrs_vis[i, ] / wave)
}

# Calcular NDI
index_492 <- which.min(abs(wave - 492))
index_665 <- which.min(abs(wave - 665))
NDI <- (Rrs_vis[, index_665] - Rrs_vis[, index_492]) / 
  (Rrs_vis[, index_665] + Rrs_vis[, index_492])

# Polinômio para ajuste
p <- c(-8.399884740300151e-09, 1.715532100780679e-05, -1.301670056641901e-02,
       4.357837742180596e+00, -5.449532021524279e+02)
avw_poly <- 400:640

# Função equivalente à polyval do MATLAB
polyval <- function(p, x) {
  y <- rep(0, length(x))
  n <- length(p)
  for (i in 1:n) {
    y <- y + p[i] * x^(n - i)
  }
  return(y)
}

fit1 <- polyval(p, avw_poly)  # criar função abaixo


# Previsão de NDI e cálculo de QWIP
NDI_pred <- polyval(p, AVW)
QWIP_score <- NDI - NDI_pred
abs_QWIP_score <- abs(QWIP_score)
QWIP_flag <- abs_QWIP_score >= 0.2

# Classificação do tipo de água
Rrs_665b <- Rrs_vis[, 266]
Rrs_560b <- Rrs_vis[, 161]
Rrs_492b <- Rrs_vis[, 93]

Step1 <- Rrs_665b > Rrs_560b
Step2 <- Rrs_665b > 0.025
Step3 <- Rrs_560b < Rrs_492b

ind_600A <- Step1 | Step2
ind_500A <- !Step1 & !Step2 & !Step3
ind_400A <- !Step1 & !Step2 & Step3

# Plotando
library(ggplot2)

df_plot <- data.frame(
  AVW = AVW,
  NDI = NDI,
  class = factor(
    ifelse(ind_600A, "600A",
           ifelse(ind_500A, "500A",
                  ifelse(ind_400A, "400A", NA)))
  )
)

# Geração das curvas de tolerância
fit_df <- data.frame(
  avw = avw_poly,
  fit = fit1,
  fit1a = fit1 + 0.1,
  fit1b = fit1 - 0.1,
  fit2a = fit1 + 0.2,
  fit2b = fit1 - 0.2,
  fit3a = fit1 + 0.3,
  fit3b = fit1 - 0.3,
  fit4a = fit1 + 0.4,
  fit4b = fit1 - 0.4
)

# Plot principal com ggplot2
a = ggplot(df_plot, aes(x = AVW, y = NDI, color = class)) +
  geom_point(size = 0.7, alpha = 0.6) +
  geom_line(data = fit_df, aes(x = avw, y = fit), color = "black", size = 1.2) +
  geom_line(data = fit_df, aes(x = avw, y = fit1a), linetype = "dashed", color = "green") +
  geom_line(data = fit_df, aes(x = avw, y = fit1b), linetype = "dashed", color = "green") +
  geom_line(data = fit_df, aes(x = avw, y = fit2a), linetype = "dashed", color = "#E6B800") +
  geom_line(data = fit_df, aes(x = avw, y = fit2b), linetype = "dashed", color = "#E6B800") +
  geom_line(data = fit_df, aes(x = avw, y = fit3a), linetype = "dashed", color = "#D95319") +
  geom_line(data = fit_df, aes(x = avw, y = fit3b), linetype = "dashed", color = "#D95319") +
  geom_line(data = fit_df, aes(x = avw, y = fit4a), color = "red", size = 1) +
  geom_line(data = fit_df, aes(x = avw, y = fit4b), color = "red", size = 1) +
  scale_color_manual(values = c("400A" = "blue", "500A" = "green", "600A" = "red")) +
  labs(
    x = "AVW (nm)",
    y = paste0("NDI (", wave[index_492], ",", wave[index_665], ")"),
    color = "Classe"
  ) +
  xlim(440, 630) +
  ylim(-2.5, 2) +
  theme_minimal(base_size = 14)

ggplotly(a)



estacoes$QWIP = QWIP_score
#
#filter(estacoes, abs(QWIP) > 0.2) %>% dim()


return(estacoes)

}