#####################################################################
# R Script to conduct sediment yield magnitude-frequency analysis
# and plots for each site
#####################################################################

# Run Import_Save_Flow_Data1.R to import and format all of the flow files
rm(list = ls())
library(zoo) # rollmean to calc "mids"
#library(Hmisc) # HD Quantile estimator for AFDC
library(MASS) # RLM
source("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/minor.ticks.axis.R")
# this is for susp sed site only, change site file and wd to run on bedload sites.
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload")

# set number of bins to use
nbins <- 25

# Create Plot?
plot = F

# create trapezoidal rule function
trapz <- function (x, y) {
  idx = 2:length(x)
  return (as.double( (x[idx] - x[idx-1]) %*% (y[idx] + y[idx-1])) / 2)
}

# Function: opposite of %in% or match
'%!in%' <- function(x,y)!('%in%'(x,y)) 

# import flow file data
# Somehow I ended up with some sites in this file that have too few sed measurements...
path <- "site_summary.csv"
flow.sites <- read.table(path, header = TRUE, sep = ",", 
                    stringsAsFactors = F)
flow.sites <- flow.sites[order(flow.sites$site.no),]

# import sed regression data table and match with most up to date list of sites from flow.site file
sed.data <- read.table("results_bed_regression.csv", header = T, stringsAsFactors = F, sep = ",")
# sed.data <- subset(sed.data, sed.data$siteID %in% flow.sites$site.no)
sed.data <- sed.data[order(sed.data$site.no),]
peak.table <- read.csv("results_peaks.csv", header = T, stringsAsFactors = F)
peak.table <- peak.table[order(peak.table$site.no),]

# # culled list from bed being too sandy from USGS GSD data 
# cull.list <- c(11254000,13336500,13120500,09415000)
# 
# keep.rows <- which(flow.sites$site.no %!in% cull.list)
# flow.sites <- flow.sites[keep.rows,]
# peak.table <- peak.table[keep.rows,]
# sed.data <- sed.data[keep.rows,]

sum(sum(flow.sites$site.no == sed.data$site.no) == nrow(flow.sites),
    sum(flow.sites$site.no == peak.table$site.no) == nrow(flow.sites)) == 2

cnames <- c("site.no", "site.nm","Qbf", "Qbf.yield", "RL", "beta", "Qeff.b", "Qeff.e", "Qeff.l","Qeff.RI", 
            "Q50.RI", "Qbf.RI","Qeff.ptile", "Qeff.P", "Q50.P","f.plus", "Qeff.yield", "Qeff.25", "Qeff.75", "Qeff.spread", 
            "Q.25", "Q.50", "Q.75", "yield.spread", "Q25.ptile","Q50.ptile","Q75.ptile",
            "Q.25.l", "Q.50.l", "Q.75.l", "yield.spread.l", 
            "Qeff.25.l", "Qeff.yield.l", "Qeff.75.l", "Qeff.spread.l")
metrics <- data.frame(matrix(data = 0, nrow = nrow(flow.sites), ncol = length(cnames)),
                      row.names = NULL, stringsAsFactors = F)
names(metrics) <- cnames  # Assigns col names to data frame

## Create empty data frame for comparing Qeff methods and assign col names to it
cnames <- c("site.no", "site.nm","Qeff.b", "Qeff.e", "Qeff.l","area.e", "RL", "CV", 
            "mean", "median", "beta")
eff.compare <- data.frame(matrix(data = 0, nrow = nrow(flow.sites), ncol = length(cnames)),
                      row.names = NULL, stringsAsFactors = F)
names(eff.compare) <- cnames  # Assigns col names to data frame

######################################################################
## Create and print (.PDF) sediment yield histograms and calc metrics
######################################################################
if(plot) {
  pdf("Figures/BedSedYieldCurves4.pdf", width = 8.5, height = 11, paper = "special") # turn on PDF creator
  # set plot to be square and 4x4 plots
  par(mfrow = c(2,2), pty = "s", oma = c(3,3,3,3), mar = c(4,5,2,2))
  options0 <- options()
  options(scipen = -1)
}

