## load packages for analyses
#load required packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(readxl, ggplot2)

## load bifurcation runs for PCLake output
df.tpsummer=read_excel('data/rawdata_fig_5.xlsx')
#first remove range 0.001-0.005 from first model runs (as you have additional more detailed runs on them)
sub.tpsummer=df.tpsummer[-c(1:3, 40:42),]

## figure 5a  -----
## plot tp summer bifurcation
plot.summer=ggplot(sub.tpsummer, aes(x=ParamBifurValue, y=AuxilBifurValue*1000))+
  geom_rect(aes(xmin=0.0013, xmax=0.0028, ymin=0, ymax=1456), fill='lightgray')+
  geom_rect(aes(xmin=0, xmax=0.0365, ymin=4.4, ymax=1456), fill=NA, col='goldenrod1', linewidth=1)+
  geom_line()+
  geom_point(shape=21, aes(fill=as.factor(RunId)), size=4)+
  scale_fill_manual(values=c('black', 'white'))+
  theme_classic()+
  theme(axis.text = element_text(size=20),
        axis.title = element_text(size=20),
        legend.position = 'none')+
  labs(x=expression(paste('Nutrient load (g P m'^-2,' d'^-1,')')), 
       y=expression(paste('TP concentration (', mu,'g L'^-1,')')),fill='Initial state')

#save the plot
ggsave('scratch/fig5a.png', plot.summer, width=8, height=6)


#figure 5b ----
fig.zoom=ggplot(sub.tpsummer, aes(x=ParamBifurValue, y=AuxilBifurValue*1000))+
  geom_rect(aes(xmin=0.0013, xmax=0.0028, ymin=0, ymax=250), fill='lightgray')+
  geom_rect(aes(xmin=0, xmax=0.00337, ymin=4.4, ymax=250), fill=NA, col='red4', linewidth=1)+
  xlim(0,0.006)+
  ylim(0,250)+
  geom_point(shape=21, aes(fill=as.factor(RunId)), size=4)+
  scale_fill_manual(values=c('black', 'white'))+
  geom_segment(aes(x = 0.00157, y = 57.311, xend = 0.0051, yend = 238.89), linewidth=1)+
  geom_segment(aes(x = 0.0001, y = 0, xend = 0.00255, yend = 27.054), linewidth=1)+
  theme_classic()+
  theme(axis.text = element_text(size=20),
        axis.title = element_text(size=20),
        legend.position = 'none')+
  labs(x=expression(paste('Nutrient load (g P m'^-2,' d'^-1,')')), 
       y=expression(paste('TP concentration (', mu,'g L'^-1,')')),fill='Initial state')

ggsave('scratch/fig5b.png', fig.zoom, width=8, height=6)