###################### INTER - SITE COMPARISON #############################
# Import all MFA analysis results and create box and scatter plots
###########################################################################

rm(list = ls())
library(vioplot) # Kernel Density Plotting Function
library(MASS)
library(ggplot2)

source("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/minor.ticks.axis.R")
'%!in%' <- function(x,y)!('%in%'(x,y)) # opposite of %in% or match

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/SuspLoad")
susp.reg.results <- read.table("results_susp_reg.csv", stringsAsFactors = F, header = TRUE, sep = ",")
susp.flow <- read.table("results_flow_metrics.csv", stringsAsFactors = F, header = TRUE, sep = ",")
susp.site.info <- read.table("susp_sites_v10.csv", stringsAsFactors = F, header = TRUE, sep = ",")
susp.yield <- read.table("results_yield_metrics.csv", stringsAsFactors = F, header = TRUE, sep = ",")
susp.peaks <- read.table("results_peaks.csv", stringsAsFactors = F, header = TRUE, sep = ",")

# Check to make sure data frame rows lines up
susp.reg.results <- susp.reg.results[order(susp.reg.results$site.no),]
susp.flow <- susp.flow[order(susp.flow$site.no),]
susp.site.info <- susp.site.info[order(susp.site.info$site.no),]
susp.yield <- susp.yield[order(susp.yield$site.no),]
susp.peaks <- susp.peaks[order(susp.peaks$site.no),]

# # remove sites culled after analysis 
# # (D50 too large from bed sed sample OR regulation, river training, etc from aerial photo inspection)
# cull.list <- c(7097000,7164500,11047300,11472900,1481500,11106550,8158000,11465200,9474000,9261000,5454500,4193500,5325000,5420500,5389500,5587455,6088500,6810000,8383000,8447410,5455000,5079000,11481500,8276500,50055000,8354900,8358400,9368000,9379500,5583000,11113000,8066500,9095500,9180500,11473900,5584685,11473900,11482120,11525655,11153900,9251000)
# 
# keep.rows <- which(susp.site.info$site.no %!in% cull.list)
# susp.reg.results <- susp.reg.results[keep.rows,]
# susp.flow <- susp.flow[keep.rows,]
# susp.site.info <- susp.site.info[keep.rows,]
# susp.yield <- susp.yield[keep.rows,]
# susp.peaks <- susp.peaks[keep.rows,]

# Check to see if all sites line up
sum(c(
  sum(susp.reg.results$site.no == susp.flow$site.no) == nrow(susp.site.info),
  sum(susp.reg.results$site.no == susp.peaks$site.no) == nrow(susp.site.info),
  sum(susp.reg.results$site.no == susp.yield$site.no) == nrow(susp.site.info),
  sum(susp.reg.results$site.no == susp.site.info$site.no) == nrow(susp.site.info))) == 4

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload/")
bed.flow <- read.table("results_flow_metrics.csv", stringsAsFactors = F, header = TRUE, sep = ",")
bed.site.info <- read.table("site_summary.csv", stringsAsFactors = F, header = TRUE, sep = ",")
bed.reg.results <- read.table("results_bed_regression.csv", stringsAsFactors = F, header = TRUE, sep = ",")
bed.yield <- read.table("results_yield_metrics.csv", stringsAsFactors = F, header = TRUE, sep = ",")
bed.peaks <- read.table("results_peaks.csv", stringsAsFactors = F, header = TRUE, sep = ",")
# Replace Q1.5 and Q2 with annual max series data where available
# peaks tables are in (cfs)
bed.peaks$Q1.5[which(is.na(bed.peaks$Q1.5))] <- bed.flow$Q1.5[which(is.na(bed.peaks$Q1.5))] / 0.028316847
bed.peaks$Q2[which(is.na(bed.peaks$Q2))] <- bed.flow$Q2[which(is.na(bed.peaks$Q2))] / 0.028316847

# bed.flow <- bed.flow[which(bed.flow$site.no %in% bed.site.info$site.no),]
# bed.flow <- bed.flow[order(bed.flow$site.no),]
# write.csv(bed.flow, "results_flow_metrics.csv", row.names = F)

# Check to see if all sites line up
sum(c(
  sum(bed.reg.results$site.no == bed.flow$site.no) == nrow(bed.site.info),
  sum(bed.reg.results$site.no == bed.peaks$site.no) == nrow(bed.site.info),
  sum(bed.reg.results$site.no == bed.yield$site.no) == nrow(bed.site.info),
  sum(bed.reg.results$site.no == bed.site.info$site.no) == nrow(bed.site.info))) == 4

Q2.s <- susp.peaks$Q2 * 0.028316847
Q2.b <- bed.peaks$Q2 * 0.028316847
Q1.5.s <- susp.peaks$Q1.5 * 0.028316847
Q1.5.b <- bed.peaks$Q1.5 * 0.028316847
Qavg.s <- susp.flow$mean
Qavg.b <- bed.flow$mean

fplus.s <- susp.yield$f.plus
fplus.b <- bed.yield$f.plus

Qeff.s <- susp.yield$Qeff.e
Qeff.b <- bed.yield$Qeff.e
Qeff.crs <- bed.yield$Qeff.e[which(bed.site.info$d50 >= 64)]
Qs50.s <- susp.yield$Q.50
Qs50.b <- bed.yield$Q.50

spread.s <- susp.yield$yield.spread
spread.b <- bed.yield$yield.spread

Qbf.s <- susp.site.info$Qbf.cms
Qbf.b <- bed.site.info$Qbf

Qbf.yield.s <- susp.yield$Qbf.yield
Qbf.yield.b <- bed.yield$Qbf.yield

beta.s <- susp.reg.results$b1_RLM
beta.b <- bed.reg.results$b1_RLM

da.s <- susp.site.info$da.km2
da.b <- bed.site.info$da.km2

cv.s <- susp.flow$cv
cv.b <- bed.flow$cv
cv.crs <- bed.flow$cv[which(bed.site.info$d50 >= 64)]

skew.s <- susp.flow$skewness
skew.b <- bed.flow$skewness
skew.crs <- bed.flow$skewness[which(bed.site.info$d50 >= 64)]

mean.b <- bed.flow$mean
mean.s <- susp.flow$mean

Qeff.mean.s <- Qeff.s / mean.s
Qeff.mean.b <- Qeff.b / mean.b
Q50.mean.s <- Qs50.s / mean.s
Q50.mean.b <- Qs50.b / mean.b