for ( i in 1:nrow(flow.sites) ) { # nrow(flow.sites)
  #import individual flow data files - note, placing leading zero for some below
  path <- paste("DV_Flows/", flow.sites$site.no[i], ".csv", sep = "")
  flow.data <- read.table(path, header = T, stringsAsFactors = F, sep = ",")
  
  flow0 <- flow.data$val*0.028316847 # Convert to cms
  min.reg <- sed.data$min.flow[i]
  max.reg <- sed.data$max.flow[i]
  Qbf <- flow.sites$Qbf[i]
  Q1.5 <- flow.sites$Q1.5[i]
  
  #Import Peak flow series (annual maximum) if we have a USGS gage for the site
  if(!is.na(as.numeric(flow.sites$site.no[i]))){
    path.pk <- paste("peaks/", flow.sites$site.no[i], ".csv", sep = "")
    peaks <- read.csv(path.pk, header = T, stringsAsFactors = F)
    peaks <- peaks$peak.cfs*0.028316847
  } else if(flow.sites$yrs[i] >= 5){
    qmat <- matrix(data = flow0, nrow = 365, ncol = length(flow0)/365)
    peaks <- apply(qmat, MARGIN = 2, max)
  } else{ peaks <- NA}
  
  # remove flows below the minimum measured bedload value
  #flow <- flow[which(flow >= min.reg - 0.1*min.reg)]
  #flow.mat <- flow # some flow files have NA values for some reason but we need to keep them divisible by 365 for the AFDC calc
  #flow <- flow[!is.na(flow)] # remove NA's
  
  ## Remove flow below threshold, if one exists
  if (!is.na(sed.data$Q0[i])) { 
    flow <- flow0[which(flow0 > sed.data$Q0[i])] 
  } else{ flow <- flow0}
  
  max.flow <- max(flow, na.rm = T)
  min.flow <- min(flow,  na.rm = T)
  cv <- round( sd(flow) / mean(flow), 2 )
  median <- round( median(flow), 2 )
  avg <- round( mean(flow), 2 )
  # I add 1 to max.flow to make sure all flows are accounted for
#   breaks <- seq(0, max.flow + 1, by = (max.flow + 1) / nbins) 
  # midpoint of histogram bins
#   mids <- seq(from = max.flow / (nbins*2), by = (max.flow + 1) / nbins, length.out = nbins)
  nyears <- floor(length(flow)/365)
  
  ###########################################################################################
  # Calc Qeff using histograms
  ##########################################################################################
  
  # create histogram object but suppress plot
  bins <- seq(signif(floor(min(flow)),1), ceiling(max(flow)), length.out = nbins+1 )
  hist <- hist(flow, breaks = bins, plot = FALSE) 
  # extract density from histo object (frequency/n)
  flow.freq <- hist$density
  breaks <- hist$breaks
  # mid value of bin
  mids <- rollmean(hist$breaks,2)
#   # extract bin values (either side of mids)
  breaks <- hist$breaks
  nbreaks <- length(breaks)
  # bin width
  width <- breaks[2] - breaks[1]
  
  # extract rating curve params from "sed.data" table and create rating curve
  b0 <- sed.data$b0_RLM[i]
  b1 <- sed.data$b1_RLM[i]
  BCF <- sed.data$BCF_qmle_RLM[i]
  qs <- b0*BCF*mids^b1

  # create effectiveness histogram
  eff.b <- qs*flow.freq
  # First rough calc of Qeff using bin midpoint
  Qeff.b <- mids[which.max(eff.b)]
  # calc difference between Qeff.b and bin values to ID bin containing Qeff
  bin.diff = abs(Qeff.b - breaks)
  # Returns the values of the bins that stradle Qeff (smallest difference)
  Qeff.bins <- sort(breaks[sort(bin.diff, index.return = T)$ix[1:2]])
  # subsets the flows on record within the bins and calcs the median of these values
  # for final Qeff.b calculation
  Qeff.b <- round(median(subset(flow, Qeff.bins[1] <= flow & flow <= Qeff.bins[2])),3)
  #Qeff.b1 <- Qeff.b

#   # Check if Qeff.b is in first bin, and if so, divide first two bins up until peak not in first bin
#   # Note that first two original bins are divided into 3, then 5, and so on until we reach 11.
#   if (Qeff.b <= breaks[2]){
#   new.bins <- 4
#   new.breaks <- breaks
#   
#     while (Qeff.b < new.breaks[2]) {
#       # sub-divide first two bins of the histogram 
#       new.breaks = c( round( seq(from = 0, to = breaks[2], length.out = new.bins), 2),
#                       breaks[3:nbreaks] )
#       # create new Histo with subdivided first two bins
#       hist <- hist(flow, breaks = new.breaks, plot = FALSE)
#       # extract density (frequency/n)
#       flow.freq <- hist$density 
#       # mid value of bin
#       mids <- rollmean(hist$breaks,2) 
#       
#       # re-calc rating curve with new bins
#       qs <- b0*BCF*mids^b1
#       # create effectiveness histogram
#       eff.b <- qs*flow.freq
#       # First rough calc of Qeff using bin midpoint
#       Qeff.b <- mids[which.max(eff.b)]
#       # calc difference between Qeff.b and bin values to ID bin containing Qeff
#       bin.diff = abs(Qeff.b - new.breaks)
#       # Returns the values of the bins that stradle Qeff
#       Qeff.bins <- sort(new.breaks[sort(bin.diff, index.return = T)$ix[1:2]])
#       # subsets the flows on record within the bins and calcs the median of these values
#       Qeff.b <- round(median(subset(flow, Qeff.bins[1] <= flow & flow <= Qeff.bins[2])),3)
#       
#       # continue to subdivide bins until Qeff.b > breaks[2] (greater than the first bin)
#       new.bins = new.bins + 2
#       
#       if(new.bins >= 11) break
#     }
#   }
#   
#   Qeff.b2 <- Qeff.b
  
  ############################################################################
  # Calc empirical flow CDF and then use a defined window to differentiate it
  # based on number of points to differentiate
  ############################################################################
  
  # number of points to differentiate
  npts <- 50 
  
  # create empirical CDF using flow record
  ecdf.prob <- ( 1:length(flow) ) / length(flow)
  flow.sort <- sort(flow)

  # plot(flow.sort, ecdf.prob, type = "s")
  
  # create logarithmic sequence of Q's at which to differentiate
  xpts1 <- exp( seq( log(min(flow) + 0.001), 
                    log(max(flow) - 0.001), length = npts) )

  # create sequence of interpolated eCDF values evaluated at 'xpts1'
  ypts1 <- approx(flow.sort, ecdf.prob, xout = xpts1)$y
  
  # add extreme from sorted flow and eCDF to vectors for differentiation
  xpts1 <- c(min(flow), xpts1, max(flow))
  ypts1 <- c(min(ecdf.prob), ypts1, max(ecdf.prob))
  
  #plot(flow.sort, ecdf.prob, type = "l")
  #points(xpts1,ypts1)
  
  # Numerical differentiation of eCDF.
  # 1) Eval eCDF at logarithmically-spaced intervals using central difference (j+1) - (j-1)
  # 2) Assign central diff to PDF at mid point flow value (j)
  pdf.e1 <- numeric(length(xpts1) - 2)
  n <- length(xpts1) - 1
  q.e1 <- xpts1[2:n]
  for (j in 2:n) { # j starts at 2, but assigns central diff to PDF at j-1
    pdf.e1[j-1] <- ( ypts1[j+1] - ypts1[j-1] ) / ( xpts1[j+1] - xpts1[j-1] )
  }
  plot(q.e, pdf.e, type = "l",log = "y")
  
  # Calc area under ePDF
  area.e <- signif(trapz(q.e1,pdf.e1), 2)
  
  # Effectiveness calc
  eff.e <- b0*BCF*q.e1^b1*pdf.e1
  #plot(q.e1, eff.e, type = "l") 
  Qeff.e <- round(q.e1[which.max(eff.e)], 3)
  
#   ## Calc Qs50 using numerical integration of effectiveness curve
#   idx = 2:length(eff.e)
#   cum.int <- cumsum( ( (q.e1[idx] - q.e1[idx-1]) * 
#                        (eff.e[idx] + eff.e[idx-1]) / 2 ) ) / trapz(q.e1,eff.e)
#   q.mid <- rollmean(q.e1,2)
#   qs50 <- approx(cum.int, q.mid, xout = .5)$y
  

  ####################################################################################
  # Qeff with Lognormal Daily Flow Approx
  ####################################################################################
  ## Calc lognormal Flwo PDF params - (Yevjevich, 2010)
  u.y <- 0.5 * log(avg^2 / (1 + cv^2))
  s.y <- sqrt(log(1 + cv^2))
  min <- ifelse(min.flow == 0, flow.sort[min(which(flow.sort > 0))], min.flow)
  x <- exp( seq(log(min), log(max.flow*2), by = .01))
  x.c <- rollmean(x, 2)
  
  lnorm.PDF <- dlnorm(x, u.y, s.y)
  
  Qeff.l <- round(exp(u.y + s.y^2*(b1 - 1)), 2)
  eff.l <- b0*BCF*x^b1*lnorm.PDF
  
  f.plus.l <- round(1 - pnorm(-sqrt(log(1 + cv^2))), 2)
  area <- trapz(x, eff.l)
  idx = 2:length(x)
  cum.p.l <- cumsum( (x[idx] - x[idx-1]) * (eff.l[idx] + eff.l[idx-1]) / 2 )/area
  
  Q.25.l <- round( approx(cum.p.l, x.c, xout = 0.25)$y, 2)
  Q.50.l <- round( approx(cum.p.l, x.c, xout = 0.50)$y, 2)
  Q.75.l <- round( approx(cum.p.l, x.c, xout = 0.75)$y, 2)
  yield.spread.l <- round( (Q.75.l - Q.25.l)/Q.50.l, 2)
  
  min.25 <- ifelse(1 - f.plus.l - 0.25 > 0, 1 - f.plus.l - 0.25 > 0, min(cum.p.l))
  max.75 <- ifelse(1 - f.plus.l + 0.25 < max(x), 1 - f.plus.l + 0.25, 1)
  Qeff.25.l <- round( approx(cum.p.l, x.c, xout = min.25)$y, 2)
  Qeff.yield.l <- 1 - f.plus.l
  Qeff.75.l <- round( approx(cum.p.l, x.c, xout = max.75)$y, 2)
  Qeff.spread.l <- round((Qeff.75.l - Qeff.25.l) / Qeff.l, 2)
      
  ####################################################################################
  # Calc Yield Metrics
  ####################################################################################
  
  # Calc event-based cumultive sediment transport curve
  yield <- b0*BCF*flow^b1/1000*86400 # tonnes/day  <----------- DBL check if regressions are in tonnes/sec or kg/sec
  yield.sort <- b0*BCF*flow.sort^b1/1000*86400 # tonnes/day
  yield.cum <- cumsum(yield.sort)/sum(yield)
  yield.start <- 1 + min( which(yield.cum > 0) ) # idx of first non-zero value in yield.cum
  
  # Idx of Qeff on ascending sorted flow vector
  Qeff.sort.idx <- which.min(abs(flow.sort- Qeff.e))

  yield.tot <- sum(yield)
  yield.plus <- sum(yield[which(flow > Qeff.e)])
  f.plus <- round( yield.plus / yield.tot, 2)
  Qeff.yield <- 1 - f.plus
  
  if(!is.na(Qbf)) Qbf.yield <- round(approx(flow.sort, yield.cum, xout = Qbf)$y,2) else Qbf.yield = NA
  
  Q.25 <- round( approx(yield.cum, flow.sort, xout = 0.25)$y, 2) # discharge below which 25% sed transported
  Q.50 <- round( approx(yield.cum, flow.sort, xout = 0.5)$y, 2)  # discharge below which 50% sed transported
  Q.75 <- round( approx(yield.cum, flow.sort, xout = 0.75)$y, 2) # discharge below which 75% sed transported

  yield.spread <- round( (Q.75 - Q.25) / Q.50, 2)
  
  # Yield Spread around Qeff
  # Check if less than 25% of sed yield is above Qeff
  if (f.plus <= 0.25) {
    Qeff.75 <- round( max(flow), 2)
  } else{
    Qeff.75 <- round( approx(yield.cum, flow.sort, xout = Qeff.yield + 0.25)$y, 2)
  }

  if (Qeff.yield <= 0.25) {
    Qeff.25 <- round( min(flow), 2)
  } else{
    Qeff.25 <- round( approx(yield.cum, flow.sort, xout = Qeff.yield - 0.25)$y, 2)
  }


  yield.25 <- approx(flow.sort, yield.cum, xout = Qeff.25)$y
  yield.75 <- approx(flow.sort, yield.cum, xout = Qeff.75)$y
  
  Qeff.spread <- round((Qeff.75 - Qeff.25)/Qeff.e, 2)
   
  
#   #### Fplus calc'd from sediment yield curve and not flow record ###################
#   # Calc f.plus from flow record vs histograms vs kernal density function
#   if( Qeff.b > max(mids)) {
#     fplus <- 0
#   }else{
#   fplus <- trapz(mids[which(mids > Qeff.b)], eff.b[which(mids > Qeff.b)])/ 
#     trapz(append(mids,0,after = 0), append(eff.b,0,after = 0))
#   fplus <- trapz(d$x[which(d$x > Qeff.b)], b0*BCF*d$x[which(d$x > Qeff.b)]^b1*d$y[which(d$x > Qeff.b)]) / 
#     trapz(d$x, b0*BCF*d$x^b1*d$y)
#  }
#   yield.int <- round(trapz(d$x, b0*BCF*d$x^b1*d$y),3)
    
  # Calc RI of metrics using either inst. annual max flow series OR max series from daily avg data
  sort <- sort(flow0, decreasing = T)
  rank <- seq(1, length(sort), 1)
  P <- rank / (length(sort) + 1) 
  P.Qeff <- signif( approx(x = sort, y = P, xout = Qeff.e, rule = 2)$y, 3)
  P.Q50 <- signif( approx(x = sort, y = P, xout = Q.50, rule = 2)$y, 3)
  
  if (!is.na(peaks[1])) {
    sort <- sort(peaks, decreasing = T)
    rank <- 1:length(sort)
    RI <- (length(sort) + 1) / rank
    RI.Qeff <- signif( approx(x = sort, y = RI, xout = Qeff.e, rule = 2)$y, 3)
    RI.Q50 <- signif( approx(x = sort, y = RI, xout = Q.50, rule = 2)$y, 3)
    
    if(Qeff.e < min(peaks, na.rm = T) & nyears > 5) RI.Qeff <- 0.99
    if(Q.50 < min(peaks, na.rm = T) & nyears > 5) RI.Q50 <- 0.99
    
    if(!is.na(Qbf)){
      RI.Qbf <- signif( approx(x = sort, y = RI, xout = Qbf, rule = 2)$y, 3)
      if(Qbf < min(peaks, na.rm = T)){RI.Qbf <- 0.99}
    } else{ RI.Qbf <- NA }
    
  } else{
    RI.Qeff <- NA
    RI.Q50 <- NA
    RI.Qbf <- NA
  }
  
#   RI = 1/365 / (1 - pnorm( ( log(Qeff.e) - mean(log(flow[which(flow>0)])) ) 
#                            / sqrt(var(log(flow[which(flow>0)]))) ) )

  #### Calc percent exceedence of some yield metrics##################
  Qeff.ptile <- signif( mean(flow <= Qeff.e), 4)
  Q25.ptile <- signif( mean(flow <= Q.25), 4)
  Q50.ptile <- signif( mean(flow <= Q.50), 4)
  Q75.ptile <- signif( mean(flow <= Q.75), 4)

  ###########################################################################################
  # Plot Stuff
  ###########################################################################################
  if (plot) {
    
    ################################
    ## Plot Sediment Rating Curve
    ################################
    cex.ax = 1
    cex.lab = 1.15
    
    data_pairs <- read.csv(paste("sediment/", sed.data$site.no[i], ".csv", sep = ""))
    
    ymax <- max( c(log10(sed.data$b0_RLM[i]) + 
                     log10(max(data_pairs$cms))*sed.data$b1_RLM[i],
                   log10(max(data_pairs$kg.sec)) ) )
    ymin <- min( c(log10(sed.data$b0_RLM[i]) + 
                     log10(min(data_pairs$cms))*sed.data$b1_RLM[i],
                   log10(min(data_pairs$kg.sec[which(data_pairs$kg.sec > 0)])) ) )
    
    plot(log10(data_pairs$cms), log10(data_pairs$kg.sec), 
         type = "n",
         xlab = bquote(italic(Q) ~~ (m^3/s)), 
         ylab = bquote(italic(Q)[s] ~~ (kg/s)), 
         main = "Sediment Rating Curve",
         axes = F,
         ylim = c(ymin, ymax),
         cex.lab = cex.lab
    )
    # Add blue RLM curve
    curve(log10(sed.data$b0_RLM[i]) + x*sed.data$b1_RLM[i],
          from = log10(min(data_pairs$cms)), 
          to = log10(max(data_pairs$cms)), n = 2, add = T, col = "royalblue", lwd = 2)
    # Add green OLS Curve
    curve(log10(sed.data$b0_OLS[i]) + x*sed.data$b1_OLS[i],
          from = log10(min(data_pairs$cms)), 
          to = log10(max(data_pairs$cms)), n = 2, add = T, col = "seagreen", lwd = 2)
    
    #Re-draw poinrs over line
    points(log10(data_pairs$cms), log10(data_pairs$kg.sec), col = "gray20")
    #   # Add vertical line where threshold exists (if it exists)
    #   if( !is.nan(sed.data$Q0[j])) {
    #     abline(v = sed.data$Q0[j], lty = 2, col = "gray20")
    #   }
    
    xmin <- log10( min(data_pairs$cms) )
    xmax <- log10( max(data_pairs$cms) )
    #   ymin <- log10( min(data_pairs$kg.sec[which(!data_pairs$kg.sec == 0)]) )
    #   ymax <- log10( max(data_pairs$kg.sec) )
    xrange = xmax - xmin 
    yrange = ymax - ymin 
    
    minor.ticks.axis(1,9, mn = floor(xmin), mx = ceiling(xmax), cex = cex.ax)
    minor.ticks.axis(2,9, mn = floor(ymin), mx = ceiling(ymax), cex = cex.ax)
    box()
    
    r2 = sed.data$r2[i]
    rs = sed.data$cor_S[i]
    b0 = format( signif(sed.data$b0_RLM[i],2), scientific = T ) 
    b1 = round(sed.data$b1_RLM[i],2)
    label1 = bquote(italic(R)^2 == .(r2))
    label2 = bquote(italic(R)[s] == .(rs))
    label3 = bquote(italic(Q)[s] ~ "=" ~ .(b0)*italic(Q)^.(b1) )
    
    text(x = (xmin + 0.05*xrange),
         y = (ymax - 0.15*yrange),
         labels = label1,  
         cex = .8, adj = c(0,0))
    text(x = (xmin + 0.05*xrange), 
         y = (ymax - 0.25*yrange),
         labels = label2,  
         cex = .8, adj = c(0,0))
    text( x = (xmin + 0.05*xrange),
          y = (ymax - 0.06*yrange),
          labels = label3,
          cex = .8, adj = c(0,0),
          col = "royalblue")
    
    site.no <- ifelse(nchar(flow.sites$site.no[i]) > 7, flow.sites$site.no[i], 
                      paste(0,flow.sites$site.no[i], sep = ""))
    
    mtext(paste(flow.sites$site.nm[i], ": ", site.no, sep = ""),
          side = 3,
          line = -2,
          cex = 1.25, 
          outer = T,
          font = 2)
    mtext("Bed Load Sites",
          side = 3, 
          line = 0, 
          cex = 1.25, 
          outer = T, 
          font = 2, 
          col = rgb(165/255,42/255,42/255, alpha = 0.75))
    
    legend( "bottomright", #x = xmax*.8, y = ymin - ymin*.05,
            c("RLM", "OLS"),
            lwd = 2,
            col = c("royalblue", "seagreen"),
            bty = "n",
            cex = 0.8,
            seg.len = 2,
            x.intersp = 1,
            y.intersp = 1
    )
    
    ##########################
    ## Plot Flow Frequencies
    ##########################
    # strip zeros for plot positioning purposes
    y.b <- flow.freq[ which(flow.freq != 0) ]
    y.e <- pdf.e1[ which(pdf.e1 != 0) ]
    
    y.min <- 10^floor( log10( min( c(min(y.b), min(y.e)) ) ) )
    y.max <- 10^ceiling( log10(max( c(max(flow.freq), max(pdf.e1)) ) ) )  
    
    plot(mids, flow.freq,
         axes = F,
         log = "y",
         type = "h",
         lend = 2,
         lwd = 6.5,
         xlab = bquote(italic(Q) ~~ (m^3/s)),
         ylab = "Probability Density",
         main = "Flow Density Distribution",
         ylim = c(y.min, y.max),
         cex.lab = cex.lab)
    
    lines(mids, flow.freq,
          type = "h",
          lend = 1,
          lwd = 3.5,
          col = "gray80")
    
    # Plot ePDF
    lines(q.e1, pdf.e1,
          lwd = 3,
          lty = 1,
          col = rgb(165/255,42/255,42/255, alpha = 0.75)
          )
    
    axis(side = 1, cex.axis = cex.ax)
    
    cycles <- log10(y.max) + abs(log10(y.min)) + 1
    ticks <- 10^seq(log10(y.min), log10(y.max), by = 1)
    axis(side = 2, at = ticks, labels = NA, cex.axis = cex.ax)
    
    at <- seq(log10(y.min), log10(y.max), by = 2)
    labels <- sapply(at, function(i) as.expression(bquote(10^ .(i))) )
    axis(side = 2, at = 10^at, labels = labels, cex.axis = cex.ax, tick = F)
    box()
    
    # Add lognormal PDF
    lines(x, lnorm.PDF, lty = 6, lwd = 2)
    
    legend("topright", #x = max(mids)*.35, y = y.max*8,
           legend = c("Histogram","Empirical PDF", "Lognormal PDF"),
           lty=c(1,1,6), # gives the legend appropriate symbols (lines)
           lwd=c(6,3,3),
           col=c("gray80","brown", "black"), # gives the legend lines the correct color and width
           cex = 0.8,
           seg.len = 2.5,
           x.intersp = 1,
           y.intersp = 1,
           bty = "n"
    )
    
    
    ##############################
    ## Plot Effectiveness Curves
    ##############################
    # create space at top of plot for text and legend
    y.max = max( c(max(eff.b), max(eff.e)) ) + 0.3*max( c(max(eff.b), max(eff.e)) )
    
    # Plot eff histo: blue vertical bars
    plot(mids, eff.b, 
         type = "h",
         lend = 2,
         lwd = 6.5,
         col = "black",
         xlab = bquote(italic(Q) ~~ (m^3/s)),
         ylab = "Sediment Yield Density",
         main = "Sediment Yield Curves: Effectiveness",
         ylim = c(0, y.max),
         cex.lab = cex.lab,
         cex.axis = cex.ax
    )
    
    lines(mids, eff.b,
          type = "h",
          lend = 1,
          lwd = 3.4,
          col = "royalblue")
    
    # Plot ePDF: red line
    lines(q.e1, eff.e, lwd = 3, col = rgb(165/255,42/255,42/255, alpha = 0.75))
    
    # Plot Lognormal Effectiveness curve
    lines(x, eff.l, lwd = 2, lty = 6)
    
    # add circles at each Qeff
    points(mids[which.max(eff.b)], max(eff.b), lwd = 2, cex = 1.5, col = "gray50")
    points(Qeff.e, max(eff.e), lwd = 2, cex = 1.5, col = "gray50")
    
    # Add verticals for regression extents and Qbf
    abline(v = min.reg, lty = 2, col = "gray40")
    abline(v = max.reg, lty = 2, col = "gray40")
    abline(v = Qbf, lty = 6)
    
    # for text and legend
    usr0 <- par("usr")       # get user coordinates
    par(usr = c(0, 1, 0, 1)) # new relative user coordinates
    
    #   a = format( signif(sed.data$b0_RLM[i]*sed.data$BCF_qmle_RLM[i],2), scientific = T ) 
    #   b = round(sed.data$b1_RLM[i],2)
    #   
    #   label1 = bquote(italic(Q)[s] ~ "=" ~ .(a)*italic(Q)^.(b))
    label2 = bquote(italic(RL) ~ "=" ~ .(nyears) ~ "yrs")
    
    #   text(x =  0.1,
    #        y = 0.95,
    #        labels = label1,  
    #        cex = .7, 
    #        adj = 0,
    #        col = "gray20")
    legend("topright", #x = 0.3, y = 1.05,
           legend = c("Yield Histogram","Empirical Yield Curve", "Lognormal Yield Curve",
                      "Rating Curve Limits", "Bankfull Discharge"),
           lty=c(1,1,6,2,6), # gives the legend appropriate symbols (lines)
           lwd=c(6,3,2,1,1),
           col=c("royalblue","brown", "black", "gray40", "black"), # gives the legend lines the correct color and width
           cex = 0.7, # changed legend text size
           #bty = "n",
           bg = rgb(1,1,1, alpha = 0.65),
           seg.len = 2.5,
           x.intersp = 1,
           y.intersp = 1,
           box.col = rgb(1,1,1, alpha = 0.65))
    box()
    
    text(x =  0.7,
         y = 0.7,
         labels = label2,  
         cex = .8, 
         adj = 0,
         col = "gray20")
    
    par(usr = usr0) # reset plot coords to original values
    
    ########################
    ## Sediment Yield Plot
    ########################
    plot(flow.sort,yield.cum, 
         type = "n", 
         main = "Cumulative Sediment Yield",
         xlab = bquote(italic(Q) ~~ (m^3/s)), 
         ylab = "Cumulative Sediment Yield",
         cex.axis = cex.ax,
         cex.lab = cex.lab)
    
    # Plot IQR centered around Qeff.yield with line segments
    flow.IQR <- c( seq(Qeff.25, max(flow.sort[which(flow.sort < Qeff.75)]), 
                       by = (Qeff.75 - Qeff.25) / 200), min( c(Qeff.75, max(flow) ) ) )
    yield.IQR <- approx(flow.sort, yield.cum, xout = flow.IQR)$y
    xx <- c(flow.IQR, flow.IQR[length(flow.IQR)], flow.IQR[1])
    yy <- c(yield.IQR, 0, 0)
    
    poly.col <- rgb(165/255,42/255,42/255, alpha = 0.75)
    polygon(x = xx, y = yy, 
            col = poly.col, 
            border = NA,
            lwd = 3)
    
    segments(Qeff.e, 0, Qeff.e, Qeff.yield, lwd = 4, col = "gray40", lend = 1)
    segments(0, Qeff.yield, Qeff.e, Qeff.yield, lwd = 4, col = "gray40", lend = 1)
    segments(0, yield.25, Qeff.25, yield.25, lwd = 2, col = poly.col, lty = 1, lend = 1)
    segments(0, yield.75, Qeff.75, yield.75, lwd = 2, col = poly.col, lty = 1, lend = 1)
    
    # Plot IQR centered around q.50 (discharge assoc. with 50% cum sed transport)
    seg.col <- rgb(red = 100/255, green = 149/255, blue = 237/255, alpha = 0.7)
    segments(Q.50, 0, Q.50, 0.50, lwd = 2, lty = 6, lend = 1, col = "cornflowerblue")
    segments(Q.25, 0, Q.25, 0.25, lwd = 3, lty = 3, col = "cornflowerblue")
    segments(Q.75, 0, Q.75, 0.75, lwd = 3, lty = 3, col = "cornflowerblue")
    
    # plot cum sed yield
    lines(flow.sort, yield.cum,
          lwd = 2.5)
    
    text(0,0.95,
         labels = paste( "Yield @ Qeff = ", round(Qeff.yield, 2)*100, "%", sep = ""),
         cex = 0.8,
         col = "gray20",
         adj = c(0,0))
    
    usr0 <- par("usr")       # get user coordinates
    par(usr = c(0, 1, 0, 1)) # new relative user coordinates
    
    legend("bottomright", #x = 0.32, y = 0.42,
           legend = c("Cumulative Yield",expression(italic(Q[eff])~"Yield"), 
                      expression(italic(Q[eff])~"Yield IQR"),
                      "Half Yield", "Yield IQR"),
           lty=c(1,1,1,6,3), # gives the legend appropriate symbols (lines)
           lwd=c(2,4,6,2,3),
           col=c("black", "gray40","brown", "cornflowerblue", "cornflowerblue"), # gives the legend lines the correct color and width
           cex = 0.7, # changed legend text size
           #bty = "n",
           seg.len = 3,
           x.intersp = 1,
           y.intersp = 1,
           box.col = rgb(1,1,1, alpha = 0.65),
           bg = rgb(1,1,1, alpha = 0.65))
    box()
    par(usr = usr0) # reset plot coords to original values
    
    mtext(bquote( "Page"~.(i)~"of"~.(nrow(flow.sites)) ), side = 1, outer = T, cex = .7)
    
  }
  
  
  ###########################################################################################
  # Write Results
  ###########################################################################################

  # write results to sed yield metrics table
  metrics[i,] <- c(flow.sites$site.no[i], flow.sites$site.nm[i], Qbf, Qbf.yield, nyears, b1, 
                   Qeff.b, Qeff.e, Qeff.l, RI.Qeff, RI.Q50, RI.Qbf, 
                   Qeff.ptile, P.Qeff, P.Q50, f.plus, Qeff.yield,
                   Qeff.25, Qeff.75, Qeff.spread, Q.25, Q.50, Q.75, yield.spread,
                   Q25.ptile, Q50.ptile, Q75.ptile,
                   Q.25.l, Q.50.l, Q.75.l, yield.spread.l, 
                   Qeff.25.l, Qeff.yield.l, Qeff.75.l, Qeff.spread.l)

  # write results to eff.compare table
  eff.compare[i,] <- c(flow.sites$site.no[i], flow.sites$site.nm[i], round(Qeff.b,2), 
                       round(Qeff.e,2), Qeff.l, area.e, nyears, cv, avg, median, b1)
}

