library(psychrolib)
SetUnitSystem("SI")


# Step 1: Create a dataframe with hourly temperature data for a full year

# Let's assume the dataframe is named "temp_data" and has two columns: "datetime" (with hourly timestamps) and "temperature" (with hourly temperature readings).
df <- read.csv("data.csv",col.names=c("DateTime", "DB", "RH","DP","P","WIND","WIND-DIR","SUN","RAIN","START"))

df=df[complete.cases(df), ]
df$DB=(df$DB/10)
df$RH=(df$RH/100)
df$P=(df$P*100)
df$HR=GetHumRatioFromRelHum(df$DB, df$RH, df$P)
df$H=GetMoistAirEnthalpy(df$DB, df$HR)
# create a sample data frame

#print(df)
# use the aggregate function to calculate the average of y values for each x value
df=aggregate(H ~ DB, data = df, FUN = mean)



# Step 1: Create a dataframe with hourly temperature and relative humidity data for a full year
# Let's assume the dataframe is named "df" and has three columns: "datetime" (with hourly timestamps), "DB" (with hourly dry bulb temperature readings), and "RH" (with hourly relative humidity readings).
# Step 2: Calculate the total number of hours in a year
# Step 2: Calculate the total number of hours in a year
total_hours <- nrow(df)

print(total_hours)
nrow(unique(df[,c("DB","H")]))
# Step 3: Create a new dataframe named "temp_pct" with four columns: "DB", "RH", "pct_exceeded", and "RH_pct_exceeded"
temp_pct <- data.frame(DB = df$DB, H = df$H, pct_exceeded = numeric(length(df$DB)))

# Step 4: For each temperature value in the "DB" column of the "df" dataframe, calculate the percentage of time that temperature is exceeded and the relative humidity for the first occurrence of that temperature, and store them in the corresponding rows of the "temp_pct" dataframe
for (i in seq_along(df$H)) {
  temp <- df$H[i]
  hours_exceeded <- sum(df$H > temp)
  pct_exceeded <- hours_exceeded / total_hours * 100
  rh <- df$DB[df$H == temp][1]
  temp_pct$pct_exceeded[i] <- pct_exceeded
#  temp_pct$RH_pct_exceeded[i] <- rh
}

# Step 5: Calculate the average relative humidity for each unique value of "pct_exceeded" and store it in the "RH_pct_exceeded" column of the "temp_pct" dataframe
#for (pct in unique(temp_pct$pct_exceeded)) {
#  avg_rh <- mean(temp_pct$RH_pct_exceeded[temp_pct$pct_exceeded == pct])
#  temp_pct$RH_pct_exceeded[temp_pct$pct_exceeded == pct] <- avg_rh
#}

# Step 6: Order the "temp_pct" dataframe by the "pct_exceeded" column, with lowest percentage first
temp_pct <- temp_pct[order(temp_pct$pct_exceeded), ]

# Step 7: Print the first 20 rows of the "temp_pct" dataframe
temp_pct$HR=GetHumRatioFromEnthalpyAndTDryBulb(temp_pct$H, temp_pct$DB)

head(temp_pct,5)
temp_pct=head(temp_pct, 5)


library(ggplot2)
library(ggpsychro)


p = ggpsychro(tdb_lim = c(-15, 40)) +
    geom_grid_relhum(alpha = 1.0, label.alpha = 1.0, label.size = 6, label.fontface = 2) +
    scale_relhum(minor_breaks = NULL) +
    geom_grid_wetbulb(size = 1.0, color = "black", alpha = 1.0, label_loc = NA) +
    scale_wetbulb(breaks = seq(25, 30, by = 5), minor_breaks = NULL) +
    geom_grid_vappres(label.size = 5) +
    scale_vappres(breaks = seq(6000, 7000, by = 500), limits = c(6000, 7000)) +
    geom_grid_specvol() +
    scale_specvol(labels = NULL) +
#    geom_line(aes(x = 20:30,  = 64000, y=1,color="orange"), stat = "wetbulb") +
    geom_point(aes(temp_pct$DB, relhum = GetRelHumFromHumRatio(temp_pct$DB, temp_pct$HR, 101300),color="orange"), stat = "relhum")+
#    geom_line(aes(x = 20:30, enthalpy =64000,y=1), stat = "enthalpy")+
#    geom_point(aes(x=25, enthalpy=75,y=1),stat = "enthalpy")+
    geom_grid_enthalpy()
ggsave("chart5.pdf", plot = p, width = 6, height = 4)


# Step 8: Create a list of unique combinations of "RH_pct_exceeded" and "pct_exceeded"
#combo_list <- unique(temp_pct[, c("RH_pct_exceeded", "pct_exceeded")])