wd.s <- susp.site.info$Wbf.m / susp.site.info$Dbf.m
wd.b <- bed.site.info$Wbf / bed.site.info$Dbf

d50.s <- susp.site.info$d50
d50.b <- bed.site.info$d50

# For Vioplot
# d50.s <- d50.s[which(!is.na(d50.s))]
# d50.b <- d50.b[which(!is.na(d50.b))]
# d50.s <- log(d50.s, base = 2)
# d50.b <- log(d50.b, base = 2)
d84.s <- susp.site.info$d84
d84.b <- bed.site.info$d84

slope.s <- susp.site.info$slope
slope.b <- bed.site.info$slope

tau.star.s <- susp.site.info$Dbf*susp.site.info$slope / (1.65 * susp.site.info$d50/1000)
tau.star.b <- bed.site.info$Dbf*bed.site.info$slope / (1.65 * bed.site.info$d50/1000)
tau.star.s <- tau.star.s[which(!is.na(tau.star.s))]
tau.star.b <- tau.star.b[which(!is.na(tau.star.b))]

RL.s <-susp.flow$yrs
RL.b <-bed.flow$yrs

RIeff.s <- susp.yield$Qeff.RI
RIeff.b <- bed.yield$Qeff.RI
RIeff.crs <- bed.yield$Qeff.RI[which(bed.site.info$d50 >= 64)]
# RI.s <- RI.s[which(!is.na(RI.s))]
# RI.b <- RI.b[which(!is.na(RI.b))]
RIs50.s <- susp.yield$Q50.RI
RIs50.b <- bed.yield$Q50.RI
RIs50.crs <- bed.yield$Q50.RI[which(bed.site.info$d50 >= 64)]

Qeff.ptile.s <- susp.yield$Qeff.ptile*100
Qeff.ptile.b <- bed.yield$Qeff.ptile*100
Q25.ptile.s <- susp.yield$Q25.ptile*100
Q25.ptile.b <- bed.yield$Q25.ptile*100
Q50.ptile.s <- susp.yield$Q50.ptile*100
Q50.ptile.b <- bed.yield$Q50.ptile*100
Q75.ptile.s <- susp.yield$Q75.ptile*100
Q75.ptile.b <- bed.yield$Q75.ptile*100

Cvb.b <- cv.b*beta.b
Cvb.s <- cv.s*beta.s

Cv_b.b <- cv.b/beta.b
Cv_b.s <- cv.s/beta.s

r2.s <- susp.reg.results$r2
r2.b <- bed.reg.results$r2

iatAvg.s <- 1/ susp.flow$lambda50 / 365
iatAvg.b <- 1/ bed.flow$lambda50 / 365

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")

### Box Plots of Qbf.yield ###
#pdf("Qbf_Yield_Percentiles.pdf", width = 6, height = 5)
postscript("Qbf_Yield_Percentiles.eps", width = 5, height = 4, paper = "special")
  par(pty = "s", cex.axis = 1, mar = c(3.1, 3.1, 0, 0), oma = c(2,1,1,1))
  boxplot(Qbf.yield.s, at = 1, col = "seagreen", xlim = c(0.5, 4.5), ylim = c(0,1), 
          ylab = expression("Percent Cumulative Yield"), 
          width = 2,
          cex.lab = 1.2)
  boxplot(1-fplus.s, at = 2, col = "darkseagreen", add = T, width = 2)
  boxplot(Qbf.yield.b, at = 3, col = "royalblue", add = T, width = 2)
  boxplot(1-fplus.b, at = 4, col = "cadetblue3", add = T, width = 2)
 
  axis(side=1, at=c(1,2,3,4),
       labels= c(expression(italic(Q[bf])),expression(italic(Q[eff])),
                 expression(italic(Q[bf])),expression(italic(Q[eff]))))
  axis(2, at = seq(0.1, 0.9, by = 0.1), labels = F, tck = -.015)
  
  mtext(side = 1, at = 1.5, "Fine Bed Sites", line = 3, cex = 1.2)
  mtext(side = 1, at = 3.5, "Coarse Bed Sites", line = 3, cex = 1.2)

dev.off()


### Box Plots of RI's ###############
pdf("Return_Intervals.pdf", width = 6, height = 6)
  par(pty = "m", cex.axis = 1, mar = c(.2, 3, 0, 0), mfrow = c(2,1), oma = c(2,2,1,1))
  layout(matrix(c(1,2), 2, 1, byrow = TRUE), 
       widths=c(1), heights=c(1,4))

  boxplot(RIeff.s, at = 0.5, col = "seagreen3", xlim = c(0, 4.5), ylim = c(5,35),#range(RIeff.s, RIs50.b, RIs50.s, na.rm = T), 
          ylab = "", width = 0.5)
  boxplot(RIeff.b, at = 1.5, col = "royalblue4", add = T, width = 0.5)
  boxplot(RIs50.s, at = 3, col = "seagreen3", add = T, width = 0.5)
  boxplot(RIs50.b, at = 4, col = "royalblue4", add = T, width = 0.5)

  boxplot(RIeff.s, at = 0.5, col = "seagreen3", xlim = c(0, 4.5), ylim = c(1,4),#range(RIeff.s, RIs50.b, RIs50.s, na.rm = T), 
          ylab = "", width = 0.5, yaxt="n", ylab="")
  boxplot(RIeff.b, at = 1.5, col = "royalblue4", add = T, width = 0.5,yaxt="n", ylab="")
  boxplot(RIs50.s, at = 3, col = "seagreen3", add = T, width = 0.5,yaxt="n", ylab="")
  boxplot(RIs50.b, at = 4, col = "royalblue4", add = T, width = 0.5,yaxt="n", ylab="")
  axis(side=1, at=c(0.5, 1.5, 3, 4),
       labels= c("Qeff.s","Qeff.b", "Qs50.s", "Qs50.b"))
  axis(side = 2, at = c(1,2,3,4))
  axis(side = 2, at = c(1.5,2.5,3.5,4.5), lab = rep("",4))
  mtext (side = 2, "Return Interval (yr)", outer = T, line = 0, cex=1.2)
dev.off()

#######################################################

### Box Plots of P-tiles ###############
pdf("Percentiles.pdf", width = 6, height = 5)
par(pty = "m", cex.axis = 0.75, mar = c(5.1, 5.1, 4.1, 2.1))
boxplot(Qeff.ptile.s, at = 0.5, col = "seagreen3", xlim = c(0, 9.5), ylim = c(0,100), 
        ylab = "Percentile of Daily Flow", width = 0.5)
