################################################################
## Comparison of Measured v. Modeled Sediment Transport
## for coarse bed sites using bed load transport models.
##
## Joel Sholtes, PhD Candidate
## Colorado State University, Dept. of Civil and Env. Engr.
## April 2015
###############################################################

rm(list = ls())
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/transport")
library(segmented) # Piece-wise regression
library(MASS)

# Specify row index of site
k = 57 # Trapper Ck
k = 48 # SF Salmon
name <- "sfsalmon"# "trapper" 

# Import Sed Data
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload")

## Import site summary data
path <- "site_summary.csv"
sites <- read.table(path, header = TRUE, sep = ",", 
                         stringsAsFactors = F)
sites <- sites[order(flow.sites$site.no),]

# Read in Bed Load Sed Data
sed.data <- read.csv("BedLoad_R_Data.csv", stringsAsFactors = F, header = TRUE)

# Extract Sed Data
site.data <- subset(sed.data, site.id == sites$site.no[k])
kg.sec <- site.data$kg.sec
cms <- site.data$cms

# Import Flow data
path <- paste("DV_Flows/", sites$site.no[k], ".csv", sep = "")
flow.data <- read.csv(path, header = T, stringsAsFactors = F)
flow <- flow.data$val*(.3048)^3
Qmin <- min(flow)
Qmax <- max(flow)

############ Qeff Analysis ################
npts = 50
xpts <- exp( seq( log(min(flow) + 0.001), 
                  log(max(flow) - 0.001), length = npts) )

# create empirical CDF using flow record
ecdf.prob <- ( 1:length(flow) ) / length(flow)
flow.sort <- sort(flow)

# create sequence of interpolated eCDF values evaluated at 'xpts1'
ypts <- approx(flow.sort, ecdf.prob, xout = xpts)$y

# add extreme from sorted flow and eCDF to vectors for differentiation
xpts <- c(min(flow), xpts, max(flow))
ypts <- c(min(ecdf.prob), ypts, max(ecdf.prob))

par(mfrow = c(1,1))
plot(flow.sort, ecdf.prob, type = "l")
lines(xpts, ypts, col = "blue")
points(xpts, ypts, col = "darkblue")

# Numerical differentiation of eCDF for Qeff analysis
# 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.e <- numeric(length(xpts) - 2)
n <- length(xpts) - 1
q.e <- xpts[2:n]
for (j in 2:n) { # j starts at 2, but assigns central diff to PDF at j-1
  pdf.e[j-1] <- ( ypts[j+1] - ypts[j-1] ) / ( xpts[j+1] - xpts[j-1] )
}

############ Hydraulic Analysis #############
# Specify Manning's n (one for entire XS)
# could also be based on grain size...
n = 0.04 #0.04 %0.055 trapper
# Specify Slope
s <- sites$slope[k]
#s <- .03

# Read in HG data and XS survey data
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/transport")
xs <- read.csv(paste(name, "XS.csv", sep = ""), header = T, stringsAsFactors = F)
hg <- read.csv(paste(name, "HG.csv", sep = ""), header = T, stringsAsFactors = F)

plot(xs$x.m, xs$y.m, type = "l")
xs$x.m <- xs$x.m + 1
y.range <- c(min(xs$y.m), min(xs$y.m[1], xs$y.m[nrow(xs)]))

########################################################
## Optional: Add verticla wall to lower side
y.range <- range(xs$y.m)

# Add vertical wall to XS if one side has lower elevation than other
if(xs$y.m[1] < xs$y.m[nrow(xs)]){
  xs <- rbind(c(min(xs$x.m), xs$y.m[nrow(xs)]),xs)
}else if(xs$y.m[nrow(xs)] < xs$y.m[1]){
  xs <- rbind(xs,c(max(xs$x.m), xs$y.m[1]))
}
       
plot(xs$x.m, xs$y.m, type = "l")
###########################################################

# create vector of depths to develop a rating curve
d.tot <- y.range[2] - y.range[1]
depth <- seq(.01,d.tot-d.tot*.001, by = 0.05)

# set up vectors and matrices for Manning's n rating curve
idx = 2:length(xs$y.m)
Q <- numeric(length(depth))
d.avg <- numeric(length(depth))
r.avg <- numeric(length(depth))
h <- numeric(length(depth))
d.trap <- matrix(ncol = nrow(xs)-1, nrow = length(depth))
w.trap <- matrix(ncol = nrow(xs)-1, nrow = length(depth))
a.avg <- numeric(length(depth))
xi <- matrix(ncol = nrow(xs)-1, nrow = length(depth) )

