################################################################
## Comparison of Measured v. Modeled Sediment Transport
## for fine bed sites using total 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)

# Import Susp Site Info, Sed, and Flow Data
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Suspload")

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

# Specify row index of site
k = 34 # indian ck
name <-  "yampa" # "indian" # 
s = 0.003
# Read in Bed Load Sed Data
sed.data <- read.table("susp_data_clean_outliers.csv", stringsAsFactors = F,header = TRUE, sep = ",",
                       colClasses = c("character", "character", "character", "numeric", 
                                      "numeric", "numeric","numeric")
)

sed.data$SITE_NO <- as.integer(sed.data$SITE_NO)

# Yampa sediment load data imported separately. For some reason not part of 
# original NWIS suspended load database
yampa.sed <- read.csv("yampaSed.csv", header = T, stringsAsFactors = F)
sed.site <- yampa.sed

# Extract Sed Data
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Suspload")
site <- sites.flow[k]
site.data <- subset(sed.data, SITE_NO == site)
# kg.sec <- site.data$kg.sec
# cms <- site.data$cms

x <- site.data$INSTFLOW * 0.0283168466 # convert CFS to CMS
q.min.sed <- signif(min(x, na.rm= T),2)
q.max.sed <- round(max(x, na.rm= T),2)

# here I am converting from [mg/L] SSC to mg/L of susp. sand (P63 is % finer than 0.0625 mm)
# 0.001 converts to kg/s: [mg/L*m^3/s * 1000L/m^3 * 1kg/1E*6mg = kg/s] 
y <- site.data$SSC * (1 - site.data$P63/100) * x * 0.001

# Write data to table along with CI for regression line
sed.site <- data.frame(cms = x, kg.sec = y, log.q = log10(x), log.qs = log10(y))

# now remove zero values for log-linear regression
if (min(y) == 0) {
  sed.site <- sed.site[which(y != 0),]
}


# Import Flow data
path <- paste("DV_Flows2/0", site, ".csv", sep = "")
#path <- paste("DV_Flows2/0","9260050", ".csv", sep = "") # Yampa
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.035
# Specify Slope
# calculated by bathymetric DEM from Cory Williams USGS
#s <- 0.003 #sites$slope[k]

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

xs$x.m <- xs$x.m + 1

plot(xs$x.m, xs$y.m, type = "l")

y.range <- c(min(xs$y.m), min(xs$y.m[1], xs$y.m[nrow(xs)]))

# Optional:
# 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]))
}
y.range <- range(xs$y.m)

plot(xs$x.m, xs$y.m, type = "l")

# create vector of depths to develop a rating curve
depth <- seq(.01,y.range[2]-y.range[1], 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
}

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)

############## XS Plots #############

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, 4000),
     ylim = c(.1, 10),
     xlab = expression("Discharge ("*m^3/s*")"),
     ylab = "Stage (m)",
     cex.lab = 1.2)
points(hg$Q.cms, hg$stage.m,  col ="blue")
txt <- "Yampa River at Deerlodge" #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 for depth 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}


  # 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
plot(lnQ,lnh)

seg.lm <- segmented( hQ.fit, seg.Z =~ lnQ, psi = list(lnQ = -1) )
  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), col = "green", lwd = 2)

# Create at-a-station HG relation for area using data
plot(hg$Q.cms, hg$a.m2, log = "xy")
lnQ <- log(hg$Q.cms)
lna <- log(hg$a.m2)
aQ.fit <- lm(lna ~ lnQ)
bcf.a <- exp(sum( residuals(hQ.fit)^2 )/(n.data-2) /2)
k <- exp(coef(aQ.fit)[1])
m <- coef(aQ.fit)[2]
Q.range <- seq(min(hg$Q.cms), max(hg$Q.cms), length.out = 100)
a.average <- k*Q.range^m
lines(Q.range, a.average)
a.fun <- function(Q){k*bcf.a*Q^m}

## 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 > 6)
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)
a.w <- exp(coef(wQ.fit)[1])
b.w <- coef(wQ.fit)[2]
plot(lnQ, lnW)
abline(coef(wQ.fit)) 
Wbcf <- exp(sum( residuals(wQ.fit)^2 )/(n.data-2) /2)
w.fun <- function(Q){a.w*Wbcf*Q^b.w}
plot(hg$Q.cms, hg$w.m, log = "xy")
lines(Qrange, w.fun(Qrange))

#### Segmented Regression for Width if needed
plot(lnQ,lnW)
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)
}

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

