# Sensitivity analyses of varying n=c(0.02:0.1)

#_________________________________________________
# CURRENTLY WORKING ESTUARY: Limpopo
#_________________________________________________

# 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., Limpopo)
current_estuary <- "Limpopo"  # 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) #dx resolution is 100 m
n <- seq(0.02, 0.1, by = 0.001)

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

# Calculate h 
h <- h0 * exp(x/Lh)

# Calculate db/dx to be independent of dx
dbdx <- data.frame((b0/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 <- cbind(x, b,h,dbdx,dhdx)

# Rename the column names
colnames(table_data) <- c("x", 
                          "b",
                          "h", 
                          "dbdx",
                          "dhdx")

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

# Create an empty matrix to store C values
C <- matrix(nrow = length(x), ncol = length(n))

# Calculate C for each n value and h value
for (i in 1:length(n)) {
  C[, i] <- (1/n[i])*(h^(1/6))
}

# Calculate F for each C value
F <- matrix(nrow = length(x), ncol = length(n))
for (i in 1:length(n)) {
  F[, i] <- 8*g/(C[, i]^2)
}

# Calculate Term 3
Term3 <- matrix(nrow = length(x), ncol = length(n))
for (i in 1:length(n)) {
  Term3[,i] <- (-F[, i] * (u^2)) / (3*pi*g*table_data$h*cos(52.5*(pi/180)))
}

# Calculate Term 3 SLR
Term3_SLR <- matrix(nrow = length(x), ncol = length(n))
for (i in 1:length(n)) {
  Term3_SLR[,i] <- (-F[, i] * (u^2)) / (3*pi*g*(table_data$h + SLR)*cos(52.5*(pi/180)))
}

# Integrate them to table_data
table_data$Term3 <- Term3
table_data$Term3_SLR <- Term3_SLR

# Initialize a list to store Sum values
Sum_list <- list()

# Calculate Sum columns for each Term3 column
for (i in 1:length(n)) {
  Sum_list[[i]] <- table_data$Term1 + table_data$Term2 + table_data$Term3[, i] + H
  col_name <- paste0("Sum", i)
  table_data[[col_name]] <- Sum_list[[i]]
}

# Calculate Sum columns for each Term3_SLR column
for (i in 1:length(n)) {
  Sum_list[[i]] <- table_data$Term1 + table_data$Term2 + table_data$Term3_SLR[, i] + H
  col_name <- paste0("Sum_SLR", i)
  table_data[[col_name]] <- Sum_list[[i]]
}

# 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(n)))]
plot_data_SLR <- table_data[, c("x",
                                paste0("Sum_SLR", 1:length(n)))]

#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 as n
colnames(plot_data_long)[2]  <- "n"
colnames(plot_data_long_SLR)[2]  <- "n"

# Calculate the normalised tidal range
plot_data_long$TRnorm <- plot_data_long$Sum_value/H
plot_data_long_SLR$TRnorm <- plot_data_long_SLR$Sum_value/H

ggplot(plot_data_long, aes(x = x, y = TRnorm, color = n)) + 
  geom_line() +
  geom_line(data = plot_data_long_SLR, color = "brown", alpha = 0.2) +
  labs(x = "Distance from the river mouth, x [m]") +
  ylab(expression(paste(italic(η), "/", italic(η[0]), " [-]"))) +
  scale_color_gradient(low="#D3DCF2", high="#48577D") +
  scale_fill_brewer() + 
  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, n ==  4), aes(x = x, y = TRnorm, group = n), color = "yellow") #Put n value from TableS1 for this estuary

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

# Exporting data as dataframe instead of csv files
# Limpopo_n <- plot_data_long 
# Limpopo_n_SLR <- plot_data_long_SLR  

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