boxplot(Qeff.ptile.b, at = 1.5, col = "royalblue4", add = T, width = 0.5)
boxplot(Q25.ptile.s, at = 3, col = "seagreen3", add = T, width = 0.5)
boxplot(Q25.ptile.b, at = 4, col = "royalblue4", add = T, width = 0.5)
boxplot(Q50.ptile.s, at = 5.5, col = "seagreen3", add = T, width = 0.5)
boxplot(Q50.ptile.b, at = 6.5, col = "royalblue4", add = T, width = 0.5)
boxplot(Q75.ptile.s, at = 8, col = "seagreen3", add = T, width = 0.5)
boxplot(Q75.ptile.b, at = 9, col = "royalblue4", add = T, width = 0.5)
axis(side=1, at=c(0.5, 1.5, 3, 4, 5.5, 6.5, 8, 9),
     labels= c("Qeff.s","Qeff.b", "Q25.s", "Q25.b", "Q50.s", "Q50.b", "Q75.s", "Q75.b"))
axis(2, at = seq(0.1, 0.9, by = 0.1), labels = F, tck = -.02)

legend("bottomright", c("Fine Bed", "Coarse Bed"), 
       fill = c("seagreen3", "royalblue4"), 
       border = "black",
       bty = "n")
dev.off()

####### Site Specific Metrics ################

pdf( "Site_Compare_VioPlots2.pdf", width = 4, height = 4.5)
  par(pty = "s", mfrow = c(3,2), 
      mar = c(0,5.1,0,0), oma = c(5.1,0.5,1,1), xpd = T)
  cex.lab = 0.9
  cex.ax = 0.8
  x.txt <- 0.1
  y.txt <- 0.9
  ##### Beta #####
  plot(1, 1, xlim=c(1.5,3.5), ylim = range(c(beta.s, beta.b)), type="n",
       xlab = NA, 
       ylab = bquote(italic(beta)), 
       cex.lab = cex.lab,
       cex.axis = cex.ax,
       axes = F) #main = "Sediment Rating Curve Exponent"
  vioplot(beta.s, at = 2, col = "forestgreen", add=T)
  vioplot(beta.b, at = 3, col = "royalblue", add = T)
  
  #axis(side=1,at=2:3,labels=c("Suspended Load","Bed Load"))
  axis(side=2, cex.axis = cex.ax)
  text(x=grconvertX(x.txt, from='npc'), y = grconvertY(y.txt, from='npc'), labels = "(a)", cex = 0.9)
  
  ##### DA #####
  plot(1, 1, xlim=c(1.5,3.5), 
       ylim = c(floor(min(log10(range(c(da.s, da.b))))), ceiling(max(log10(range(c(da.s, da.b))))) ),
       type="n",
       xlab = "", 
       ylab = bquote("Drainage Area  ("~km^2~")"),
       axes = F,
       cex.lab = cex.lab,
       cex.axis = cex.ax) #main = "Drainage Area"
  vioplot(log10(da.s), at = 2, col = "forestgreen", add=T)
  vioplot(log10(da.b), at = 3, col = "royalblue", add = T)
  
  at1 <- floor(min(log10(range(c(da.s, da.b))))):ceiling(max(log10(range(c(da.s, da.b)))))

  labels <- sapply(at1, function(i) as.expression(bquote(10^ .(i))) )
  at1[which(at1 %% 2 == 0)] <- NA

  at2 <- floor(min(log10(range(c(da.s, da.b))))):ceiling(max(log10(range(c(da.s, da.b)))))
  at2 <- at2[which(at2 %% 2 == 0)]
  #axis(side=1,at=2:3,labels=c("Suspended Load","Bed Load"))
  axis(side = 2, at = at1, labels = labels, cex.axis = cex.ax)
  axis(side = 2, labels = F, at = at2, tck = -0.03, cex.axis = cex.ax)
  text(x=grconvertX(x.txt, from='npc'), y = grconvertY(y.txt, from='npc'), labels = "(b)", cex = 0.9)

  
  ##### CV #####
  plot(1, 1, xlim=c(1.5,3.5), ylim = range(c(cv.s, cv.b)), type="n",
       xlab = "", ylab = expression(italic(C[v])), axes = F,
       cex.lab = cex.lab,
       cex.axis = cex.ax) #main = "Flow Coefficient of Variation"
  vioplot(cv.s, at = 2, col = "forestgreen", add=T)
  vioplot(cv.b, at = 3, col = "royalblue", add = T)
  
  #axis(side=1,at=2:3,labels=c("Suspended Load","Bed Load"))
  axis(side=2, cex.axis = cex.ax)
  text(x=grconvertX(x.txt, from='npc'), y = grconvertY(y.txt, from='npc'), labels = "(c)", cex = 0.9)

  
  #### RL #####
  plot(1, 1, xlim=c(1.5,3.5), ylim = range(c(RL.s, RL.b)), type="n",
       xlab = "", ylab = "Flow Record Length (yrs)", axes = F,
       cex.lab = cex.lab,
       cex.axis = cex.ax)
  vioplot(RL.s, at = 2, col = "forestgreen", add=T)
  vioplot(RL.b, at = 3, col = "royalblue", add = T)
  
  #axis(side=1,at=2:3,labels=c("Suspended Load","Bed Load"))
  axis(side = 2, cex.axis = cex.ax, at = c(0,20,40,60,80,100),
       labels = c(NA, 20, 40, 60, 80, 100))
  axis(side = 2, labels = F, at = c(10, 30, 50, 70, 90, 110), tck = -0.03, 
       cex.axis = cex.ax)
  text(x=grconvertX(x.txt, from='npc'), y = grconvertY(y.txt, from='npc'), labels = "(d)", cex = 0.9)

  ##### d50 #####
  plot(1, 1, xlim=c(1.5,3.5), ylim = range(c(d50.s, d50.b), na.rm = T), type="n",
       xlab = "", 
       ylab =  expression(italic(D)[50]~~"(mm)"), 
       axes = F,
       cex.lab = cex.lab,
       cex.axis = cex.ax)
  vioplot(d50.s, at = 2, col = "forestgreen", add=T, na.rm = T)
  vioplot(d50.b, at = 3, col = "royalblue", add = T, na.rm = T)
  
  at <- floor(min(c(d50.s, d50.b))):ceiling(max(c(d50.s, d50.b)))
  labels <- 2^at
  
  axis(side=1,at=2:3, cex.axis = cex.ax, labels = c("Suspended \nLoad", "Bed Load"), tck = -0.04)
  axis(side=2, at = at, labels = labels, cex.axis = cex.ax)
  text(x=grconvertX(x.txt, from='npc'), y = grconvertY(y.txt, from='npc'), labels = "(e)", cex = 0.9)

  
  ##### tau.star #####
  ymin <- -2 #floor(min(log10(range(c(tau.star.s, tau.star.b)))))
  ymax <- ceiling(max(log10(range(c(tau.star.s, tau.star.b)))))
  plot(1, 1, xlim=c(1.5,3.5), 
       ylim = c(ymin, ymax),
       type="n",
       xlab = "", 
       ylab = expression(italic(tau*"*")~~"@ Bankfull"), 
       axes = F, 
       cex.lab = cex.lab,
       cex.axis = cex.ax)
  vioplot(log10(tau.star.s), at = 2, col = "forestgreen", add=T)
  vioplot(log10(tau.star.b), at = 3, col = "royalblue", add = T)
  
  at <- floor(min(log10(range(c(tau.star.s, tau.star.b))))):
        ceiling(max(log10(range(c(tau.star.s, tau.star.b)))))
  
  labels <- sapply(at, function(i)
    as.expression(bquote(10^ .(i))) )
  
  axis(side=1,at=2:3,cex.axis = cex.ax, labels = c("Suspended \nLoad", "Bed Load"), tck = -0.04)
  axis(side=2, at = at, labels = labels, cex.axis = cex.ax)
  minor.ticks.axis(2,9, mn = ymin, mx = ymax, cex = cex.ax)
  text(x=grconvertX(x.txt, from='npc'), y = grconvertY(y.txt, from='npc'), labels = "(f)", cex = 0.9)