for(i in length(depth):1){
  
  elv <- y.range[1]+depth[i]
  d <- elv - xs$y.m
  
  # Calc extent of top width at elv
  x.i <- ifelse(d[idx-1] < 0, # if previous depth < 0
             ifelse(d[idx] < 0, # AND if current depth < 0
                    0, 
                    # if current depth > 0 AND previous depth < 0, 
                    # calc x using interp prev x value 
                   xs$x.m[idx-1] + (xs$x.m[idx]-xs$x.m[idx-1]) / (xs$y.m[idx]-xs$y.m[idx-1]) *
                   (elv - xs$y.m[idx-1])
                   ),
             ifelse(d[idx] < 0, # if prev depth > 0 and current depth < 0
                   xs$x.m[idx-1] + (xs$x.m[idx] - xs$x.m[idx-1])/(xs$y.m[idx]-xs$y.m[idx-1]) *
                          (elv-xs$y.m[idx-1]),
                    # if prev depth > 0 and current depth > 0 use trapezoid area
                   xs$x.m[idx]
             )
            )
  # Calc flow area at elv
  a <-  sum(ifelse(d[idx-1] < 0, # if previous depth < 0
              ifelse(d[idx] < 0, # AND if current depth < 0
                0, 
                # if current depth > 0 AND previous depth < 0, 
                # calc area using interp prev x value (triangle area)
                0.5*(elv - xs$y.m[idx])*
                    (xs$x.m[idx] - (xs$x.m[idx-1] + 
                      (xs$x.m[idx]-xs$x.m[idx-1]) / (xs$y.m[idx]-xs$y.m[idx-1]) *
                      (elv - xs$y.m[idx-1]))
                    )
              ),
              ifelse(d[idx] < 0, # if prev depth > 0 and current depth < 0
                0.5*(elv - xs$y.m[idx-1])*
                    ((xs$x.m[idx-1] + (xs$x.m[idx] - xs$x.m[idx-1])/(xs$y.m[idx]-xs$y.m[idx-1]) *
                     (elv-xs$y.m[idx-1]))-xs$x.m[idx-1]),
                # if prev depth > 0 and current depth > 0 use trapezoid area
                (xs$x.m[idx]-xs$x.m[idx-1])*(d[idx] + d[idx-1])/2
              )
            )
          )
  
  # Calc wetted perimeter at elv
  p <-  sum(ifelse(d[idx-1] < 0, # if previous depth < 0
                   ifelse(d[idx] < 0, # AND if current depth < 0
                          0, 
                          # if current depth > 0 AND previous depth < 0, 
                          # calc dist using interp prev x value using Pythagoreom
                          sqrt(( xs$x.m[idx] - (xs$x.m[idx-1] + (xs$x.m[idx]-xs$x.m[idx-1]) / 
                                                  (xs$y.m[idx]-xs$y.m[idx-1]) * (elv-xs$y.m[idx-1])) )^2 +
                                 (elv-xs$y.m[idx])^2
                          )
                   ),
                   ifelse(d[idx] < 0, # if prev depth > 0 and current depth < 0
                          # calc dist using interp prev x value using Pythagoreom
                          sqrt(((xs$x.m[idx]-xs$x.m[idx-1]) / (xs$y.m[idx]-xs$y.m[idx-1]) *
                                  (elv-xs$y.m[idx-1]))^2 + (elv-xs$y.m[idx-1])^2
                          ),
                          # if prev depth > 0 and current depth > 0 use pythag with XS pts only
                          sqrt((xs$x.m[idx]-xs$x.m[idx-1])^2 + (d[idx] - d[idx-1])^2)
                   )
                )
              )
  # Calc avg. depth at each midpoint at elv
  d.trap[i,] <- ifelse(d[idx-1] < 0, 
                 ifelse(d[idx] < 0, 
                    # if previous depth < 0 AND current depth < 0
                    0, 
                    # if current depth > 0 AND previous depth < 0, 
                    d[idx]/2
                  ),
                  ifelse(d[idx] < 0, 
                    # if prev depth > 0 and current depth < 0
                    d[idx-1]/2,
                    # if prev depth > 0 and current depth > 0 
                    (d[idx] + d[idx-1])/2
                  )
                )
  
  # Calc width over interval with depth d.trap at elv
  w.trap[i,] <-  ifelse(d[idx-1] < 0, # if previous depth < 0
                   ifelse(d[idx] < 0, # AND if current depth < 0
                          0, 
                          # if current depth > 0 AND previous depth < 0, 
                          # calc width using interp prev x value 
                          xs$x.m[idx] - (xs$x.m[idx-1] + (xs$x.m[idx]-xs$x.m[idx-1]) / 
                          (xs$y.m[idx]-xs$y.m[idx-1])*(elv - xs$y.m[idx-1])) 
                        ),
                   ifelse(d[idx] < 0, 
                          # if prev depth > 0 and current depth < 0
                          # calc width using interp prev x value 
                          (xs$x.m[idx] - xs$x.m[idx-1]) / 
                          (xs$y.m[idx] - xs$y.m[idx-1])*(elv - xs$y.m[idx-1]),
                          # if prev depth > 0 and current depth > 0 
                          xs$x.m[idx] - xs$x.m[idx-1]
                        )
                  )
  r <- a/p
  a.avg[i] <- a
  r.avg[i] <- r
  Q[i] <- sqrt(s)*a*r^(2/3)/n
  d.avg[i] <- a/(max(x.i)-min(x.i[which(x.i>0)]))
  h[i] <- elv
  xi[i,] <- x.i
}

