# Sensitivity analyses of varying b0=c(100:65700)
# b0 = 100:65700 (550m in Clyde, 1000 m from Kojhasteh et al., 2021, 25000 m from van Rijn, 2011, Glasgow b0 = 3000 m, b0 = 65700 from Outer Bay of Fundy)

#_________________________________________________
# 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("~/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., 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
b0 <- seq(100, 65700, by = 100)
x <- seq(0, xmax, by = -100) #dx resolution is 100 m

# Calculate b for every b0 value
b <- sapply(1:length(b0), 
            function(i) b0[i] * exp(x / Lb))

# Calculate h for every b0 and Lb value
h <- h0 * exp(x / Lh)

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

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

# Calculate dh/dx to be independent of dx
dhdx <- data.frame(-(h0/Lh)*exp(x/Lh))

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

# Calculate Term 1 values and add new columns to table_data
for (i in 1:length(b0)) {
  column_name <- paste0("dbdx", i)
  t_column <- 0.5 * H * table_data[[column_name]] / table_data[[paste0("b", i)]]
  new_column_name <- paste0("Term1_", i)
  table_data[[new_column_name]] <- t_column
}

# Calculate Term 2 values and add new columns to table_data
table_data$Term2 <- 0.5 * H * table_data$dhdx / table_data$h

# Calculate Term 2 values and add new columns to table_data with SLR
table_data$Term2_SLR <- 0.5 * H * table_data$dhdx / (table_data$h + SLR)

# Calculate Term 3
# Other assumed constants
# 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
table_data$Term3 <- -F*(u^2)/(3*pi*9.81*table_data$h* cos(52.5*(pi/180)))

# Calculate Term 3 values and add new columns to table_data with SLR
table_data$Term3_SLR <- -F*(u^2)/(3*pi*9.81*(table_data$h + SLR)* cos(52.5*(pi/180)))

# Calculate sum of Term 1, 2 and 3
for (i in 1:length(b0)) {
  sum_col <- paste0("Sum", i)
  Term1_col <- paste0("Term1_", i)
  table_data[[sum_col]] <- table_data[[Term1_col]] + table_data$Term2 + table_data$Term3 + H
}

# Calculate sum of Term 1, 2 and 3 with SLR
for (i in 1:length(b0)) {
  sum_col <- paste0("Sum_SLR", i)
  Term1_col <- paste0("Term1_", i)
  table_data[[sum_col]] <- table_data[[Term1_col]] + table_data$Term2_SLR + table_data$Term3_SLR + 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(b0)))]
plot_data_SLR <- table_data[, c("x",
                                paste0("Sum_SLR", 1:length(b0)))]

# 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 b0 values as x resolution
num_repetitions <- length(b0) * length(x)

# Repeat each value of Lb length(x) times
repeated_b0 <- rep(b0, each = length(x))

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

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

# 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 = b0)) +
#   geom_tile(data = plot_data_long_SLR, aes(x = x, y = TRnorm, color = b0)) +
#   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, b0 ==  3000), aes(x = x, y = TRnorm, group = b0), color = "red") #Put b0 value from TableS1 for this estuary

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

# Exporting data as dataframe instead of csv files
# Clyde_b0 <- plot_data_long 
# Clyde_b0_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_b0.csv", row.names=FALSE)
write.csv(plot_data_long_SLR, "Clyde_b0_SLR.csv", row.names=FALSE)