# plot(hg$Q.cms, hg$v.ms, log = "xy")

########## HG Plots ##################
pdf("YampaHG2.pdf", width = 6, height = 6)

  par(mfrow = c(2,2), 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 <- "Yampa River at Deerlodge, CO" #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(.005, 10))
  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(hg$Q.cms, hg$w.m, 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.avg, col = "brown3", lwd = 2, lty = 6)
  
  plot(hg$Q.cms, hg$a.m2, log = "xy",
       ylab = expression("Flow Area ("*m^2*")"),
       xlab = expression("Discharge ("*m^3/s*")"),
       cex.lab = 1.2)
  lines(Q.range,a.fun(Q.range), col = "slateblue", lwd = 2)
  lines(xs.hydr$Q, xs.hydr$a.avg, col = "slateblue", lwd = 2, lty = 6)

dev.off()

#### Now Calc Sediment Transport ####
FlowVel = xs.hydr$v.avg
q = xs.hydr$Q/xs.hydr$w.avg
  
## 1) Brownlie Bed Material Load Concentration

# Yampa Median Values
D16 = 0.27/1000 # m
D50 = 0.41/1000 # m
D84 = 0.68/1000 # m
slope = 0.003

# Indian Ck Median Values
D16 = 0.27/1000 # m
D50 = 1/1000 # m
D84 = 5.6/1000 # m
slope = 0.001

rho = 998.2 #'kg/m3 at T=20°C
t0 = 273.16 #'K
mu0 = 0.001792 #'kg/(m-s) = Ns/m2 = FLT-2
t = 20 # temp in C
t = t + t0
a = -1.94
b = -4.8
C = 6.74
# dynamic viscosity
mu = exp(a + b * (t0 / t) + C * (t0 / t) ^ 2 + log(mu0))

# kinematic viscosity
nu = mu / rho

SG = 2.65
g = 9.81 # m/s^2
BrownlieCoef = 1.268
GrReynolds = sqrt(g * D50 ^ 3) / nu
Y = (GrReynolds * sqrt(SG - 1)) ^ -0.6
DmCritShear = 0.22 * Y + 0.06 * 10 ^ (-7.7 * Y)
SigmaG = sqrt(D84 / D16)

# put an upper limit on dbSigmaG
if(SigmaG > 5) SigmaG = 5

GrFroudeO = 4.596 * DmCritShear ^ 0.5293 * slope ^ -0.1405 * SigmaG ^ -0.1606

GrFroude = FlowVel / sqrt((SG - 1) * g * D50)
GrFroudePrime = 1.74 * slope ^ -(1 / 3)
RatCoef = GrFroude / GrFroudePrime

## Need to choose one here... ########################

#If in stand alone configuration, relative roughness
#can be calculated with the Depth Predictor equations as
#follows:
RelRoughness <- ifelse(RatCoef < 1.25 | slope > 0.006,
                       0.3742 * ( q / sqrt(g * D50 ^ 3)) ^ 0.6539 * slope ^ -0.2542 * SigmaG ^ 0.105,
                       0.2836 * ( q / sqrt(g * D50 ^ 3)) ^ 0.6248 * slope ^ -0.2877 * SigmaG ^ 0.08013
                       )
#Or else relative roughness can be calculated as follows if
#the user has explicitly entered a hydraulic radius
  HydRad = xs.hydr$r.avg
  RelRoughness = HydRad / D50
######################################################

GrFroudeDiff <- GrFroude - GrFroudeO

Cppm <- ifelse(GrFroudeDiff <= 0, 
               0,
               round(7115 * BrownlieCoef * (GrFroudeDiff) ^ 1.978 * slope ^ 0.6601 * RelRoughness ^ -0.3301, 4)
              )

# convert ppm to mg/L
Cmg.l <- SG * Cppm / (SG + (1 - SG)*10^-6 * Cppm)

#Yield in kg/sec: 0.001 converts to kg/s: [mg/L*m^3/s * 1000L/m^3 * 1kg/1E*6mg = kg/s] 
Qs.br = Cmg.l * Q * 0.001

plot(Q,Qs.br, type = "l", log = "xy")
# points(sed.site$cms, sed.site$kg.sec, log = "xy", ylim = c(1E-5, 1), xlim = c(1E-3, 10))
points(sed.site$cms, sed.site$kg.sec, col = "blue")
# Brownlie over-predicts by quite a lot... data line up much better if D50 = 10 * D50

## 2) Yang D50 sand

#Calculate the terminal particle fall velocity via Rubey's formula
fRubey = (2 / 3 + 36 * nu ^ 2 / (g * D50 ^ 3 * (SG - 1))) ^ 0.5 - (36 * nu ^ 2 / (g * D50 ^ 3 * (SG - 1))) ^ 0.5
FallVel = fRubey * (D50 * g * (SG - 1)) ^ 0.5

# # Eqn 5.23d from Julien 2010
# dStar <- D50 * ( (SG - 1)*g / nu^2)^(1/3)
# FallVel = 8 * nu / D50 *( (1 + 0.0139*dStar^3)^0.05 - 1)

#Calculate parameters for Yang's Equation
UStar = sqrt(g * HydRad * slope)
UnitSP = FlowVel * slope / FallVel

#Determine Critical Unit Stream power based on previously calculated
#parameters.
CritUnitSP <- ifelse(UStar * D50 / nu > 1.02 & UStar * D50 / nu < 70,
                     slope *( 2.5 / (log10(UStar * D50 / nu) - 0.06) + 0.66),
                     ifelse(UStar * D50 / nu >= 70,
                            2.05 * slope,
                            0
                        )
                      )

#The following code calculates the log of the fine sediment concentration.
#A base 10 exponentiation is then used to determine the fine sediment
#concentration in PPM.
CritUnitSPdiff <- UnitSP - CritUnitSP

logCppm <- ifelse(CritUnitSPdiff > 0,
                  5.435 - 0.286 * log10(FallVel * D50 / nu) - 
                  0.457 * log10(UStar / FallVel) + (1.799 - 0.409 * log10(FallVel * D50 / nu) - 
                  0.314 * log10(UStar / FallVel)) * log10(UnitSP - CritUnitSP),
                  0
                  )
Cppm <- 10^logCppm
# convert ppm to mg/L
Cmg.l <- SG * Cppm / (SG + (1 - SG)*10^-6 * logCppm)

#Yield in kg/sec: 0.001 converts to kg/s: [mg/L*m^3/s * 1000L/m^3 * 1kg/1E*6mg = kg/s] 
Qs.yang = Cmg.l * Q * 0.001

lines(Q, Qs.yang, lty = 2)

###################################################
# Range of flows and average depths and widths
Qrange <- seq(Qmin+0.00001, Qmax+.5, by = 0.1)
depths <- d.fun(Qrange)
widths <- w.fun(Qrange)
areas <- a.fun(Qrange)
FlowVel = Qrange/areas
q <- Qrange / widths

## 3) HG with Brownlie
GrFroude = FlowVel / sqrt((SG - 1) * g * D50)
GrFroudePrime = 1.74 * slope ^ -(1 / 3)
RatCoef = GrFroude / GrFroudePrime

#If in stand alone configuration, relative roughness
#can be calculated with the Depth Predictor equations as
#follows:
RelRoughness <- ifelse(RatCoef < 1.25 | slope > 0.006,
                       0.3742 * ( q / sqrt(g * D50 ^ 3)) ^ 0.6539 * slope ^ -0.2542 * SigmaG ^ 0.105,
                       0.2836 * ( q / sqrt(g * D50 ^ 3)) ^ 0.6248 * slope ^ -0.2877 * SigmaG ^ 0.08013
)

GrFroudeDiff <- GrFroude - GrFroudeO

Cppm <- ifelse(GrFroudeDiff <= 0, 
               0,
               round(7115 * BrownlieCoef * (GrFroudeDiff) ^ 1.978 * slope ^ 0.6601 * RelRoughness ^ -0.3301, 4)
)

# convert ppm to mg/L
Cmg.l <- SG * Cppm / (SG + (1 - SG)*10^-6 * Cppm)

#Yield in kg/sec: 0.001 converts to kg/s: [mg/L*m^3/s * 1000L/m^3 * 1kg/1E*6mg = kg/s] 
Qs.br.hg = Cmg.l * Qrange * 0.001

plot(sed.site$cms, sed.site$kg.sec, log = "xy")

lines(Qrange, Qs.br.hg, lty = 1, col = "darkcyan")

## 4) Yang with HG
#Calculate parameters for Yang's Equation
UStar = sqrt(g * depths * slope)
UnitSP = FlowVel * slope / FallVel

#Determine Critical Unit Stream power based on previously calculated
#parameters.
CritUnitSP <- ifelse(UStar * D50 / nu > 1.02 & UStar * D50 / nu < 70,
                     slope *( 2.5 / (log10(UStar * D50 / nu) - 0.06) + 0.66),
                     ifelse(UStar * D50 / nu >= 70,
                            2.05 * slope,
                            0
                     )
)

#The following code calculates the log of the fine sediment concentration.
#A base 10 exponentiation is then used to determine the fine sediment
#concentration in PPM.
CritUnitSPdiff <- UnitSP - CritUnitSP

logCppm <- ifelse(CritUnitSPdiff > 0,
                  5.435 - 0.286 * log10(FallVel * D50 / nu) - 
                    0.457 * log10(UStar / FallVel) + (1.799 - 0.409 * log10(FallVel * D50 / nu) - 
                                                        0.314 * log10(UStar / FallVel)) * log10(UnitSP - CritUnitSP),
                  0
)
Cppm <- 10^logCppm
# convert ppm to mg/L
Cmg.l <- SG * Cppm / (SG + (1 - SG)*10^-6 * logCppm)

#Yield in kg/sec: 0.001 converts to kg/s: [mg/L*m^3/s * 1000L/m^3 * 1kg/1E*6mg = kg/s] 
Qs.yang.hg = Cmg.l * Qrange * 0.001

lines(Qrange, Qs.yang.hg, lty = 2, col = "darkgreen")


## 5) Data Model
Q.sort <- sort(flow)
cms <- sed.site$cms
kg.sec <- sed.site$kg.sec

# Original RLM empirical model
# x-vector for linear regression mean response
x.range = range(log(flow[which(flow != 0)]))
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)
lines(Qrange, exp(bhat[1])*Qrange^bhat[2], col = "red")