w.avg <- a.avg/d.avg
w.top <- apply(w.trap, MARGIN = 1, FUN = sum)

xs.hydr <- data.frame(d.avg = d.avg, maxd = h-y.range[1], Q = Q, v.avg = Q/a.avg, a.avg = a.avg, p = p, r.avg = r.avg, w.avg = a.avg/d.avg, w.top = w.top)
#xs.hydr <- xs.hydr[-nrow(xs.hydr),] 

### PLOT XS ####
par(mfrow = c(1,1), pty = "m")
plot(xs$x.m, xs$y.m, type = "l", col = "grey30", lwd = 2,
     xlab = "Station (m)",
     ylab = "Elevation (m)")
for (i in 1:length(Q)){
  x.wse <- xi[i,]
  x.wse <- x.wse[which(x.wse !=0)]
  lines(x.wse, rep(xs.hydr$maxd[i], length(x.wse)) + y.range[1], col = "royalblue")
}

polygon( c(xs$x.m, rev(xs$x.m)), c(xs$y.m,rep(0,nrow(xs))), col="white", border=F) 
lines(xs$x.m, xs$y.m, col = "grey30", lwd = 2)
box()

par(mfrow = c(1,1), mar = c(4.2, 4.5, 1, 1))
plot(Q,depth, log = "xy", type = "l", lwd = 2, col = "gray30", 
#      xlim = c(1, 500),
#      ylim = c(.1, 10),
     xlab = expression("Discharge ("*m^3/s*")"),
     ylab = "Average Depth (m)",
     cex.lab = 1.2)
points(hg$Q.cms, hg$h.m,  col ="blue")
txt <- sites$site.nm[k]
text(x=grconvertX(.05, from='npc'), y = grconvertY(.90, from='npc'), labels = txt, cex = 0.9, adj = c(0,0))

# Create at-a-station HG relation using data
n.data <- nrow(hg)
plot(hg$Q.cms, hg$h.m,log = "xy")
lnQ <- log(hg$Q.cms)
lnh <- log(hg$h.m)
hQ.fit <- lm(lnh ~ lnQ)
bcf <- exp(sum( residuals(hQ.fit)^2 )/(n.data-2) /2)
c <- exp(coef(hQ.fit)[1])
f <- coef(hQ.fit)[2]
Q.range <- seq(min(hg$Q.cms), max(hg$Q.cms), length.out = 100)
d.average <- c*Q.range^f
lines(Q.range, d.average)
d.fun <- function(Q){c*bcf*Q^f}
lines(Qrange, d.fun(Qrange), col = "blue")
#   # Conduct segmented regression using best guess values of upper break point
#   # NOTE: This is conducted site-by-site and not in the loop. 
#   # Need to first find the slopes and breaks, add them to the seg.data table and use the 
#   # seg.data table in the loop to calc the stats and extend the flow records (2 steps)
#   # Need to repeat the segmented regression until a visually good piecewise regression results
# 
# seg.lm <- segmented( hQ.fit, seg.Z =~ lnQ, psi = list(lnQ = log(150)) )
#   plot.segmented(seg.lm)
#   points(lnQ, lnh)
#   int1 <- intercept(seg.lm)$lnQ[1]
#   int2 <- intercept(seg.lm)$lnQ[2]
#   slp1 <- slope(seg.lm)$lnQ[1]
#   slp2 <- slope(seg.lm)$lnQ[2]
#   break.Q <- exp(summary.segmented(seg.lm)$psi[2])
#   break.h <- exp(int1 + slp1*log(break.Q))
# 
# plot(hg$Q.cms, hg$h.m,log = "xy")
# lines(Q.range[which(Q.range<=break.Q)],exp(int1)*Q.range[which(Q.range<=break.Q)]^slp1, col = "blue", lwd = 2)
# lines(Q.range[which(Q.range>break.Q)], exp(int2)*Q.range[which(Q.range>break.Q)]^slp2, col = "blue", lwd = 2)
# 
# # create function for piece-wise rating curve
# # int1 and int2 have not been exponentiated
# d.fun <- function(Q) {
#   ifelse(Q <= break.Q, exp(int1)*Q^slp1*bcf, exp(int2)*Q^slp2*bcf)
# }
# 
# lines(Q.range,d.fun(Q.range))