dev.off()



#############################Inter-Arrival Time v Skew and CV #################

par(mfrow = c(1,2), mar = c(3,0,0,0), oma = c(0,5,1,1), pty = "s")
plot(cv.s,iatAvg.s, ylim = c(0,20),
     xlab = expression(italic(C[v])))
#points(cv.b, iatAvg.b, pch = "*", cex = 1.2)
lines(loess.smooth(cv.s, iatAvg.s), 
      lwd = 2, 
      col = rgb( 46/255,139/255,87/255, alpha = 0.75),
      lty = 1)

#plot(cv.b, bed.flow$lambda50, ylim = c(0,0.02))

plot(skew.s, iatAvg.s,
     xlab = "Skewness",
     ylab = "",
     yaxt = "n",
     ylim = c(0,20))
lines(loess.smooth(skew.s, iatAvg.s), 
      lwd = 2, 
      col = rgb( 46/255,139/255,87/255, alpha = 0.75),
      lty = 1)
#points(skew.b, iatAvg.b, pch = "*", cex = 1.2)
mtext("Interarrival Time (yr)", side = 2, outer = T, line = 3)

###################################################################
plot(cv.b,iatAvg.b, #ylim = c(0,20),
     xlab = expression(italic(C[v])),
     pch = "*", cex = 1.2)
#points(cv.b, iatAvg.b, pch = "*", cex = 1.2)
lines(loess.smooth(cv.b, iatAvg.b), 
      lwd = 2, 
      col = rgb( 65/255,105/255,225/255, alpha = 0.75),
      lty = 1)

#plot(cv.b, bed.flow$lambda50, ylim = c(0,0.02))

plot(skew.b, iatAvg.b,
     xlab = "Skewness",
     ylab = "",
     yaxt = "n",
     #ylim = c(0,20),
     pch = "*",
     cex = 1.2
     )
lines(loess.smooth(skew.b, iatAvg.b), 
      lwd = 2, 
      col = rgb( 65/255,105/255,225/255, alpha = 0.75),
      lty = 1)
#points(skew.b, iatAvg.b, pch = "*", cex = 1.2)
mtext("Interarrival Time (yr)", side = 2, outer = T, line = 3)


################### CV Relations ###############

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")
pdf( "Cv_fplus_spread.pdf", width = 4, height = 6)

par(mfrow = c(2,1), pty = "s", mar = c(0, 0, 0, 0), oma = c(4,4,1,1))
#   par(mfrow = c(2,2), pty = "s", 
#       mar = c(0,5,0,0), 
#       oma = c(5,1,1,1), 
#       cex.axis = 0.9, 
#       cex.lab = 1, 
#       mgp = c(2, 0.6, 0))
cex.axis = 1

#   ### Qs50/Q1.5 = f(Cv)
#   plot(cv.s, log10(Qs50.s/mean.s), 
#        xlim = range(cv.s),
#        ylim = c(-0.5, 2.5),
#        axes = F,
#        col = "gray40",
#        ylab = "",
#        xlab = "")
#   points(cv.b, log10(Qs50.b/mean.b), pch = "*", cex = 1.5)
#   minor.ticks.axis(2,9, mn = -1, mx = 3, cex = cex.axis)
# 
#   box()
#   mtext(expression(italic(Q)[s50]/italic(bar(X)[flow])), side = 2, line = 2.2)

## f+ vs. Cv
plot(cv.s, fplus.s,
     xaxt = "n",
     yaxt = "n",
     col = "gray40",
     ylim = c(0,1),
     cex.axis = cex.axis,
     ylab = "",
     xlab = "")
lines(loess.smooth(cv.s, fplus.s), lwd = 2, col = "seagreen", lty = 6)
curve(1 - pnorm( - sqrt( log(1 + x^2) ) ), add = T, lwd = 2)
points(cv.b, fplus.b,  pch = "*", cex = 1.5)
axis(side = 2, at = c(0, 0.25, 0.50, 0.75, 1), labels = c(0, 0.25, 0.50, 0.75, 1), cex = cex.axis)
mtext(expression(italic(f)*"+"), side = 2, line = 2.2, cex = 1.2)
legend("bottomright",
       legend = c("Lognormal \nRelation", "LOESS \nFine"),
       lty = c(1,6),
       lwd = c(1.5,1.5),
       cex = 0.6,
       col = c("black", "seagreen"),
       bty = "n",
       y.intersp = 1.5,
       seg.len = 3)
box()

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(b)", cex = 1.25)
par(usr = usr0)