plot(Q.meas, Qs.meas, log = "xy")
lines(Qrange, exp(bhat[1])*Qrange^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.br <- approx(x = Q, y = Qs.br, xout = q.e)$y*pdf.e
eff.br.hg <- approx(x = Qrange, y = Qs.br.hg, xout = q.e)$y*pdf.e
eff.yang <- approx(x = Q, y = Qs.yang, xout = q.e)$y*pdf.e
eff.yang.hg <- approx(x = Qrange, y = Qs.yang.hg, xout = q.e)$y*pdf.e

Qeff.rlm <- q.e[which.max(eff.rlm)]
Qeff.br <- q.e[which.max(eff.br)]
Qeff.br.hg <- q.e[which.max(eff.br.hg)]
Qeff.yang <- q.e[which.max(eff.yang)]
Qeff.yang.hg <- q.e[which.max(eff.yang.hg)]

## Modeled sediment transport yield calc - by Brownlie and Manning's Vel
yield <- approx(x = Q, y = Qs.br, xout = flow)$y
plot(yield, type = "l", col = "chocolate")
plot(flow, type = "l", col = "blue")
yield.br.sort <- sort(yield)
yield.cum.br <- cumsum(yield.br.sort)/sum(yield)

## Modeled sediment transport yield calc - by Brownlie with HG Vel
yield <- approx(x = Qrange, y = Qs.br.hg, xout = flow)$y
yield[which(is.na(yield))] <- 0
plot(yield, type = "l", col = "chocolate")
plot(flow, type = "l", col = "blue")
yield.brHG.sort <- sort(yield)
yield.cum.brHG <- cumsum(yield.brHG.sort)/sum(yield)

## Modeled sediment transport yield calc - by Yang with Manning's Vel and Depth
yield <- approx(x = Q, y = Qs.yang, xout = flow)$y
plot(yield, type = "l", col = "chocolate")
plot(flow, type = "l", col = "blue")
yield.yang.sort <- sort(yield)
yield.cum.yang <- cumsum(yield.yang.sort)/sum(yield)

## Modeled sediment transport yield calc - by Yang with HG
yield <- approx(x = Qrange, y = Qs.yang.hg, xout = flow)$y
yield[which(is.na(yield))] <- 0
plot(yield, type = "l", col = "chocolate")
plot(flow, type = "l", col = "blue")
yield.yangHG.sort <- sort(yield)
yield.cum.yangHG <- cumsum(yield.yangHG.sort)/sum(yield)

q50.rlm <- Q.sort[which.min(abs(yield.cum.rlm-0.5))]
q50.br <- Q.sort[which.min(abs(yield.cum.br-0.5))]
q50.brGH <- Q.sort[which.min(abs(yield.cum.brHG-0.5))]
q50.yang <- Q.sort[which.min(abs(yield.cum.yang-0.5))]
q50.yangHG <- Q.sort[which.min(abs(yield.cum.yangHG-0.5))]

## Plot it

# Sediment Rating Curves
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Uncertainty/transport/suspended/figures")
pdf("YampaSed.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 = "darkcyan",
       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 = 2)
  lines(Q, Qs.br, lty = 2, lwd = 1.5, col = "grey20") 
  lines(Qrange, Qs.br.hg, lty = 4, lwd = 1.5, col = "grey30")  
  lines(Q, Qs.yang, lty = 5, lwd = 1.5, col = "grey40") 
  lines(Qrange, Qs.yang.hg, lty = 6, lwd = 1.5, col = "grey50") 
  lines(Qrange, exp(bhat[1])*Qrange^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",
           "Brownlie / Manning's",
           "Brownlie / Hydr Geom",
           "Yang / Manning's",
           "Yang / Hydr Geom"),
         lwd = 1.5,
         col = c("black", "grey20", "grey30", "grey40", "grey50"),
         lty = c(1,2,4,5,6),
         y.intersp = 0.9,
         bty = "n",
         cex = 0.7)