## Create segmented model for width to dimensionalize Qs
plot(hg$Q.cms, hg$w.m, log = "xy")
# For SF salmon, need to remove some spurios data from either U/S or D/S
idx.w <- which(hg$w.m > 31.5 & hg$Q.cms < 19 ) # SF Salmon
#idx.w <- which(hg$w.m > 6) # Trapper
w.hg <- hg$w.m[-idx.w]
q.hg <- hg$Q.cms[-idx.w]
points(q.hg,w.hg, col = "green")

lnQ <- log(q.hg)
lnW <- log(w.hg)
wQ.fit <- lm(lnW ~ lnQ)
plot(lnQ, lnW)
abline(coef(wQ.fit)) 
Wbcf <- exp(sum( residuals(wQ.fit)^2 )/(n.data-2) /2)

seg.lm <- segmented( wQ.fit, seg.Z =~ lnQ, psi = list(lnQ = 1) )
plot.segmented(seg.lm)
points(lnQ, lnW)
Wint1 <- intercept(seg.lm)$lnQ[1]
Wint2 <- intercept(seg.lm)$lnQ[2]
Wslp1 <- slope(seg.lm)$lnQ[1]
Wslp2 <- slope(seg.lm)$lnQ[2]
Wbreak.Q <- exp(summary.segmented(seg.lm)$psi[2])
break.W <- exp(Wint1 + Wslp1*log(Wbreak.Q))

# create function for piece-wise rating curve for width
# Wint1 and Wint2 have not been exponentiated
w.fun <- function(Q) {
  ifelse(Q <= Wbreak.Q, exp(Wint1)*Q^Wslp1*Wbcf, exp(Wint2)*Q^Wslp2*Wbcf) # Trapper
}

plot(hg$Q.cms, hg$w.m)
lines(Q.range,w.fun(Q.range), col = "green")
lines(loess.smooth(q.hg, w.hg), col = "blue")


########## HG Plots ##################
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/transport/figures")

pdf("sfsalmonHG2.pdf", width = 8, height = 3.5)

par(mfrow = c(1,3), mar = c(5,5,1,1), pty = "s")

plot(hg$Q.cms, hg$stage.m, log = "xy", 
     xlab = expression("Discharge ("*m^3/s*")"),
     ylab = "Stage (m)",
     cex.lab = 1.2)
lines(Q, depth, lwd = 2, lty = 6, col = "royalblue")
txt <- "S.F. Salmon River"# Trapper Creek, ID"" #Indian Ck near Wyoming, IL"
text(x=grconvertX(.05, from='npc'), y = grconvertY(.9, from='npc'), labels = txt, cex = 0.9, adj = c(0,0))

plot(hg$Q.cms, hg$h.m,log = "xy",
     ylab = "Average Depth (m)",
     xlab = expression("Discharge ("*m^3/s*")"),
     cex.lab = 1.2,
     ylim = c(.05, 5),
     xlim = c(1, 200))
lines(Q.range,d.fun(Q.range), col = "darkcyan", lwd = 2)
lines(xs.hydr$Q, xs.hydr$d.avg, col = "darkcyan", lwd = 2, lty = 6)

legend("bottomright",
       c("Log-Linear Regression", "Manning's XS"),
       lty = c(1,6),
       lwd = 2, 
       col = "darkcyan",
       bty = "n",
       cex = 0.75)

plot(q.hg, w.hg, log = "xy",
     ylab = "Top Width (m)",
     xlab = expression("Discharge ("*m^3/s*")"),
     cex.lab = 1.2)
lines(Q.range,w.fun(Q.range), col = "brown3", lwd = 2)
lines(xs.hydr$Q, xs.hydr$w.top, col = "brown3", lwd = 2, lty = 6)

dev.off()

#### Now Calc Sediment Transport ####

#Define sediment sizes (D50)
d50 <- sites$d50[k]/1000

# Range of flows and average depths and widths
Qrange <- seq(Qmin-.01, Qmax+.5, by = 0.1)
depths <- d.fun(Qrange)
widths <- w.fun(Qrange)