if(plot) dev.off()

# embed_fonts("FlowHists_VaryBin.pdf", outfile="FlowHists_VaryBin_embed.pdf")
############################################################################

# save table
write.csv(metrics, file = "results_yield_metrics.csv", row.names = F)
write.csv(eff.compare, file = "results_eff_compare.csv", row.names = F)


######################################################################################
# Compare different methods for calculating Qeff with flow metrics and beta
######################################################################################
rm(list = ls())

# this is for susp sed site only, change site file and wd to run on bedload sites.
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/BedLoad/")

eff.compare <- read.table("results_eff_compare.csv", header = T, stringsAsFactors = F, sep = ",")
reg.results <- read.table("results_bed_regression.csv", header = T, stringsAsFactors = F, sep = ",")

# make sure the records line up
sum(eff.compare$site.no != reg.results$site.no) == 0

# 
# compare <- data.frame(apply(eff.compare[,3:11], 2, as.numeric))
# compare <- compare[-39,]
# compare <- compare[-49,]

# compare <- data.frame(apply(eff.compare[,3:10], 2, as.numeric))
# compare <- compare[-39,]
# compare <- compare[-49,]

# rm(Qeff.e, mean, Qeff.e, Qeff.l, Qeff.b, cv, avg, median)

# plot stuff: 
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload/Figures")
attach(eff.compare)