#   ## RI.Qs50
#   plot(cv.s, Q50.ptile.s,
#        xaxt = "n",
#        col = "gray40",
#        cex.axis = cex.axis,
#        ylab = "",
#        xlab = "")
#   points(cv.b, Q50.ptile.b,  pch = "*", cex = 1.5)
#   axis(side = 1, cex.axis = cex.axis)
#   box()
#   mtext(expression("Percentile" ~~ italic(Q)[s50]), side = 2, line = 2.2)
#   mtext(expression(italic(C[v])), side = 1, line = 2.2)

## yield Spread
plot(cv.s, spread.s,
     col = "gray40",
     cex.axis = cex.axis,
     yaxt = "n",
     ylim = c(0,4.1),
     ylab = "",
     xlab = "")
points(cv.b, spread.b, pch = "*", cex = 1.5)
axis(side = 2, at = c(0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4), 
     labels = c(NA, 1, NA, 2, NA, 3, NA, 4), cex.axis = cex.axis)
mtext("Yield Spread", side = 2, line = 2.2, cex = 1.2)
axis(side = 1, cex.axis = cex.axis)
mtext(expression(italic(C[v])), side = 1, line = 2.2, cex = 1.2)

legend("topright",
       legend = c("Fine Bed", "Coarse Bed"),
       pch = c(1,8),
       cex = 0.6,
       y.intersp = 0.8,
       yjust = 1,
       bty = "n")

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(d)", cex = 1.25)
par(usr = usr0)

dev.off()


############  Cv - Qeff - Qs50 Relations ##################
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")

pdf("Cv_Qeff_Qs50.pdf", height = 6, width = 4)

par(mfrow = c(2,1), pty = "s", mar = c(0, 0, 0, 0), oma = c(4,4,1,1))

# Qeff / mean vs. Cv
plot(cv.b, log10(Qeff.b/mean.b), pch = "*", cex.lab = 1.25, cex.axis = 1,
     cex = 1.5,
     ylim = c(-1.5, 2.5),#log10(range(c(Qeff.b/mean.b, Qeff.s/mean.s))),
     xlim = range(c(cv.b, cv.s)),
     axes = F,
     ylab = "",
     xlab = expression(italic(C[v]))
)
points(cv.s, log10(Qeff.s/mean.s))