## 1) HG relation
tauC <- 0.03
tauStar <- s*depths / (1.65*d50)
qStar <- 11.2*(tauStar - tauC)^4.5 / tauStar^3
# Remove NaN's
qStar <- ifelse(is.nan(qStar), 0, qStar)
# Calc dimensionalizing constant
k1 <- widths*sqrt(1.65*9.81)*d50^1.5
Qs.hg <- qStar*k1

plot(Qrange, Qs.hg, log = "xy", ylim = c(1E-6, max(kg.sec)))
points(cms,kg.sec, pch = "*", col = "orange")

## 2) Avg Depth with Manning's n w/ Parker 1979
tauC <- 0.03
tauStar <- d.avg*s / (1.65*d50)
qStar <- 11.2*(tauStar - tauC)^4.5 / tauStar^3
qStar <- ifelse(is.nan(qStar), 0, qStar)
Qs.avg <- qStar*w.avg*sqrt(1.65*9.81)*d50^1.5
plot(Q, Qs.avg, log = "xy", ylim = c(1E-4,max(range(Qs.avg),range(kg.sec))), typ = "l")
points(cms,kg.sec, pch = "*", col = "orange")


## 3) Non-Avg Depth Manning's n w/ Parker 1979
tauC <- 0.03
tauStar <- d.trap*s / (1.65*d50)
qStar <- 11.2*(tauStar - tauC)^4.5 / tauStar^3
qStar <- ifelse(is.nan(qStar), 0, qStar)
Qs.trap <- qStar*w.trap*sqrt(1.65*9.81)*d50^1.5
Qs.trap <- apply(Qs.trap, MARGIN = 1, sum)
plot(Q, Qs.trap, log = "xy", ylim = c(1E-4,max(range(Qs.trap),range(kg.sec))), xlim = range(cms), type = "l")
points(cms,kg.sec, pch = "*", col = "orange")

## 4) Non-Avg Depth Manning's n w/ Wilcock and Kenworthy 2002
# Surface Based TauR

# Trapper Ck
Fs <- 0.05 # Fraction of sand in bed
Dg <- 82.5/1000
Ds <- 2/1000

# SF Salmon R
Fs <- .46 # Fraction of sand in bed
Dg <- 45/1000
Ds <- 1/1000

tauStarRg0 <- 0.035
tauStarRg1 <- 0.011
tauStarRs0 <- tauStarRg0*(Dg/Ds)
tauStarRs1 <- 0.065
tauStarRg <- tauStarRg1 + (tauStarRg0 - tauStarRg1)*exp(-14*Fs)
tauStarRs <- tauStarRs1 + (tauStarRs0 - tauStarRs1)*exp(-14*Fs)

tauRg <- tauStarRg*(1650)*9.81*Dg
tauRs <- tauStarRs*(1650)*9.81*Ds

# # Sub-surface-based
# Fs <- .12 # Fraction of sand in bed
# Dg <- 17
# Ds <- 1.1
# tauRg0 <- 0.043
# tauRg1 <- 0.008
# tauRs0 <- 1.8*tauRg0*(Dg/Ds)
# tauRs1 <- 0.045
# tauRg <- tauRg0 - (tauRg0 - tauRg1)/(1+exp(-25*(Fs-0.16)))
# tauRs <- tauRs0 - (tauRs0 - tauRs1)/(1+exp(-25*(Fs-0.16)))

# Convert depth to d' or the depth of flow acting on the bed grains
D90 <- 295/1000
ks <- 0.84*D90
U <- Q/a.avg
tol = 1E-7
N = 100
dPrime <- numeric(length(d.avg))
dp <- function(hStar) {2.5*log(11*hStar/ks) - U[i]/sqrt(9.81*hStar*s)}

for(i in 1:length(d.avg)){
  dPrime[i] <- uniroot(dp, c(0,max(d.avg)))$root
}

tau <- 9810*dPrime*s
uStar <- sqrt(tau/1000)

phi.g <- tau / tauRg
phi.s <- tau / tauRs

phiPrime <- 1.27

A <- 115 # coeff for WK 2002
chi <- 0.923 # numerator for WK 2002

wStar.g <- ifelse(phi.g < phiPrime,
                0.002*phi.g^7.5,
                A*(1-chi/phi.g^0.25)^4.5
               )
wStar.s <- ifelse(phi.s < phiPrime,
                  0.002*phi.s^7.5,
                  A*(1-chi/phi.s^0.25)^4.5
                )

