# Sensitivity analyses of varying Lh=c(5000:300000)=Lb

#_________________________________________________
# CURRENTLY WORKING ESTUARY: Clyde
#_________________________________________________

# 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("/Volumes/Samsung_T5/analyses_v2")

# 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., Clyde)
current_estuary <- "Clyde"  # 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
Lh <- seq(5000, 300000, by = 100)
x <- seq(0, xmax, by = -100)

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

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

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

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

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

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

# 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(Lh)) {
  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(Lh)) {
  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)

#Calculate Term 3 values and add new columns to table_data
for (i in 1:length(Lh)) {
  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 with SLR and add new columns to table_data
for (i in 1:length(Lh)) {
  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(Lh)) {
  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 with SLR
for (i in 1:length(Lh)) {
  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(Lh)))]

plot_data_SLR <- table_data[, c("x",
                                paste0("Sum_SLR", 1:length(Lh)))]

# 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)

# Calculate the total number of repetitions to replace column Sum_column as Lh values as x resolution
num_repetitions <- length(Lh) * length(x)

# Repeat each value of Lh length(x) times
repeated_Lh <- rep(Lh, each = length(x))

# Update the second column of the table with the repeated values
plot_data_long[, 2] <- repeated_Lh[1:num_repetitions]
plot_data_long_SLR[, 2] <- repeated_Lh[1:num_repetitions]

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

# 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 = Lh)) +
#   geom_tile(data = plot_data_long_SLR, aes(x = x, y = TRnorm, color = Lh),size = 2) +
#   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)) +
#   geom_line(data = subset(plot_data_long, Lh ==  17000), aes(x = x, y = TRnorm, group = Lh), color = "red") #Put Lh value from TableS1 for this estuary

# 
# # Plot without SLR
# ggplot(plot_data_long, aes(x = x, y = Sum_value, color = Lh)) +
#   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, Lh ==  66300), aes(group = Lh), color = "red") #Put Lh value from TableS1 for this estuary

# Exporting data as dataframe instead of csv files
# Clyde_Lh <- plot_data_long 
# Clyde_Lh_SLR <- plot_data_long_SLR  

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