beta <- as.numeric(reg.results$b1_RLM)

pdf("Qeff_Pairs.pdf", width = 8, height = 6)
  par(mfrow = c(1,3), pty = "s", mar = c(5.1,5.1,4.1,2.1), 
      cex.lab = 1.1, cex.axis = 1, oma = c(0,3,0,0))
  
  plot(Qeff.b, Qeff.l, log = "xy",
       xlab = bquote(italic(Q[eff.b]) ~~~ "(" ~ m^3/s ~ ")"),
       ylab = bquote(italic(Q[eff.l]) ~~~ "(" ~ m^3/s ~ ")"),
       ylim = c(1E-1, 6E3),
       xlim = c(1E-1, 6E3))
  abline(0,1)
  
  plot(Qeff.b, Qeff.e, log = "xy",
       xlab = bquote(italic(Q[eff.b]) ~~~ "(" ~ m^3/s ~ ")"),
       ylab = bquote(italic(Q[eff.e]) ~~~ "(" ~ m^3/s ~ ")"),
       main = "Bed Load Sites",
       ylim = c(1E-1, 6E3),
       xlim = c(1E-1, 6E3))
  abline(0,1)
  
  plot(Qeff.l, Qeff.e, log = "xy",
       xlab = bquote(italic(Q[eff.l]) ~~~ "(" ~ m^3/s ~ ")"),
       ylab = bquote(italic(Q[eff.e]) ~~~ "(" ~ m^3/s ~ ")"),
       ylim = c(1E-1, 6E3),
       xlim = c(1E-1, 6E3))
  abline(0,1)