# Sum gravel and sand Wi's
# unit width volumetric flux (m^3/s/m)
qs.s <- wStar.s * Fs * uStar^3  / (1.65 * 9.81)
qs.g <- wStar.g * (1-Fs) * uStar^3  / (1.65 * 9.81)

# mass flux: kg/s
Qs.s <- qs.s * 2650 * w.avg
Qs.g <- qs.g * 2650 * w.avg

Qs.wk <- Qs.s + Qs.g

plot(Q, Qs.wk, log = "xy", ylim = c(1E-5,max(range(Qs.wk),range(kg.sec))), xlim = range(cms), type = "l")
points(cms,kg.sec, pch = "*", col = "orange")
lines(Q, Qs.s, lty = 2, col = "blue")
lines(Q, Qs.g, lty = 2, col = "green")

## 5) Barry et al. 2004 general power law function and Manning's flow resitance

# SF Salmon 
DA = 853*1E6 # m^2
s = 0.0025 #m/m
d50ss = 14/1000 #m
d50s = 38/1000 #m
Q2 = 96.3 #cms

# Trapper
DA = 21*1E6 # m^2
s = 0.0414 #m/m
d50ss = 17/1000 #m
d50s = 75/1000 #m
Q2 = 2.21 #cms

tauC <- 0.03 # After Barry et al 2003
hQ2 <- d.avg[which.min(abs(Q-Q2))]
tauQ2 <- 9810*hQ2*s
tauD50s <- tauC*1650*9.81*d50s # tauC*(rho_s - rho)*g*d50s
tauD50ss <- tauC*1650*9.81*d50ss # tauC*(rho_s - rho)*g*d50ss

qStar <- (tauQ2 - tauD50s) / (tauQ2 - tauD50ss)

alpha = 8.13E-7*DA^0.49
beta <- -2.45*qStar+3.56
Qs.b <- alpha*(Q/Q2)^beta*w.avg # kg/s  width.trap ????
Qs.b.trap <- alpha*(Q/Q2)^beta*w.trap # kg/s  width.trap ????
Qs.b.trap <- apply(Qs.b.trap,1,sum)

plot(Q, Qs.b.trap, log = "xy", ylim = c(1E-5,max(range(Qs.b),range(kg.sec))), xlim = range(cms), type = "l")
points(cms,kg.sec, pch = "*", col = "orange")

Qs.b <- Qs.b.trap

## 6) Data Model
Q.sort <- sort(flow)

# Original RLM empirical model
# x-vector for linear regression mean response
x.range = range(log(flow))
x <- seq(x.range[1], x.range[2], length.out = 100) 

Qc <- mean(cms[which(kg.sec == 0)])

Qs.meas <- kg.sec[which(kg.sec != 0)]
Q.meas <- cms[which(kg.sec != 0)]
n.data = length(Qs.meas)

model.rlm <- rlm(log(Qs.meas) ~ log(Q.meas), method = "MM", maxit = 100)
bhat = coef(model.rlm)
plot(Q.meas, Qs.meas, log = "xy")
lines(Q, exp(bhat[1])*Q^bhat[2])
bcf <- signif(exp(sum( residuals(model.rlm)^2 )/(n.data-2) /2), 3)
yield.rlm <- exp(bhat[1])*bcf*flow^bhat[2]
plot(yield.rlm, type = "l", col = "chocolate")

yield.rlm.sort <- sort(yield.rlm)
yield.cum.rlm <- cumsum(yield.rlm.sort)/sum(yield.rlm.sort)

#### Sed Yield Calcs #####
eff.rlm <- exp(bhat[1])*bcf*q.e^bhat[2]*pdf.e
eff.trap <- approx(x = Q, y = Qs.trap, xout = q.e)$y*pdf.e
eff.avg <- approx(x = Q, y = Qs.avg, xout = q.e)$y*pdf.e
eff.hg <- approx(x = Qrange, y = Qs.hg, xout = q.e)$y*pdf.e
eff.b <- approx(x = Q, y = Qs.b, xout = q.e)$y*pdf.e
eff.wk <- approx(x = Q, y = Qs.wk, xout = q.e)$y*pdf.e

Qeff.rlm <- q.e[which.max(eff.rlm)]
Qeff.trap <- q.e[which.max(eff.trap)]
Qeff.avg <- q.e[which.max(eff.avg)]
Qeff.hg <- q.e[which.max(eff.hg)]
Qeff.b <- q.e[which.max(eff.b)]
Qeff.wk <- q.e[which.max(eff.wk)]

## Modeled sediment transport yield calc - by avg depth
yield <- approx(x = Q, y = Qs.avg, xout = flow)$y
plot(yield, type = "l", col = "chocolate")
plot(flow, type = "l", col = "blue")
yield.avg.sort <- sort(yield)
yield.cum.avg <- cumsum(yield.avg.sort)/sum(yield)

