#Sensitivity analyses of varying h0=c(2:60)
#h0 = 2:60 (van Rijn, 2011; Khojasteh et al., 2021)

#_________________________________________________
# CURRENTLY WORKING ESTUARY: Weser 
#_________________________________________________

# Import or install libraries needed. If installing is needed, please uncomment line 3-4 and install those libraries first.
# install.packages("tidyr")
# install.packages("ggplot2")
library(tidyr)
library(ggplot2)
library(readxl)

# # Set working directory
# setwd("~/OneDrive - University of Glasgow/Post-doc/GALLANT/Fairhurst_Model/CURRENT DESIGN MODEL/129384_M02_008/results/analyses")

# Import Table_S1_data
# Table_S1_data <- read_excel("/Volumes/Macintosh HD/Users/oprasojo/OneDrive - Postdoc/OneDrive - University of Glasgow/Post-doc/GALLANT/Paper/Table_S1.xlsx")

# Identify the row for the currently working estuary (e.g., Weser)
current_estuary <- "Weser"  # Change this to the actual estuary name in your Table S1
estuary_row <- Table_S1_data[Table_S1_data$Name == current_estuary, ]

# Extract constants for the currently working estuary
H <- estuary_row$`H [m]`
b0 <- estuary_row$`b0 [m]`
h0 <- estuary_row$`h0 [ m ]`
Lb <- estuary_row$`Lb [ m ]`
Lh <- estuary_row$`Lh[m]`
SLR <- 0.84
g <- 9.81
u <- estuary_row$`u [m/s]`
F <- estuary_row$`F [-]`
xmax <- -(estuary_row$`x [m]`)
#x shall have negative value upstream of the estuary mouth

# Variables [m]
# Replace the value 'by' by any value to do sensitivity of the resolution of variables
x <- seq(0, xmax, by = -100)
# h0 <- seq(2, 60, length.out = abs(xmax/100) + 1)
h0 <- seq(2, 60, by = 0.1)

# Calculate b 
b <- b0 * exp(x / Lb)

# Calculate h for every h0 value
h <- sapply(1:length(h0), 
            function(i) h0[i] * exp(x / Lh))

# Calculate db/dx
dbdx <- data.frame((b0/Lb)*exp(x/Lb))

# Create an empty matrix to store db/dx values
dhdx <- matrix(nrow = length(x), ncol = length(h0))

# Calculate db/dx for each b0 value
for (i in 1:length(h0)) {
  dhdx[, i] <- -(h0[i] / Lh) * exp(x / Lh)
}

# Make them as a table
table_data <- data.frame(x,b,h,dbdx,dhdx)
colnames(table_data) <- c("x", 
                          "b",
                          paste0("h", 1:length(h0)), 
                          "dbdx",
                          paste0("dhdx", 1:length(h0)))

# Calculate Term 1 values and add new columns to table_data
table_data$Term1 <- 0.5 * H * table_data$dbdx / table_data$b

# Calculate Term 2 values and add new columns to table_data
for (i in 1:length(h0)) {
  column_name <- paste0("dhdx", i)
  t_column <- 0.5 * H * table_data[[column_name]] / table_data[[paste0("h", i)]]
  new_column_name <- paste0("Term2_", i)
  table_data[[new_column_name]] <- t_column
}

# Calculate Term 2 values and add new columns to table_data with SLR
for (i in 1:length(h0)) {
  column_name <- paste0("dhdx", i)
  t_column <- 0.5 * H * table_data[[column_name]] / (table_data[[paste0("h", i)]] + SLR)
  new_column_name <- paste0("Term2_SLR_", i)
  table_data[[new_column_name]] <- t_column
}

# Calculate Term 3
# F is assumed to be constant since we want to do sensitivity analyses of estuary geometry
# F is gathered from Khojasteh et al., 2021 for each estuary
# Other assumed constants
# pi = 3.141593
# g = 9.81 m2/s
# theta = 52.5 deg (Leuven et al., 2023)z

# Calculate Term 3 values and add new columns to table_data
for (i in 1:length(h0)) {
  t_column <- (-F*(u^2))/(3*pi*9.81*table_data[[paste0("h", i)]]* cos(52.5*(pi/180)))
  new_column_name <- paste0("Term3_", i)
  table_data[[new_column_name]] <- t_column
}

