rm(list=ls())
library(tidyverse)
library(ggpubr)
source("fun.R")

# species vectors for plotting:
cols <- c("AI" = "#FDE725FF", "KO" = "#5DC863FF","AC" = "#21908CFF",  "AV" = "#3B528BFF", "SA" = "#440154FF")
speciesnames <- c("AI" = "Acanthus ilicifolius", "KO" = "Kandelia obovata", "AC" = "Aegiceras corniculatum", "AV" = "Avicennia marina", "SA" = "Sonneratia apetala")

# load data
df.branch <- read.csv("dfbranch.csv")
df.leaf <- read.csv("dfleaf.csv")
df.drag <- read.csv("dfdrag.csv")
df.size <- read.csv("dfModelledDimensions.csv")

HATO_ms <- 0.5661871 #max velocity during typhoon hato
RHO = 1000 #fluid density, kg/m3

### --- get parameter values for branch of 1 cm diameter -----
diameters <- c(1) #cm
u_ms <- seq(0,1.5,0.01) #m/s

df.MOR_spp <- df.branch %>%
  group_by(species) %>%
  drop_na(meanMOR_MPa) %>%
  summarise(meanMOR = mean(meanMOR_MPa))

df.Cd_spp <- df.drag %>%
  group_by(species, leaves) %>%
  drop_na(meanCd) %>%
  summarise(meanCd = mean(meanCd))

### --- combine parameters value into one dataset -----
df.hat <-
  full_join(df.MOR_spp, df.Cd_spp, by = "species") %>%
  full_join(df.size, by = c("species", "leaves")) 

df.hat <- df.hat %>% slice(rep(1:n(), each = length(u_ms)))
df.hat <- df.hat %>% mutate(u_ms = rep(u_ms, times = nrow(df.hat)/length(u_ms)))

df.hat$species<-factor(df.hat$species, levels = c("AI", "KO", "AC", "AV","SA"))

### --- calculate branch strength -----
df.calc <- df.hat %>%
  mutate(
    R = diam_cm/100/2,
    L = length_m, #* Lfac,
    Fd_N = (0.5 * RHO * area_m2 * meanCd * u_ms^2),
    Fd_N = (0.5 * RHO * area_m2 * meanCd * u_ms^2),
    Fmax_N = ((meanMOR * 1000000 * pi * R^3)/L))

### --- calculate branch breakage point (at which u_ms) -----
df.breakingpoint <- 
  df.calc %>%
  group_by(species, leaves, diam_cm) %>% 
  filter(Fd_N > Fmax_N) %>% 
  select(species, leaves, diam_cm, break_ms = u_ms, break_N = Fmax_N) %>% 
  filter(row_number() == 1)

### --- plot estimated drag per species, with and without leaves -----
p.storm_wl <- 
  ggplot(filter(df.calc, leaves == "with"), 
         aes(x = u_ms, y = Fd_N, group = species)) +
  geom_line(aes(colour = species), size = 1) +
  geom_point(data = filter(df.breakingpoint, leaves == "with"), 
             aes(x = break_ms, y = break_N, colour = species), size = 3) +
  geom_hline(aes(yintercept = Fmax_N, colour = species), linetype = "dashed") +
  geom_vline(xintercept = HATO_ms, colour = "red", linetype = 2) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  theme_classic() +
  labs(x = expression(paste("Orbital velocity u (m ", s^-1,")")), 
       y = expression(paste("Estimated drag force ", hat(F)[D], " (N)")),
       colour = "Species",
       title = "(a) Branches with leaves") +
  scale_x_continuous(limits=c(0, 1)) + scale_y_continuous(limits = c(0,100)) +
  theme(legend.position="bottom") +
  guides(col = guide_legend(ncol = 3))

p.storm_nl <- 
  ggplot(filter(df.calc, leaves == "without"), 
         aes(x = u_ms, y = Fd_N, group = species)) +
  geom_line(aes(colour = species), size = 1) +
  geom_point(data = filter(df.breakingpoint, leaves == "without"), 
             aes(x = break_ms, y = break_N, colour = species), size = 3) +
  geom_hline(aes(yintercept = Fmax_N, colour = species), linetype = "dashed") +
  geom_vline(xintercept = HATO_ms, colour = "red", linetype = 2) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  theme_classic() +
  labs(x = expression(paste("Orbital velocity u (m ", s^-1,")")), 
       y = expression(paste("Estimated drag force ", hat(F)[D], " (N)")),
       colour = "Species",
       title = "(b) Branches without leaves") +
  xlim(0,1) +
  theme(legend.position="bottom") 

ggarrange(p.storm_wl, p.storm_nl, ncol = 2, common.legend = TRUE, legend = "bottom")