## Modeled sediment transport yield calc - by trapezoid
yield <- approx(x = Q, y = Qs.trap, xout = flow)$y
plot(yield, type = "l", col = "chocolate")
plot(flow, type = "l", col = "blue")
yield.trap.sort <- sort(yield)
yield.cum.trap <- cumsum(yield.trap.sort)/sum(yield)

## Modeled sediment transport yield calc - by HG relation
yield <- approx(x = Qrange, y = Qs.hg, xout = flow)$y
plot(yield, type = "l", col = "chocolate")
plot(flow, type = "l", col = "blue")
yield.hg.sort <- sort(yield)
yield.cum.hg <- cumsum(yield.hg.sort)/sum(yield)

## Modeled sediment transport yield calc - by Mannings and Barry et al. Power Law Qs
yield <- approx(x = Q, y = Qs.b, xout = flow)$y
plot(yield, type = "l", col = "chocolate")
plot(flow, type = "l", col = "blue")
yield.b.sort <- sort(yield)
yield.cum.b <- cumsum(yield.b.sort)/sum(yield)

## Modeled sediment transport yield calc - by Mannings and Wilcock and Kenworthy
yield <- approx(x = Q, y = Qs.wk, xout = flow)$y
plot(yield, type = "l", col = "chocolate")
plot(flow, type = "l", col = "blue")
yield.wk.sort <- sort(yield)
yield.cum.wk <- cumsum(yield.wk.sort)/sum(yield)

q50.rlm <- Q.sort[which.min(abs(yield.cum.rlm-0.5))]
q50.avg <- Q.sort[which.min(abs(yield.cum.avg-0.5))]
q50.trap <- Q.sort[which.min(abs(yield.cum.trap-0.5))]
q50.hg <- Q.sort[which.min(abs(yield.cum.hg-0.5))]
q50.b <- Q.sort[which.min(abs(yield.cum.b-0.5))]
q50.wk <- Q.sort[which.min(abs(yield.cum.wk-0.5))]

######################################################
## Plot it

# Sediment Rating Curves
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/transport/figures")
pdf("sfsalmon2.pdf", width = 8, height = 3)
  par(pty = "s", mar = c(4.1,4.6,.5,.5), mfrow = c(1,3))
  plot(cms,kg.sec, pch = "*", col = "brown2",
       ylab = "Sediment Load (kg/s)",
       xlab = expression("Discharge ("*m^3/s*")"),
       cex.lab = 1.2,
       log = "xy",
       yaxt = "n",
       xaxt = "n",
       ylim = c(min(kg.sec[which(kg.sec !=0)])*.1, max(kg.sec)*5),
       xlim = c(min(cms)*.1, max(cms)*5),
       cex = 1.5)
  lines(Q, Qs.avg, lty = 2, lwd = 1.5, col = "grey20") # XS-avg. Manning's depth-discharge 
  lines(Q, Qs.trap, lty = 3, lwd = 1.5, col = "grey30") # Non-avg. Manning's depth-discharge 
  lines(Qrange, Qs.hg, lty = 4, lwd = 1.5, col = "grey40") # HG relation
  lines(Q, Qs.b, lty = 5, lwd = 1.5, col = "grey50") # Barry 2004 + Q
  lines(Q, Qs.wk, lty = 6, lwd = 1.5, col = "grey60") # Non-avg. Manning's depth-discharge 
  lines(Q, exp(bhat[1])*Q^bhat[2], lty = 1, lwd = 1.5) # Regression Model

  yat <- floor(log10(min(kg.sec[which(kg.sec !=0)])*.1)):ceiling(log10(max(kg.sec)*5))
  labels <- sapply(yat, function(i) as.expression(bquote(10^ .(i))) )
  axis(side = 2, at = 10^yat, labels = labels)
  
  xat <- floor(log10(min(cms)*.1)):ceiling(log10(max(cms)*10))
  labels <- sapply(xat, function(i) as.expression(bquote(10^ .(i))) )
  axis(side = 1, at = 10^xat, labels = labels)
  