# Calculate Term 3 values and add new columns to table_data with SLR
for (i in 1:length(h0)) {
  t_column <- (-F*(u^2))/(3*pi*9.81*(table_data[[paste0("h", i)]] + SLR)* cos(52.5*(pi/180)))
  new_column_name <- paste0("Term3_SLR_", i)
  table_data[[new_column_name]] <- t_column
}

# Calculate sum of Term 1, 2 and 3
for (i in 1:length(h0)) {
  sum_col <- paste0("Sum", i)
  term2_col <- paste0("Term2_", i)
  term3_col <- paste0("Term3_", i)
  table_data[[sum_col]] <- table_data$Term1 + table_data[[term2_col]] + table_data[[term3_col]] + H
}

# Calculate sum of Term 1, 2 and 3
for (i in 1:length(h0)) {
  sum_col <- paste0("Sum_SLR", i)
  term2_col <- paste0("Term2_SLR_", i)
  term3_col <- paste0("Term3_SLR_", i)
  table_data[[sum_col]] <- table_data$Term1 + table_data[[term2_col]] + table_data[[term3_col]] + H
}

# Plotting x vs Sum1, 2, 3 and so on using ggplot
# to plot multiple Sum columns against x, it's better to reshape the data from wide format to long format.

# Create a new data frame for plotting
plot_data <- table_data[, c("x", 
                            paste0("Sum", 1:length(h0)))]
plot_data_SLR <- table_data[, c("x",
                                paste0("Sum_SLR", 1:length(h0)))]

# Convert to long data format
plot_data_long <- gather(plot_data, key = "Sum_column", value = "Sum_value", -x)
plot_data_long$Sum_column <- sub("Sum", "", plot_data_long$Sum_column)
plot_data_long$Sum_column <- as.numeric(plot_data_long$Sum_column)

plot_data_long_SLR <- gather(plot_data_SLR, key = "Sum_column", value = "Sum_value", -x)
plot_data_long_SLR$Sum_column <- sub("Sum_SLR", "", plot_data_long_SLR$Sum_column)
plot_data_long_SLR$Sum_column <- as.numeric(plot_data_long_SLR$Sum_column)

# Rename second column's name as h0
colnames(plot_data_long)[2]  <- "h0"
colnames(plot_data_long_SLR)[2]  <- "h0"

# Calculate the normalised tidal range
# Normalised tidal range = new tidal range affected by intervention/initial tidal range
plot_data_long$TRnorm <- plot_data_long$Sum_value/H
plot_data_long_SLR$TRnorm <- plot_data_long_SLR$Sum_value/H

# Plot with SLR
# ggplot() +
#   geom_tile(data = plot_data_long, aes(x = x, y = TRnorm, color = h0)) +
#   geom_tile(data = plot_data_long_SLR, aes(x = x, y = TRnorm, color = h0)) +
#   labs(x = "Distance from the river mouth, x [m]") +
#   ylab(expression(paste(italic(η), "/", italic(η[0]), " [-]"))) +
#   scale_colour_gradient(low="#D3DCF2", high="#48577D") +
#   theme_linedraw() +
#   theme(axis.ticks.length=unit(-0.25, "cm"),
#         axis.text = element_text(size = 15),
#         text = element_text(size = 15),
#         panel.border = element_rect(colour = "black", fill=NA, size=1),
#         legend.position = "none") +
#   geom_line(data = subset(plot_data_long, h0 ==  63), aes(x = x, y = TRnorm, group = h0), color = "red") #Put h0 value from TableS1 for this estuary

# 
# # Plot without SLR
# ggplot(plot_data_long, aes(x = x, y = Sum_value, color = h0)) +
#   geom_tile() +
#   labs(x = "Distance from the river mouth, x [m]", 
#        y = "Tidal range, TR [m]") +
#   scale_colour_gradient(low="#D3DCF2", high="#48577D") +
#   scale_y_continuous(limits = c(2.99970,3.0002)) +
#   theme_classic() +
#   geom_line(data = subset(plot_data_long, h0 ==  9), aes(group = h0), color = "red") #Put h0 value from TableS1 for this estuary

# Exporting data as dataframe instead of csv files
# Weser_h0 <- plot_data_long 
# Weser_h0_SLR <- plot_data_long_SLR  

# Export results (interventions with and without SLR) as csv tables to be further analysed
write.csv(plot_data_long, "Weser_h0.csv", row.names=FALSE)
write.csv(plot_data_long_SLR, "Weser_h0_SLR.csv", row.names=FALSE)