dev.off()

######################################################################
pdf("Qeff_Sens_log.pdf", width = 11, height = 6)
  par(mfrow = c(1,3), pty = "s", cex = 1, cex.axis = 0.8)
  
  plot(CV, Qeff.e/mean, #ylim = c(5E-3,500),
       log = "xy",
       ylab = "Qeff.e/mean",
       xlab = "Daily Flow Coefficient of Variation")
  
  points(CV, Qeff.b/mean, pch = 20, col = "gray40")
  points(CV, Qeff.l/mean, pch = "*", col = "blue")
  
  text(0.75, 1E-2, labels = "Qeff.e", cex = 0.8, adj = c(0,0))
  text(1.5, 1E-2, labels = "Qeff.b", cex = 0.8, adj = c(0,0))
  text(3, 1E-2, labels = "Qeff.l", cex = 0.8, adj = c(0,0))
  
  points(0.7, 1.3e-2)
  points(1.4, 1.3e-2, pch = 20, col = "gray40", cex = 1.5)
  points(2.8, 1.3e-2, pch = "*", col = "blue", cex = 1.5)
  
  plot(beta, Qeff.e/mean, log = "xy", xlab = expression(beta), ylab = "Qeff.e/mean", 
       ylim = c(5E-3,500),
       main = "Bed Load Sites")
  points(beta, Qeff.b/mean, pch = 20, col = "gray40")
  points(beta, Qeff.l/mean, pch = "*", col = "blue")
  
  plot(RL, Qeff.e/mean, log = "xy", ylab = "Qeff.e/mean", ylim = c(5E-3,500), 
       xlab = "Flow Record Length [yrs]")
  points(RL, Qeff.b/mean, pch = 20, col = "gray40")
  points(RL, Qeff.l/mean, pch = "*", col = "blue")