#   legend("topleft",
#          c("Log-Linear Regression",
#            "Parker79 Avg. XS Depth",
#            "Parker79 Discr. Depth",
#            "Parker79 HG Depth",
#            "Barry04",
#            "WK02"),
#          lwd = 1.5,
#          col = c("black", "grey20", "grey30", "grey40", "grey50", "grey60"),
#          lty = c(1,2,3,4,5,6),
#          y.intersp = 0.9,
#          #bty = "n",
#          #border = "n",
#          box.col = rgb(1,1,1,alpha = 0.75),
#          cex = 0.7,
#          seg.len = 3,
#          bg = rgb(1,1,1,alpha = 0.75))
  box()

  # Effective Discharge
  plot(q.e,eff.rlm, col = "brown2", type = "l",
       ylab = expression("Sediment Yield Density (kg/s/"*m^3/s*")"),
       xlab = expression("Discharge ("*m^3/s*")"),
       cex.lab = 1.2,
       lwd = 2,
       log = "y",
       ylim = c(1E-6, max(eff.wk)))
  lines(q.e, eff.avg, lty = 2, lwd = 1.5, col = "grey20") # XS-avg. Manning's depth-discharge 
  lines(q.e, eff.trap, lty = 3, lwd = 1.5, col = "grey30") # Non-avg. Manning's depth-discharge 
  lines(q.e, eff.hg, lty = 4, lwd = 1.5, col = "grey40") # HG relation
  lines(q.e, eff.b, lty = 5, lwd = 1.5, col = "grey50") # Barry 2004
  lines(q.e, eff.wk, lty = 6, lwd = 1.5, col = "grey60") # WK 2002
  segments(q.e[which.max(eff.rlm)], 1E-10, q.e[which.max(eff.rlm)], max(eff.rlm), col = "brown2")
  segments(q.e[which.max(eff.avg)], 1E-10, q.e[which.max(eff.avg)], max(eff.avg),col = "grey20", lty = 2)
  segments(q.e[which.max(eff.trap)], 1E-10, q.e[which.max(eff.trap)], max(eff.trap), col = "grey30", lty =3)
  segments(q.e[which.max(eff.hg)], 1E-10, q.e[which.max(eff.hg)], max(eff.hg), col = "grey40", lty = 4)
  segments(q.e[which.max(eff.b)], 1E-10, q.e[which.max(eff.b)], max(eff.b), col = "grey50", lty = 5)
  segments(q.e[which.max(eff.wk)], 1E-10, q.e[which.max(eff.wk)], max(eff.wk), col = "grey60", lty = 6)
  
# Cum Sed Yield
  plot(Q.sort, yield.cum.rlm, type = "l",
       lwd = 2,
       ylab = "Cum. Sediment Yield",
       xlab = expression("Discharge ("*m^3/s*")"),
       cex.lab = 1.2,
       col = "brown2")
  lines(Q.sort, yield.cum.avg, lty = 2, lwd = 1.5, col = "grey20") # XS-avg. Manning's depth-discharge 
  lines(Q.sort, yield.cum.trap, lty = 3, lwd = 1.5, col = "grey30") # Non-avg. Manning's depth-discharge 
  lines(Q.sort, yield.cum.hg, lty = 4, lwd = 1.5, col = "grey40") # HG relation
  lines(Q.sort, yield.cum.b, lty = 5, lwd = 1.5, col = "grey50") # Barry 2004
  lines(Q.sort, yield.cum.wk, lty = 6, lwd = 1.5, col = "grey60") # WK 2002

  segments(Q.sort[which.min(abs(yield.cum.rlm-0.5))], -10, Q.sort[which.min(abs(yield.cum.rlm-0.5))], 0.5,
           col = "brown2")
  segments(Q.sort[which.min(abs(yield.cum.avg-0.5))], -10, Q.sort[which.min(abs(yield.cum.avg-0.5))], 0.5,
          col = "grey30", lty = 2)
  segments(Q.sort[which.min(abs(yield.cum.trap-0.5))], -10, Q.sort[which.min(abs(yield.cum.trap-0.5))], 0.5,
           col = "grey50", lty = 3)
  segments(Q.sort[which.min(abs(yield.cum.hg-0.5))], -10, Q.sort[which.min(abs(yield.cum.hg-0.5))], 0.5,
           col = "grey70", lty = 4)
  segments(Q.sort[which.min(abs(yield.cum.b-0.5))], -10, Q.sort[which.min(abs(yield.cum.b-0.5))], 0.5,
           col = "grey70", lty = 5)
  segments(Q.sort[which.min(abs(yield.cum.wk-0.5))], -10, Q.sort[which.min(abs(yield.cum.wk-0.5))], 0.5,
           col = "grey70", lty = 6)
  segments(-10, 0.5, Q.sort[which.min(abs(yield.cum.rlm-0.5))], 0.5, col = "brown2")
  txt <- sites$site.nm[k]
  text(x=grconvertX(.01, from='npc'), y = grconvertY(.95, from='npc'), labels = txt, cex = 0.9, adj = c(0,0))

dev.off()