lines(loess.smooth(cv.b, log10(Qeff.b/mean.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(cv.s, log10(Qeff.s/mean.s)),
      lwd = 2, col = "seagreen", lty = 6)

minor.ticks.axis(2,9, mn = -1, mx = 6, cex = 1)
box()
mtext(expression(italic(Q[eff]/bar(X)[flow])), side = 2, cex = 1.25, line = 3)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(a)", cex = 1.25)
par(usr = usr0)

# Qs50 / mean vs. Cv
plot(cv.b, log10(Qs50.b/mean.b), pch = "*", cex.lab = 1.25,
     cex = 1.5,
     ylim = log10(range(c(Qs50.b/mean.b, Qs50.s/mean.s))),
     xlim = range(c(cv.b, cv.s)),
     axes = F,
     ylab = "", #expression(italic(Q[eff]/bar(X)[flow])),
     xlab = ""
)
points(cv.s, log10(Qs50.s/mean.s))

lines(loess.smooth(cv.b, log10(Qs50.b/mean.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(cv.s, log10(Qs50.s/mean.s)),
      lwd = 2, col = "seagreen", lty = 6)

minor.ticks.axis(2,9, mn = -1, mx = 6, cex = 1)
axis(side = 1)
box()
mtext(expression(italic(C[v])), side = 1, line = 3, cex = 1.25)
mtext(expression(italic(Q[s50]/bar(X)[flow])), side = 2, cex = 1.25, line = 3)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(c)", cex = 1.25)
par(usr = usr0)

legend("bottomright",
       c("LOESS Fine", "LOESS Coarse"),
       lty = c(6, 2),
       col = c("seagreen", "royalblue"),
       bty = "n",
       cex = 0.6,
       lwd = c(2,2),
       seg.len = 4)

dev.off()

###### Beta Relations #####################
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")
pdf( "Beta_Relations.pdf", width = 6, height = 5)
par(mfrow = c(2,2), pty = "s", 
    mar = c(0,5,0,0), 
    oma = c(5,1,1,1), 
    cex.axis = 0.9, 
    cex.lab = 1, 
    mgp = c(2, 0.6, 0))
cex.axis = 1

### Qs50/Q1.5.s = f(beta)
plot(beta.s, log10(Qs50.s/Q1.5.s), 
     xlim = range(beta.b),
     axes = F,
     ylim = c(-1.9, 1.9),
     ylab = "",
     col = "gray40")
points(beta.b, log10(Qs50.b/Q1.5.b), pch = "*", cex = 1.5)
abline(0,0)
minor.ticks.axis(2,9, mn = -2, mx = 2, cex = cex.axis)
#   at <- -2:2
#   labels <- sapply(at, function(i) as.expression(bquote(10^ .(i))) )
#   #   labels[1] <- NA
#   #   labels[length(at)] <- NA
#   axis(side = 2, at = at, labels = labels, cex.axis = cex.axis)
#     at <- floor(min(log10(range(c(Qs50.s/Q1.5.s, Qs50.b/Q1.5.b), na.rm = T)))):
#           ceiling(max(log10(range(c(Qs50.s/Q1.5.s, Qs50.b/Q1.5.b), na.rm = T))))
#     at <- -2:2
#     labels <- sapply(at, function(i) as.expression(bquote(10^ .(i))) )
#     at[which(at %% 2 != 0)] <- NA
#       axis(side = 2, at = at, labels = labels, cex.axis = cex.axis)
#     at <- -2:2    
#     at[which(at %% 2 == 0)] <- NA
#     axis(side = 2, at = at, labels = rep("", length(at)) )

lines(loess.smooth(beta.b, log10(Qs50.b/Q1.5.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(beta.s, log10(Qs50.s/Q1.5.s)),
      lwd = 2, col = "seagreen", lty = 6)
box()
mtext(expression(italic(Q)[s50]/italic(Q)[1.5]), side = 2, line = 2.2, cex = 1.25)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.9, .9, labels = "(a)", cex = 1.25)
par(usr = usr0)

### fplus = f(beta)
plot(beta.s, fplus.s, 
     xlim = range(beta.b),
     ylim = c(0,1),
     xaxt = "n",
     ylab = "",
     col = "gray40")
points(beta.b, fplus.b, pch = "*", cex = 1.5)
mtext(expression(italic(f)*"+"), side = 2, line = 2.2, cex = 1.25)
box()

lines(loess.smooth(beta.b, fplus.b),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(beta.s, fplus.s),
      lwd = 2, col = "seagreen", lty = 6)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.9, .9, labels = "(b)", cex = 1.25)
par(usr = usr0)

### RI = f(beta)
plot(beta.s, log10(RIs50.s), col = "gray40",
     ylim = c(-0.5, 1.5),
     xlim = range(beta.b),
     axes = F,
     ylab = ""
)
points(beta.b, log10(RIs50.b), pch = "*", cex = 1.5 )
box()
axis(side = 1, cex.axis = cex.axis)
minor.ticks.axis(2,9, mn = -0.5, mx = 2, cex = cex.axis)
mtext(expression(italic(RI) ~~ italic(Q)[s50] ~~~ "(yrs)"), side = 2, line = 2.2, cex = 1.25)
mtext(expression(beta), side = 1, line = 2.2, cex = 1.25)

lines(loess.smooth(beta.b, log10(RIs50.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(beta.s, log10(RIs50.s)),
      lwd = 2, col = "seagreen", lty = 6)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.9, .9, labels = "(c)", cex = 1.25)
par(usr = usr0)

legend("bottomright",
       c("Fine Bed", "Coarse Bed", "LOESS Fine", "LOESS Coarse"),
       pch = c(1,8, NA, NA),
       pt.cex = 1,
       lty = c(NA, NA, 6, 2),
       col = c("gray40", "black", "seagreen", "royalblue"),
       bty = "n",
       cex = 0.6,
       y.intersp = 0.9,
       lwd = c(1,1,2,2),
       seg.len = 4)

### spread = f(beta)
plot(beta.s, spread.s, 
     xlim = range(beta.b),
     axes = F,
     ylim = c(0,4.1),
     col = "gray40",
     ylab = "")
points(beta.b, spread.b, pch = "*", cex = 1.5)
axis(side = 1, cex.axis = cex.axis)
box()
mtext(expression(beta), side = 1, line = 2.2, cex = 1.25)
axis(side = 2, at = c(0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4), 
     labels = c(NA, 1, NA, 2, NA, 3, NA, 4), cex.axis = cex.axis)
mtext("Yield Spread", side = 2, line = 2.2, cex = 1.25)

lines(loess.smooth(beta.b, spread.b),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(beta.s, spread.s),
      lwd = 2, col = "seagreen", lty = 6)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.9, .9, labels = "(d)", cex = 1.25)
par(usr = usr0)

dev.off()

################ Cv / Beta Vs. Qeff and Qs50 #####################
pdf("Cv_B_plots3.pdf", height = 6, width = 4)
par(mfrow = c(2,1), pty = "s", mar = c(0, 0, 0, 0), oma = c(4,4,1,1))

# Cv_B vs. Qeff / Q1.5
plot(log10(Cv_b.b), log10(Qeff.b/Q1.5.b), pch = "*", cex = 1.5,
     ylim = log10(range(c(Qeff.b/mean.b, Qeff.s/mean.s))),#range(c(Qeff.b/Q1.5.b, Qeff.s/Q1.5.s), na.rm = T),
     xlim = log10(c(0.1, 10)),#log10(range(c(Cv_b.b, Cv_b.s))),
     axes = F,
     ylab = "", 
     xlab = "")
points(log10(Cv_b.s), log10(Qeff.s/Q1.5.s))

lines(loess.smooth(log10(Cv_b.b), log10(Qeff.b/Q1.5.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(log10(Cv_b.s), log10(Qeff.s/Q1.5.s)),
      lwd = 2, col = "seagreen", lty = 6)

minor.ticks.axis(2,9, mn = -3, mx = 6, cex = 1)
box()
mtext(expression(italic(Q[eff]/Q)[1.5]), side = 2, cex = 1.25, line = 3)

legend("topright",
       c("Fine Bed", "Coarse Bed", "LOESS Fine", "LOESS Coarse"),
       pch = c(1,8, NA, NA),
       pt.cex = 1.2,
       lty = c(NA, NA, 6, 2),
       col = c("black", "black", "seagreen", "royalblue"),
       bty = "n",
       cex = 0.8,
       lwd = c(1,1,2,2),
       seg.len = 4)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(b)", cex = 1.25)
par(usr = usr0)

# Cv_B vs. Qs50 / Q1.5
plot(log10(Cv_b.b), log10(Qs50.b/Q1.5.b), pch = "*", cex = 1.5,
     ylim = log10(range(c(Qs50.b/Q1.5.b, Qs50.s/Q1.5.s), na.rm = T)),
     xlim = log10(c(0.1, 10)),#log10(range(c(Cv_b.b, Cv_b.s))),
     axes = F,
     ylab = "", 
     xlab = "")
points(log10(Cv_b.s), log10(Qs50.s/Q1.5.s))

lines(loess.smooth(log10(Cv_b.b), log10(Qs50.b/Q1.5.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(log10(Cv_b.s), log10(Qs50.s/Q1.5.s)),
      lwd = 2, col = "seagreen", lty = 6)

minor.ticks.axis(2,9, mn = -3, mx = 6, cex = 1)
minor.ticks.axis(1,9, mn = -3, mx = 2, cex = 1)
mtext(expression(italic(C[v]/beta)), side = 1, line = 3, cex = 1.25)
mtext(expression(italic(Q[s50]/Q)[1.5]), side = 2, cex = 1.25, line = 3)
box()

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(d)", cex = 1.25)
par(usr = usr0)

dev.off()


### Plot Qbf/Qeff vs Cv  #################
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")
pdf( "Qbf_Qeff_v_Cv.pdf", width = 4, height = 6)

par(pty = "s", mfrow = c(2,1), mar = c(0, 4, 0, 0), oma = c(4, 1, 1, 1) )
plot(log10(cv.s), log10(Qeff.s / Qbf.s), 
     ylim = c(-2, 1.5),
     xlim = c(-0.5, 1),
     axes = F,
     ylab = expression( italic(Q[eff])/italic(Q[BF]) ), 
     col = "gray40",
     cex.lab = 1.25
)
abline(0,0, col = "gray40")
points(log10(cv.b), log10(Qeff.b / Qbf.b), pch = "*", cex = 1.5)
lines(loess.smooth(log10(cv.b), log10(Qeff.b / Qbf.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(log10(cv.s), log10(Qeff.s / Qbf.s)),
      lwd = 2, col = "seagreen", lty = 6)

minor.ticks.axis(2,9, mn = -2, mx = 6, cex = 1)
box()

legend("topleft",
       c("Fine Bed", "Coarse Bed", "LOESS Fine", "LOESS Coarse"),
       pch = c(1,8, NA, NA),
       pt.cex = 1,
       lty = c(NA, NA, 6, 2),
       col = c("gray40", "black", "seagreen", "royalblue"),
       bty = "n",
       cex = 0.6,
       y.intersp = 0.9,
       lwd = c(1,1,2,2),
       seg.len = 2.5)


### Plot Qbf/Qeff vs Cv #################
plot(log10(cv.s), log10(Qs50.s / Qbf.s ), 
     ylim = c(-1.2, 1.2),
     xlim = c(-0.5, 1),
     ylab = expression( italic(Q[s50]) / italic(Q[BF])),
     axes = F,
     col = "gray40",
     cex.lab = 1.25)
abline(0,0, col = "gray40")
points(log10(cv.b), log10(Qs50.b / Qbf.b), pch = "*", cex = 1.5)

lines(loess.smooth(log10(cv.b), log10(Qs50.b / Qbf.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(log10(cv.s), log10(Qs50.s / Qbf.s)),
      lwd = 2, col = "seagreen", lty = 6)


minor.ticks.axis(2,9, mn = -1, mx = 6, cex = 1)
minor.ticks.axis(1,9, mn = -1, mx = 2, cex = 1)
box()

mtext(side = 1, expression( italic(C[v]) ), line = 2, cex = 1.25)

dev.off()

################# beta = f(D50) ####################################
phi.s <- log(susp.site.info$d50, base = 2)
phi.b <- log(bed.site.info$d50, base = 2)

setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")
pdf( "Beta_Phi.pdf", width = 4, height = 4)
par(pty = "s", mar = c(4, 5.2, 1, 1), mfrow = c(1,1), oma = c(0, 0, 0, 0))
plot(phi.b, beta.b, pch = "*", cex = 1.5,
     xlim = range(c(phi.s, phi.b), na.rm = T),
     xlab = expression(italic(D)[50] ~~~ "(mm)"),
     ylab = expression(beta),
     cex.axis = 1,
     xaxt = "n",
     cex.lab = 1.5)
axis(side = 1, at = seq(-2, 8, by = 2),
     labels = 2^seq(-2, 8, by = 2), 
     cex.axis = 1)
points(phi.s, beta.s, col = "gray40", cex = 1.1)


legend("topleft",
       legend = c("Fine Bed", "Coarse Bed"),
       pch = c(1,8),
       cex = 0.9,
       y.intersp = 0.8,
       yjust = 1,
       bty = "n")

dev.off()

######################### PLOT Beta relationships ###############
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")
pdf( "Beta_Relations.pdf", width = 3, height = 8)
par(pty = "s", 
    mar = c(0, 3, 0, 0), 
    mfrow = c(4,1), 
    oma = c(4.5, 0.5, 0.5, 0.5))
### Qs50/Q1.5.s = f(beta)
plot(beta.s, Qs50.s/Q1.5.s, log = "y", 
     xlim = range(beta.b),
     xaxt = "n",
     ylab = expression(italic(Q)[s50]/italic(Q)[1.5]),
     col = "gray40")
points(beta.b, Qs50.b/Q1.5.b, pch = "*", cex = 1.5)
abline(0,0)

### fplus = f(beta)
plot(beta.s, fplus.s, 
     xlim = range(beta.b),
     xaxt = "n",
     ylab = expression(italic(f)*"+"),
     col = "gray40")
points(beta.b, fplus.b, pch = "*", cex = 1.5)

### spread = f(beta)
plot(beta.s, spread.s, 
     xlim = range(beta.b),
     xaxt = "n",
     ylab = "Yield Spread",
     col = "gray40")
points(beta.b, spread.b, pch = "*", cex = 1.5)

### RI = f(beta)
plot(beta.s, susp.yield$Qeff.RI, col = "gray40",
     log = "y", 
     xlim = range(beta.b),
     ylab = expression(italic(RI) ~~ "of" ~ italic(Q[eff]))
)
points(beta.b, bed.yield$Qeff.RI, pch = "*", cex = 1.5 )
mtext(expression(beta), side = 1, line = 3)
legend(x = 4.5, y = 1e-2,
       legend = c("Susp. Load", "Bed Load"),
       pch = c(1,8),
       cex = 0.8,
       y.intersp = 0.8,
       yjust = 1,
       bty = "n")
dev.off()

#############  Cv*Beta Relations ##################
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")

### Cv*Beta Vs. Qeff and Qs50 ###
pdf("CvB_plots_LOESS.pdf", height = 6, width = 4)

par(mfrow = c(2,1), pty = "s", mar = c(0, 0, 0, 0), oma = c(4,4,1,1))

# Qeff / mean vs. Cv*B
plot(log10(Cvb.b), log10(Qeff.b/mean.b), pch = "*", cex.lab = 1.25, cex.axis = 1,
     cex = 1.5,
     ylim = c(-1.5, 2.5),#log10(range(c(Qeff.b/mean.b, Qeff.s/mean.s))),
     xlim = log10(range(c(Cvb.b, Cvb.s))),
     axes = F,
     ylab = "",
     xlab = expression(italic(C[v]*beta))
)
points(log10(Cvb.s), log10(Qeff.s/mean.s))

lines(loess.smooth(log10(Cvb.b), log10(Qeff.b/mean.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(log10(Cvb.s), log10(Qeff.s/mean.s)),
      lwd = 2, col = "seagreen", lty = 6)

minor.ticks.axis(2,9, mn = -1, mx = 6, cex = 1)
box()
mtext(expression(italic(Q[eff]/bar(X)[flow])), side = 2, cex = 1.25, line = 3)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(a)", cex = 1.25)
par(usr = usr0)

# Qs50 / mean vs. Cv*B
plot(log10(Cvb.b), log10(Qs50.b/mean.b), pch = "*", cex.lab = 1.25,
     cex = 1.5,
     ylim = log10(range(c(Qs50.b/mean.b, Qs50.s/mean.s))),
     xlim = log10(range(c(Cvb.b, Cvb.s))),
     axes = F,
     ylab = "", #expression(italic(Q[eff]/bar(X)[flow])),
     xlab = ""
)
points(log10(Cvb.s), log10(Qs50.s/mean.s))

lines(loess.smooth(log10(Cvb.b), log10(Qs50.b/mean.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(log10(Cvb.s), log10(Qs50.s/mean.s)),
      lwd = 2, col = "seagreen", lty = 6)

minor.ticks.axis(2,9, mn = -1, mx = 6, cex = 1)
minor.ticks.axis(1,9, mn = -1, mx = 2, cex = 1)
box()
mtext(expression(italic(C[v]*beta)), side = 1, line = 3, cex = 1.25)
mtext(expression(italic(Q[s50]/bar(X)[flow])), side = 2, cex = 1.25, line = 3)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(d)", cex = 1.25)
par(usr = usr0)

dev.off()

### Cv / Beta Vs. Qeff and Qs50 ###
pdf("Cv_B_plots.pdf", height = 6, width = 4)
par(mfrow = c(2,1), pty = "s", mar = c(0, 0, 0, 0), oma = c(4,4,1,1))

# Cv_B vs. Qeff / Q1.5
plot(log10(Cv_b.b), log10(Qeff.b/Q1.5.b), pch = "*", cex = 1.5,
     ylim = log10(range(c(Qeff.b/mean.b, Qeff.s/mean.s))),#range(c(Qeff.b/Q1.5.b, Qeff.s/Q1.5.s), na.rm = T),
     xlim = log10(c(0.1, 10)),#log10(range(c(Cv_b.b, Cv_b.s))),
     axes = F,
     ylab = "", 
     xlab = "")
points(log10(Cv_b.s), log10(Qeff.s/Q1.5.s))

lines(loess.smooth(log10(Cv_b.b), log10(Qeff.b/Q1.5.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(log10(Cv_b.s), log10(Qeff.s/Q1.5.s)),
      lwd = 2, col = "seagreen", lty = 6)

minor.ticks.axis(2,9, mn = -3, mx = 6, cex = 1)
box()
mtext(expression(italic(Q[eff]/Q)[1.5]), side = 2, cex = 1.25, line = 3)

legend("topright",
       c("Susp. Sites", "Bed Sites", "LOESS Susp.", "LOESS Bed"),
       pch = c(1,8, NA, NA),
       pt.cex = 1.2,
       lty = c(NA, NA, 6, 2),
       col = c("black", "black", "seagreen", "royalblue"),
       bty = "n",
       cex = 0.8,
       lwd = c(1,1,2,2),
       seg.len = 4)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(b)", cex = 1.25)
par(usr = usr0)

# Cv_B vs. Qs50 / Q1.5
plot(log10(Cv_b.b), log10(Qs50.b/Q1.5.b), pch = "*", cex = 1.5,
     ylim = log10(range(c(Qs50.b/Q1.5.b, Qs50.s/Q1.5.s), na.rm = T)),
     xlim = log10(c(0.1, 10)),#log10(range(c(Cv_b.b, Cv_b.s))),
     axes = F,
     ylab = "", 
     xlab = "")
points(log10(Cv_b.s), log10(Qs50.s/Q1.5.s))

lines(loess.smooth(log10(Cv_b.b), log10(Qs50.b/Q1.5.b)),
      lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(log10(Cv_b.s), log10(Qs50.s/Q1.5.s)),
      lwd = 2, col = "seagreen", lty = 6)

minor.ticks.axis(2,9, mn = -3, mx = 6, cex = 1)
minor.ticks.axis(1,9, mn = -3, mx = 2, cex = 1)
mtext(expression(italic(C[v]/beta)), side = 1, line = 3, cex = 1.25)
mtext(expression(italic(Q[s50]/Q)[1.5]), side = 2, cex = 1.25, line = 3)
box()

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(e)", cex = 1.25)
par(usr = usr0)

dev.off()

### Cv_B vs. fplus and Yield Spread
pdf("Cv_B_plots2.pdf", height = 6, width = 4)

par(mfrow = c(2,1), pty = "s", mar = c(0, 0, 0, 0), oma = c(4,4,1,1))

# Cv_Beta vs. f.plus 
plot(Cv_b.b, fplus.b, pch = "*", cex = 1.5, 
     xlim = range(c(Cv_b.b, Cv_b.s)),
     xaxt = "n",
     xlab = "",
     ylab = "",
     cex.lab = 1.25,
     ylim = c(0,1.15))
points(Cv_b.s, fplus.s)

mtext(expression(italic(f)+""), side = 2, line = 3, cex = 1.25)

lines(loess.smooth(Cv_b.b, fplus.b), lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(Cv_b.s, fplus.s), lwd = 2, col = "seagreen", lty = 6)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(c)", cex = 1.25)
par(usr = usr0)

#   legend("bottomright",
#          c("Susp. Sites", "Bed Sites", "LOESS Susp.", "LOESS Bed"),
#          pch = c(1,8, NA, NA),
#          pt.cex = 1.2,
#          lty = c(NA, NA, 6, 2),
#          col = c("black", "black", "seagreen", "royalblue"),
#          bty = "n",
#          cex = 0.8,
#          lwd = c(1,1,2,2),
#          seg.len = 4)

# CV_Beta v.s Spread #
plot(Cv_b.b, spread.b, pch = "*", cex = 1.5,
     ylim = range(c(spread.b, spread.s)),
     xlim = range(c(Cv_b.b, Cv_b.s)),
     xlab = "",
     ylab = "",
     cex.lab = 1.25,
     yaxt = "n"
)
points(Cv_b.s, spread.s)

axis(side = 2, at = 1:3)
axis(side = 2, at = 0.5:3.5, labels = F, tcl = -0.2)

lines(loess.smooth(Cv_b.b, spread.b), lwd = 2, col = "royalblue", lty = 2)
lines(loess.smooth(Cv_b.s, spread.s), lwd = 2, col = "seagreen", lty = 6)

usr0 <- par("usr")       # get user coordinates
par(usr = c(0, 1, 0, 1))
text(.1, .9, labels = "(f)", cex = 1.25)
par(usr = usr0)

mtext(expression(italic(C[v]/beta)), side = 1, line = 3, cex = 1.25)
mtext("Yield Spread", side = 2, line = 3, cex = 1.25)

dev.off()