dev.off()

######################################################################
pdf("Qeff_Sens.pdf", width = 11, height = 6)
par(mfrow = c(1,3), pty = "s", cex = 1, cex.axis = 0.8)

plot(CV, Qeff.e/mean, ylim = c(0,200),
     ylab = "Qeff.e/mean",
     xlab = "Daily Flow Coefficient of Variation")

points(CV, Qeff.b/mean, pch = 20, col = "gray40")
points(CV, Qeff.l/mean, pch = "*", col = "blue")

plot(beta, Qeff.e/mean, xlab = expression(beta), ylab = "Qeff.e/mean", 
     ylim = c(0,200),
     main = "Bed Load Sites")
points(beta, Qeff.b/mean, pch = 20, col = "gray40")
points(beta, Qeff.l/mean, pch = "*", col = "blue")

text(1.5, 190, labels = "Qeff.e", cex = 0.8, adj = c(0,0))
text(3, 190, labels = "Qeff.b", cex = 0.8, adj = c(0,0))
text(4.5, 190, labels = "Qeff.l", cex = 0.8, adj = c(0,0))

points(1.3, 193)
points(2.8, 193, pch = 20, col = "gray40", cex = 1.5)
points(4.35, 195, pch = "*", col = "blue", cex = 1.5)

plot(RL, Qeff.e/mean, ylab = "Qeff.e/mean", ylim = c(0,200), 
     xlab = "Flow Record Length [yrs]")