#   legend("topleft",
#          c("Log-Linear Regression",
#            "Brownlie / Hydr Geom",
#            "Yang / Hydr Geom"),
#          lwd = 1.5,
#          col = c("black", "grey30", "grey50"),
#          lty = c(1,4,6),
#          y.intersp = 0.9,
#          bty = "n",
#          cex = 0.7)

  # Effective Discharge
  plot(q.e,eff.rlm, col = "darkcyan", 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(0, max(eff.br)),
       xlim = c(0,800)
       )
  lines(q.e, eff.br, lty = 2, lwd = 1.5, col = "grey20") 
  lines(q.e, eff.br.hg, lty = 4, lwd = 1.5, col = "grey30") 
  lines(q.e, eff.yang, lty = 5, lwd = 1.5, col = "grey40") 
  lines(q.e, eff.yang.hg, lty = 6, lwd = 1.5, col = "grey50") 
  segments(q.e[which.max(eff.rlm)], 1E-10, q.e[which.max(eff.rlm)], max(eff.rlm), col = "darkcyan")
  segments(q.e[which.max(eff.br)], 1E-10, q.e[which.max(eff.br)], max(eff.br),col = "grey20", lty = 2)
  segments(q.e[which.max(eff.br.hg)], 1E-10, q.e[which.max(eff.br.hg)], max(eff.br.hg), col = "grey30", lty = 4)
  segments(q.e[which.max(eff.yang)], 1E-10, q.e[which.max(eff.yang)], max(eff.yang), col = "grey40", lty = 5)
  segments(q.e[which.max(eff.yang.hg)], 1E-10, q.e[which.max(eff.yang.hg)], max(eff.yang.hg), col = "grey50", 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 = "darkcyan")
  lines(Q.sort, yield.cum.br, lty = 2, lwd = 1.5, col = "grey20")  
  lines(Q.sort, yield.cum.brHG, lty = 4, lwd = 1.5, col = "grey30")  
  lines(Q.sort, yield.cum.yang, lty = 5, lwd = 1.5, col = "grey40") 
  lines(Q.sort, yield.cum.yangHG, lty = 6, lwd = 1.5, col = "grey50") 

  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 = "darkcyan")
   segments(Q.sort[which.min(abs(yield.cum.br-0.5))], -10, Q.sort[which.min(abs(yield.cum.br-0.5))], 0.5,
           col = "grey30", lty = 2)
  segments(Q.sort[which.min(abs(yield.cum.brHG-0.5))], -10, Q.sort[which.min(abs(yield.cum.brHG-0.5))], 0.5,
           col = "grey40", lty = 4)
   segments(Q.sort[which.min(abs(yield.cum.yang-0.5))], -10, Q.sort[which.min(abs(yield.cum.yang-0.5))], 0.5,
            col = "grey50", lty = 5)
  segments(Q.sort[which.min(abs(yield.cum.yangHG-0.5))], -10, Q.sort[which.min(abs(yield.cum.yangHG-0.5))], 0.5,
           col = "grey60", lty = 6)
  segments(-100, 0.5, Q.sort[which.min(abs(yield.cum.rlm-0.5))], 0.5, col = "darkcyan")
  #txt <- sites$site.nm[34]
  text(x=grconvertX(.5, from='npc'), y = grconvertY(.05, from='npc'), labels = txt, cex = 0.75, adj = c(0,0))

dev.off()