points(RL, Qeff.b/mean, pch = 20, col = "gray40")
points(RL, Qeff.l/mean, pch = "*", col = "blue")
dev.off()

# plot(RL, Qeff.e/mean)
# points(RL, Qeff.k/mean, pch = 0)
# points(RL, Qeff.b/mean, pch = 18)
# points(RL, Qeff.m/mean, pch = 20)

############################################################################################
# Are differences between Qeff methods influenced by RL, beta, or CV?
change1 <- Qeff.b/mean - Qeff.e/mean
change2 <- Qeff.b/mean - Qeff.l/mean
change3 <- Qeff.e/mean - Qeff.l/mean

y.min <- floor(min(c(change1, change2, change3)))
y.max <- ceiling(max(c(change1, change2, change3)))

y.min <- -10
y.max <- 10

###################################################################
pdf("Qeff_Compare_Zm.pdf", width = 11, height = 6)
par (mfrow = c(1,3), pty = "s", cex = 1)

plot(RL, change1, ylab = "Qeff1/mean - Qeff2/mean",
     xlab = "Record Length [yrs]",
     main = "Comparing Normalized Qeff by Method \nas a Function of Record Length",
     cex.main = 0.75, cex.axis = 0.9,
     ylim = c(y.min, y.max))
points(RL, change2, pch = 20, col = "gray40")
points(RL, change3, pch = "*", col = "red")

abline(0,0, col = "gray60")
# for text and legend
usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1)) # new relative user coordinates
text(x = 0.6, y = 0.95, labels = "Qeff.b - Qeff.e", adj = c(0,0), cex = .7)
text(x = 0.6, y = 0.9, labels = "Qeff.b - Qeff.l", adj = c(0,0), cex = .7)
text(x = 0.6, y = 0.85, labels = "Qeff.e - Qeff.l", adj = c(0,0), cex = .7)
points(0.55, 0.96)
points(0.55, 0.91, pch = 20, col = "gray40")
points(0.55, 0.86, pch = "*", col = "red")

# legend(x = c(80, 110), y = c(400, 300),
#        c("Qeff.b - Qeff.e", "Qeff.b - Qeff.l", "Qeff.e - Qeff.l"),
#        pch=c(1,18,8), # gives the legend appropriate symbols (lines)
#        col=c("black","black", "red"), 
#        cex = 0.7,
#        bty = "n") # removes legend box

# par(usr = usr0) # reset plot coords to original values

plot(beta, change1, ylab = "Qeff1/mean - Qeff2/mean", xlab = bquote(beta),
     main = "Comparing Normalized Qeff by Method \nas a Function of Beta",
     cex.main = 0.75, cex.axis = 0.9,
     ylim = c(y.min, y.max))
points(beta, change2,pch = 20, col = "gray40")
points(beta, change3, pch = "*", col = "red")
abline(0,0, col = "gray60")
abline(mean(change1),0, lty = 6)
abline(mean(change2),0, lty = 6, col = "gray40")
abline(mean(change3),0, lty = 6, col = "red")

mtext("Bed Load Sites", line = 3, col = "seagreen", font = 4)
# legend("topright",
#          c("Qeff.b - Qeff.e", "Qeff.b - Qeff.l", "Qeff.e - Qeff.l"),
#          pch=c(1,18,8), # gives the legend appropriate symbols (lines)
#          col=c("black","black", "red"), 
#          cex = 0.7,
#          bty = "n")

plot(CV, change1, type = "n", 
     ylab = "Qeff1/mean - Qeff2/mean", 
     xlab = "Daily Flow Coefficient of Variation",
     main = "Comparing Normalized Qeff by Method \nas a Function of CV",
     cex.main = 0.75, cex.axis = 0.9,
     ylim = c(y.min, y.max))
points(CV, change1)
points(CV, change2, pch = 20, col = "gray40")
points(CV, change3, pch = "*", col = "red")
abline(0,0, col = "gray60")

dev.off()

#########################################################################################
pdf("Qeff_Compare.pdf", width = 11, height = 6)
par (mfrow = c(1,3), pty = "s", cex = 1)

plot(RL, change1, ylab = "Qeff1/mean - Qeff2/mean",
     xlab = "Record Length [yrs]",
     main = "Comparing Normalized Qeff by Method \nas a Function of Record Length",
     cex.main = 0.75, cex.axis = 0.9)
points(RL, change2, pch = 20, col = "gray40")
points(RL, change3, pch = "*", col = "red")

abline(0,0, col = "gray60")
# for text and legend
usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1)) # new relative user coordinates
text(x = 0.6, y = 0.95, labels = "Qeff.b - Qeff.e", adj = c(0,0), cex = .7)
text(x = 0.6, y = 0.9, labels = "Qeff.b - Qeff.l", adj = c(0,0), cex = .7)
text(x = 0.6, y = 0.85, labels = "Qeff.e - Qeff.l", adj = c(0,0), cex = .7)
points(0.55, 0.96)
points(0.55, 0.91, pch = 20, col = "gray40")
points(0.55, 0.86, pch = "*", col = "red")

# legend(x = c(80, 110), y = c(400, 300),
#        c("Qeff.b - Qeff.e", "Qeff.b - Qeff.l", "Qeff.e - Qeff.l"),
#        pch=c(1,18,8), # gives the legend appropriate symbols (lines)
#        col=c("black","black", "red"), 
#        cex = 0.7,
#        bty = "n") # removes legend box

# par(usr = usr0) # reset plot coords to original values

plot(beta, change1, ylab = "Qeff1/mean - Qeff2/mean", xlab = bquote(beta),
     main = "Comparing Normalized Qeff by Method \nas a Function of Beta",
     cex.main = 0.75, cex.axis = 0.9)
points(beta, change2,pch = 20, col = "gray40")
points(beta, change3, pch = "*", col = "red")
abline(0,0, col = "gray60")
abline(mean(change1),0, lty = 6)
abline(mean(change2),0, lty = 6, col = "gray40")
abline(mean(change3),0, lty = 6, col = "red")

mtext("Bed Load Sites", line = 3, col = "seagreen", font = 4)
# legend("topright",
#          c("Qeff.b - Qeff.e", "Qeff.b - Qeff.l", "Qeff.e - Qeff.l"),
#          pch=c(1,18,8), # gives the legend appropriate symbols (lines)
#          col=c("black","black", "red"), 
#          cex = 0.7,
#          bty = "n")

plot(CV, change1, type = "n", 
     ylab = "Qeff1/mean - Qeff2/mean", 
     xlab = "Daily Flow Coefficient of Variation",
     main = "Comparing Normalized Qeff by Method \nas a Function of CV",
     cex.main = 0.75, cex.axis = 0.9)
points(CV, change1)
points(CV, change2, pch = 20, col = "gray40")
points(CV, change3, pch = "*", col = "red")
abline(0,0, col = "gray60")

dev.off()

detach(compare)

x <- 1
if (x ==1){
print('same')
} else if (x > 1){
 print('bigger')
} else {
print('smaller')
}

results <- data.frame(site.no = numeric(nrow(flow.sites)), years = numeric(nrow(flow.sites)))

for(i in 1:nrow(flow.sites)){
  path <- paste("DV_Flows/", flow.sites$site.no[i], ".csv", sep = "")
  flow.data <- read.table(path, header = T, stringsAsFactors = F, sep = ",")
  flow0 <- flow.data$val*0.028316847 # Convert to cms
  years <- length(flow0)/365
  results[i,] <- c(flow.sites$site.no[i], years)
}

### Plot different representations of flow regimes ###
pdf("FlowReprSamples.PDF", width = 9, height = 6)
  par(mfrow = c(1,2), pty = "s", oma = c(0,3,0,3))
  
  plot(flow.sort, ecdf.prob, 
       type = "l",
       col = "royalblue", 
       lwd = 1,
       ylab = "Probability",
       xlab = expression("Discharge"~"("*m^3/s*")"),
       cex.lab = 1.2
       )
  lines(xpts1, ypts1, col = rgb(210/255,105/255, 30/255, alpha = 0.7), lwd = 8)
  lines(flow.sort, ecdf.prob, col = "royalblue", lwd = 2)
  lines(x, plnorm(x, u.y, s.y), col = "gray20", lwd = 2, lty = 6)
  text(0.4,1, labels = "CDF", font = 2)
  legend("bottomright",
         legend = c("ECDF", "Log-Smoothed ECDF", "Fittted Lognormal CDF"),
         lwd = c(3,6,3),
         col = c("royalblue", rgb(210/255,105/255, 30/255, alpha = 0.7), "gray20"),
         lty = c(1,1,6),
         bty = "n",
         cex = 0.8)
  
  plot(mids, flow.freq,
       type = "h", 
       lend = 2,
       lwd = 6,
       col = "royal blue",
       ylab = "Density",
       xlab = expression("Discharge"~"("*m^3/s*")"),
       ylim = c(1e-4,10),
       log = "y",
       cex.lab = 1.2
      )
  lines(q.e1, pdf.e1, col = rgb(210/255,105/255, 30/255, alpha = 0.7), lwd = 4)
  #lines(mids, flow.freq, type = "h", lend = 2, lwd = 6, col = "royalblue")
  lines(x,lnorm.PDF, lwd = 3, col = "gray20", lty = 6)
  text(0.7,10, labels = "PDF", font = 2)
  legend("topright",
         legend = c("Histogram", "EPDF", "Fittted Lognormal PDF"),
         lwd = c(6,4,3),
         lty = c(1,1,6),
         col = c("royalblue", rgb(210/255,105/255, 30/255, alpha = 0.7), "gray20"),
         bty = "n",
         cex = 0.8)
dev.off()


