#Set up the working environment

rm(list = ls())

library(reshape2)
library(xtable)
library(ggplot2)

# RUN IF YOU HAVE NOT INSTALLED THE BICALC PACKAGE
#  library(devtools)
#  install_github("bceaton/bicalc")

library(bicalc)

load("GSD_data.RData")  #load the dataset
load("GSD_simulations.RData")  #load the first dataset in the list

#Analysis in the paper The following code chunks contain the commands used to create all the figures in the paper, using the data files provided.

#sample.data = sample(all1.data$D, size = k)
d.demo = sort(sample(demo.data$D, size = 200))
cfd.demo = seq_along(d.demo)/length(d.demo)

# calculate the binomial limits for specified values n and p
n = 200
p = 0.84
firstINT = QuantBD(n, p)

#the exact binomial confidence interval, based on measurements
l = firstINT$interval[1]
u = firstINT$interval[2]
dl = d.demo[l]
du = d.demo[u]

#the equal tail approximation applied to the complete dataset
le = firstINT$equaltail[1]
ue = firstINT$equaltail[2]
dle = as.numeric(approx(x=seq_along(d.demo), y = d.demo, xout = le )[2])
due = as.numeric(approx(x=seq_along(d.demo), y = d.demo, xout = ue )[2])

#the equal tail approximation applied to the binned data
d.binned = MakeCFD(d.demo)
dlb = 2^as.numeric(approx(x=d.binned$probs, y = log2(d.binned$size), xout = le/200 )[2])
dub = 2^as.numeric(approx(x=d.binned$probs, y = log2(d.binned$size), xout = ue/200 )[2])
dest = 2^as.numeric(approx(x=d.binned$probs, y = log2(d.binned$size), xout = 0.84 )[2])
plot(d.demo, cfd.demo,
     type ="b",
     pch = 21,
     #cex = 0.75,
     log = "x",
     col = rgb(0,0,1,0.35),
     ylim = c(0, 1),
     xlim = c(0.5, 8),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
lines(MakeCFD(d.demo),
      col = rgb(1,0,0,0.75),
      lwd = 1.5,
      lty =1)

lines(c(0.4, dub), c(ue/200, ue/200), lty = 2)
lines( c(dub, dub), c(-0.1, ue/200), lty = 2)
lines(c(0.4, dlb), c(le/200, le/200), lty = 2)
lines( c(dlb, dlb), c(-0.1, le/200), lty = 2)

lines(c(0.4, dest), c(0.84, 0.84), lty = 1)
lines( c(dest, dest), c(-0.1, 0.84), lty = 1)

k = 150:190
pr = dbinom(k, n, p)
df = data.frame(k, pr)

# find upper and lower limits and output coverage probability.
qci = QuantBD(n, p)
lu_exact = qci$interval
p_c = qci$coverage


df$bar_col = ifelse(df$k >= lu_exact[1] & df$k < lu_exact[2], "in", "out")

# using ggplot since barplots in base graphics are a bit awkward
par(mar=c(4,4,1,1))
par(oma = c(0,0,0,0))

ggplot(df, aes(x = k, y = pr)) +
  geom_col(aes(fill = bar_col), na.rm = T) +
  scale_x_continuous(limits = c(150, 190)) +
  scale_y_continuous(limits = c(0, 0.08)) +
  scale_fill_manual(values = c("darkgrey", "lightgrey"), name = "") +
  labs(x = expression(italic(k)),
       y = expression(italic(P[r](k, 200, 0.84)))) + 
  theme_bw() + 
  theme(axis.text = element_text(size = 12)) +
  theme(axis.title = element_text(size = 12)) +
  geom_vline(xintercept = firstINT$equaltail+0.5, lty = 2)

par(mfrow = c(1,2))

mycolors = c(rgb(0,0,1,0.4), rgb(0,1,1,0.4), rgb(1,1,0,0.4), rgb(1,0,0,0.4))
linecolors = c(rgb(0,0,1,0.6), rgb(0,1,1,1), rgb(1,1,0,1), rgb(1,0,0,0.6))

plot(cfd1hr$size, cfd1hr$probs,
     type ="l",
     lty = 0,
     log = "x",
     col = linecolors[1],
     ylim = c(0, 1),
     xlim = c(0.5, 5),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
  polygon(polyhr1,
          col=mycolors[1],
          lty = 0)
  polygon(polyhr2,
          col=mycolors[2],
          lty = 0)
  polygon(polyhr3,
          col=mycolors[3],
          lty = 0)
  polygon(polyhr4,
          col=mycolors[4],
          lty = 0)
   lines(cfd1hr$size, cfd1hr$probs,
        col = linecolors[1],
        lty = 2)
  lines(cfd2hr$size, cfd2hr$probs,
        col = linecolors[2],
        lty = 2)
  lines(cfd3hr$size, cfd3hr$probs,
        col = linecolors[3],
        lty = 2)
  lines(cfd4hr$size, cfd4hr$probs,
        col = linecolors[4],
        lty = 2)
  #abline(h = c(0.05, 0.95), lty = 2)
  legend("topleft",
       inset = 0.02,
       cex = 0.75,
       bg = rgb(1,1,1),
       #box.lty = 0,
       legend = c("sample 1", "sample 2", "sample 3", "sample 4"),
       lty = c(0,0,0,0,2),
       pch = c(15,15, 15, 15, NA),
       pt.cex = c(2,2,2,2,1),
       col = c(mycolors, 
               rgb(0,0,0))
  )

  
plot(cfd5hr$size, cfd5hr$probs,
     type ="l",
     lty = 0,
     log = "x",
     col = rgb(0,0,1),
     ylim = c(0, 1),
     xlim = c(0.5, 5),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
  polygon(polyhr5,
          col=mycolors[1],
          lty = 0)
  polygon(polyhr6,
          col=mycolors[2],
          lty = 0)
  polygon(polyhr7,
          col=mycolors[3],
          lty = 0)
  polygon(polyhr8,
          col=mycolors[4],
          lty = 0)
  lines(cfd5hr$size, cfd5hr$probs,
        col =linecolors[1],
        lty = 2)
  lines(cfd6hr$size, cfd6hr$probs,
        col =linecolors[2],
        lty = 2)
  lines(cfd7hr$size, cfd7hr$probs,
        col = linecolors[3],
        lty = 2)
  lines(cfd8hr$size, cfd8hr$probs,
        col = linecolors[4],
        lty = 2)
  #abline(h = c(0.05, 0.95), lty = 2)
  legend("topleft",
       inset = 0.02,
       cex = 0.75,
       bg = rgb(1,1,1),
       #box.lty = 0,
       legend = c("sample 5", "sample 6", "sample 7", "sample 8"),
       lty = c(0,0,0,0,2),
       pch = c(15,15, 15, 15, NA),
       pt.cex = c(2,2,2,2,1),
       col = mycolors
  )

par(mfrow = c(1,3))

ind = 2
boxplot(results[,ind],
        ylim = c(ci$lower[ind], ci$upper[ind])*c(0.8, 1.2),
        xlab = "")
title(xlab="D16", line=1, cex.lab=1.2)
abline(h = ci$estimate[ind], lty = 2, col = rgb(0,0,0))
abline(h = c(ci$upper[ind], ci$lower[ind]), lty = 2, col = rgb(1,0,0))
abline(h = c(ci50$upper[ind], ci50$lower[ind]), lty = 2, col = rgb(0,0,1))
d16.in95ci =sum(results[,ind]>ci$lower[ind] & results[,ind]<ci$upper[ind])/length(results[,ind])
d16.in50ci =sum(results[,ind]>ci50$lower[ind] & results[,ind]<ci50$upper[ind])/length(results[,ind])

ind = 4
boxplot(results[,ind],
        ylim = c(ci$lower[ind], ci$upper[ind])*c(0.8, 1.2),
        xlab = "")
title(xlab="D50", line=1, cex.lab=1.2)
abline(h = ci$estimate[ind], lty = 2, col = rgb(0,0,0))
abline(h = c(ci$upper[ind], ci$lower[ind]), lty = 2, col = rgb(1,0,0))
abline(h = c(ci50$upper[ind], ci50$lower[ind]), lty = 2, col = rgb(0,0,1))
d50.in95ci =sum(results[,ind]>ci$lower[ind] & results[,ind]<ci$upper[ind])/length(results[,ind])
d50.in50ci =sum(results[,ind]>ci50$lower[ind] & results[,ind]<ci50$upper[ind])/length(results[,ind])

ind = 6
boxplot(results[,ind],
        ylim = c(ci$lower[ind], ci$upper[ind])*c(0.8, 1.2),
        xlab = "")
title(xlab="D84", line=1, cex.lab=1.2)
abline(h = ci$estimate[ind], lty = 2, col = rgb(0,0,0))
abline(h = c(ci$upper[ind], ci$lower[ind]), lty = 2, col = rgb(1,0,0))
abline(h = c(ci50$upper[ind], ci50$lower[ind]), lty = 2, col = rgb(0,0,1))

d84.in95ci =sum(results[,ind]>ci$lower[ind] & results[,ind]<ci$upper[ind])/length(results[,ind])
d84.in50ci =sum(results[,ind]>ci50$lower[ind] & results[,ind]<ci50$upper[ind])/length(results[,ind])
par(mfcol=c(1,2))

plot(k.surf$size, k.surf$probs,
     type ="b",
     log = "x",
     pch = 19,
     col = "blue",
     ylim = c(0, 1),
     xlim = c(2,200),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
text(x = 180, y = 0.01, labels = "A")
lines(k.bulk$size, k.bulk$probs, lty = 2, col = rgb(1,0,0))
points(k.bulk$size, k.bulk$probs, pch = 1, col = rgb(1,0,0))


plot(k.surf$size, k.surf$probs,
     type ="p",
     pch = 19,
     log = "x",
     col = "blue",
     ylim = c(0, 1),
     xlim = c(2,200),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )

polygon(poly.surf,
        col=rgb(0, 0, 1,0.3),
        lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
abline(v = c(22.6), lty = 2, col = rgb(1,0,0))

#lines(k.bulk$size, k.bulk$probs, lty = 2, col = rgb(1,0,0))
points(k.bulk$size, k.bulk$probs, pch = 1, col = rgb(1,0,0))

text(x = 180, y = 0.01, labels = "B")

legend("topleft",
       inset = 0.02,
       cex = 0.75,
       bg = rgb(1,1,1),
       #box.lty = 0,
       legend = c("pebble count", "bulk sample (> 4 mm)", "22.6 mm"),
       lty = c(0,0,2),
       pch = c(19,1, NA),
       col = c(rgb(0, 0, 1), rgb(1,0,0), rgb(1,0,0))
)

par(mfcol=c(1,2))
 
plot(cfdpool$size, cfdpool$probs,
     type ="p",
     pch = 19,
     log = "x",
     col = "blue",
     ylim = c(0, 1),
     xlim = c(2,400),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
polygon(polypool,
        col=rgb(0, 0, 1,0.3),
        lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
abline(v = c(22.6), lty = 2, col = rgb(1,0,0))
text(x = 360, y = 0.01, labels = "A")

#lines(cfdbar$size, cfdbar$probs, lty = 1, col = rgb(0.5,0,1))
points(cfdbar$size, cfdbar$probs, pch = 19, col = rgb(0.5,0,1))
polygon(polybar,
        col=rgb(0.5, 0, 1,0.3),
        lty = 0)

#lines(cfdrun$size, cfdrun$probs, lty = 1, col = rgb(1,0.25,0))
points(cfdrun$size, cfdrun$probs, pch = 19, col = rgb(1,0.25,0))
polygon(polyrun,
        col=rgb(1, 0.25, 0,0.3),
        lty = 0)

plot(cfdpool4$size, cfdpool4$probs,
     type ="p",
     pch = 19,
     log = "x",
     col = "blue",
     ylim = c(0, 1),
     xlim = c(2,400),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
polygon(polypool4,
        col=rgb(0, 0, 1,0.3),
        lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
abline(v = c(22.6), lty = 2, col = rgb(1,0,0))

text(x = 360, y = 0.01, labels = "B")

#lines(cfdbar4$size, cfdbar4$probs, lty = 1, col = rgb(0.5,0,1))
points(cfdbar4$size, cfdbar4$probs, pch = 19, col = rgb(0.5,0,1))
polygon(polybar4,
        col=rgb(0.5, 0, 1,0.3),
        lty = 0)

#lines(cfdrun4$size, cfdrun4$probs, lty = 1, col = rgb(1,0.25,0))
points(cfdrun4$size, cfdrun4$probs, pch = 19, col = rgb(1,0.25,0))
polygon(polyrun4,
        col=rgb(1, 0.25, 0,0.3),
        lty = 0)
legend("topleft",
       inset = 0.02,
       cex = 0.75,
       bg = rgb(1,1,1),
       #box.lty = 0,
       legend = c("pool", "exposed bar", "run/riffle", "22.6 mm"),
       lty = c(0,0,0,2),
       pch = c(19,19, 19, NA),
       col = c(rgb(0, 0, 1), rgb(0.5, 0, 1), rgb(1, 0.25, 0), rgb(1,0,0))
)

par(mfcol = c(1,2))

plot(cfdopA$size, cfdopA$probs,
     type ="p",
     pch = 1,
     log = "x",
     col = "blue",
     ylim = c(0, 1),
     xlim = c(2,400),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
lines(cfdopA$size, cfdopA$probs, lty = 1, col = rgb(0,0,1))

lines(cfdopB$size, cfdopB$probs, lty = 1, col = rgb(1,0,0))
points(cfdopB$size, cfdopB$probs, pch = 19, col = rgb(1,0,0))
text(x = 360, y = 0.01, labels = "A")

plot(cfdopA$size, cfdopA$probs,
     type ="p",
     pch = 1,
     log = "x",
     col = "blue",
     ylim = c(0, 1),
     xlim = c(2,400),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
polygon(polyopA,
        col=rgb(0, 0, 1,0.3),
        lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
#abline(v = c(16), lty = 1, col = rgb(1,0,0))

#lines(cfdbar$size, cfdbar$probs, lty = 1, col = rgb(1,0,0))
points(cfdopB$size, cfdopB$probs, pch = 19, col = rgb(1,0,0))
polygon(polyopB,
        col=rgb(1, 0, 0,0.3),
        lty = 0)
text(x = 360, y = 0.01, labels = "B")

legend("topleft",
       inset = 0.02,
       cex = 0.75,
       bg = rgb(1,1,1),
       #box.lty = 0,
       legend = c("Operator A", "Operator B"),
       lty = c(0,0),
       pch = c(1,19),
       col = c(rgb(0, 0, 1),rgb(1,0,0))
)

par(mfcol = c(1,2))

plot(cfdheel$size, cfdheel$probs,
     type ="p",
     pch = 1,
     log = "x",
     col = "blue",
     ylim = c(0, 1),
     xlim = c(2,400),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
lines(cfdheel$size, cfdheel$probs, lty = 1, col = rgb(0,0,1))

lines(cfdframe$size, cfdframe$probs, lty = 1, col = rgb(1,0,0))
points(cfdframe$size, cfdframe$probs, pch = 19, col = rgb(1,0,0))
text(x = 360, y = 0.01, labels = "A")

plot(cfdheel$size, cfdheel$probs,
     type ="p",
     pch = 1,
     log = "x",
     col = "blue",
     ylim = c(0, 1),
     xlim = c(2,400),
     xlab = "grain size (mm)",
     ylab = "Proportion Finer"
     )
polygon(polyheel,
        col=rgb(0, 0, 1,0.3),
        lty = 0)
#abline(h = c(0.05, 0.95), lty = 2)
#abline(v = c(16), lty = 1, col = rgb(1,0,0))

#lines(cfdbar$size, cfdbar$probs, lty = 1, col = rgb(1,0,0))
points(cfdframe$size, cfdframe$probs, pch = 19, col = rgb(1,0,0))
polygon(polyframe,
        col=rgb(1, 0, 0,0.3),
        lty = 0)
text(x = 360, y = 0.01, labels = "B")

legend("topleft",
       inset = 0.02,
       cex = 0.75,
       bg = rgb(1,1,1),
       #box.lty = 0,
       legend = c("Heel-toe", "Sampling frame"),
       lty = c(0,0),
       pch = c(1,19),
       col = c(rgb(0, 0, 1),rgb(1,0,0))
)

CalcError = function(p, X, n= c(100,200,300,400,500,600,700, 800, 900, 1000)){
  tmp = rbind(WolmanCI(X,n[1],probs = p),
                    WolmanCI(X,n[2], probs =p), 
                    WolmanCI(X,n[3], probs =p),
                    WolmanCI(X,n[4], probs =p),
                    WolmanCI(X,n[5], probs =p),
                    WolmanCI(X,n[6], probs =p),
                    WolmanCI(X,n[7], probs =p),
                    WolmanCI(X,n[8], probs =p),
                    WolmanCI(X,n[9], probs =p),
                    WolmanCI(X,n[10], probs =p)
                    )
err = as.numeric(0.5*(tmp$upper - tmp$lower)/tmp$estimate)
return(err)
}

CalcSDlog = function(X){
  tmp = WolmanCI(X, 100, p = c(0.5, 0.84))
  SDlog = diff(log2(tmp$estimate))
  return(SDlog)
}

srcs = c("C1", "C2", "C3", "C4", "C5", "C6", "C7", "C8", "cfdopA", "cfdopB", "cfdpool", "cfdpool4", "cfdrun", "cfdrun4", "cfdbar", "cfdbar4", "k.surf", "cfd1hr", "cfd2hr", "cfd3hr", "cfd4hr", "cfd5hr", "cfd6hr", "cfd6hr", "cfd7hr", "cfd8hr", "G309A", "G309B")

Err50 = matrix(data=NA, nrow = 10 , ncol =  length(srcs)+1)
Err84 = matrix(data=NA, nrow = 10, ncol =  length(srcs)+1)
SDlog = numeric(length = length(srcs))

Err50[,1] = seq(100,1000,100)
Err84[,1] = Err50[,1]

for(i in seq_along(srcs)) {
  Err50[,i+1] = CalcError(0.5, get(srcs[i]))
  Err84[,i+1] = CalcError(0.84, get(srcs[i]))
  SDlog[i] = CalcSDlog(get(srcs[i]))
}

n100 = data.frame(med50 = NA)
n100$med50 = mean(( Err50[1,2:19]))
n100$min50 = min(( Err50[1,2:19]))
n100$max50 = max(( Err50[1,2:19]))
n100$med84 = mean(( Err84[1,2:19]))
n100$min84 = min(( Err84[1,2:19]))
n100$max84 = max(( Err84[1,2:19]))


n200 = data.frame(med50 = NA)
n200$med50 = mean(( Err50[2,2:19]))
n200$min50 = min(( Err50[2,2:19]))
n200$max50 = max(( Err50[2,2:19]))
n200$med84 = mean(( Err84[2,2:19]))
n200$min84 = min(( Err84[2,2:19]))
n200$max84 = max(( Err84[2,2:19]))

n500 = data.frame(med50 = NA)
n500$med50 = mean(( Err50[5,2:19]))
n500$min50 = min(( Err50[5,2:19]))
n500$max50 = max(( Err50[5,2:19]))
n500$med84 = mean(( Err84[5,2:19]))
n500$min84 = min(( Err84[5,2:19]))
n500$max84 = max(( Err84[5,2:19]))
par(mfcol = c(1,2))
par(mar=c(4,4,1,1))
par(oma = c(0,0,0,0))

matplot(Err50[,1], Err50[,2:19], 
        type = "l", 
        col = rgb(0,0,1, 0.3), 
        lty = 1, 
        lwd = 2,
        ylim = c(0, 0.6),
        xlim = c(100,1000),
        xlab = "sample size",
        ylab = expression(epsilon[50])
)
abline(v = c(200,500), lty = 1, lwd = 4, col = rgb(0,0,0,0.3))
text(x = 950, y = 0.0, labels = "A")

matplot(Err84[,1], Err84[,2:19], 
        type = "l", 
        col = rgb(0,0,1, 0.3), 
        lty = 1, 
        lwd = 2,
        ylim = c(0, 0.6),
        xlim = c(100,1000),
        xlab = "sample size",
        ylab = expression(epsilon[84])
)
abline(v = c(200,500), lty = 1, lwd = 4, col = rgb(0,0,0,0.3))
text(x = 950, y = 0.0, labels = "B")

par(mfcol = c(1,2))

plot.cols = character(length = NSim)
mycolors = colorRampPalette(c("blue", "cyan", "yellow", "red"))(6)
filt = SErr50$logSD < 0.75
plot.cols[filt] = mycolors[1]
filt = SErr50$logSD >= 0.75 & SErr50$logSD < 1.0
plot.cols[filt] = mycolors[2]
filt = SErr50$logSD >= 1 & SErr50$logSD < 1.25
plot.cols[filt] = mycolors[3]
filt = SErr50$logSD >= 1.25 & SErr50$logSD < 1.5
plot.cols[filt] = mycolors[4]
filt = SErr50$logSD >= 1.5 & SErr50$logSD < 1.75
plot.cols[filt] = mycolors[5]
filt = SErr50$logSD >= 1.75 & SErr50$logSD < 2.0
plot.cols[filt] = mycolors[6]

plot(SErr50$N, SErr50$err,
     type = "p",
     pch = 20,
     cex = 0.75,
     col = plot.cols,
     ylim = c(0, 0.6), 
     xlim=(c(0,1000)),
     xlab = "sample size",
      ylab = expression(epsilon[50])
     )

X = seq(0,1000,1)
Y1 = exp(powfit50$coefficients[1]+powfit50$coefficients[3]*0.5)*X^powfit50$coefficients[2]
Y2 = exp(powfit50$coefficients[1]+powfit50$coefficients[3]*2.0)*X^powfit50$coefficients[2]
lines(X,Y2, lwd = 2, col = "black" )
lines(X,Y1, lwd =2, col = "black", lty = 2 )
text(x = 950, y = 0.0, labels = "A")

legend("topright",
       inset = 0.02,
       cex = 1.0,
       bg = rgb(1,1,1),
       #box.lty = 0,
       legend = c(expression(2.0~phi), 
                  expression(0.5~phi),
                  expression((0.5-0.75~phi)),
                  expression((0.75-1.0~phi)),
                  expression((1.0-1.25~phi)),
                  expression((1.25-1.5~phi)),
                  expression((1.5-1.75~phi)),
                  expression((1.75-2.0~phi))
                  ),
       lty = c(1,2,0,0,0,0,0,0),
       lwd = c(2,2,0,0,0,0,0,0),
       pch = c(NA,NA,20,20,20,20,20,20),
       col = c("black", "black", mycolors)
)

plot(SErr84$N, SErr84$err,
     type = "p",
     pch = 20,
     cex = 0.75,
     col = plot.cols,
     ylim = c(0, 0.6), 
     xlim=(c(0,1000)),
     xlab = "sample size",
      ylab = expression(epsilon[84])
     )

X = seq(0,1000,1)
Y3 = exp(powfit84$coefficients[1]+powfit84$coefficients[3]*0.5)*X^powfit84$coefficients[2]
Y4 = exp(powfit84$coefficients[1]+powfit84$coefficients[3]*2.0)*X^powfit84$coefficients[2]
lines(X,Y4, lwd = 2, col = "black" )
lines(X,Y3, lwd = 2, col = "black", lty = 2 )
text(x = 950, y = 0.0, labels = "B")

Coef = c("A","B")
CoefTable = data.frame(Coef)
CoefTable$SD0.50 = c(coeff50[1] + coeff50[3]*0.5,
                     coeff84[1] + coeff84[3]*0.5)
CoefTable$SD0.75 = c(coeff50[1] + coeff50[3]*0.75,
                    coeff84[1] + coeff84[3]*0.75)
CoefTable$SD1.00 = c(coeff50[1] + coeff50[3]*1.0,
                    coeff84[1] + coeff84[3]*1.0)
CoefTable$SD1.25 = c(coeff50[1] + coeff50[3]*1.25,
                  coeff84[1] + coeff84[3]*1.25)
CoefTable$SD1.50 = c(coeff50[1] + coeff50[3]*1.5,
                  coeff84[1] + coeff84[3]*1.5)
CoefTable$SD1.75 = c(coeff50[1] + coeff50[3]*1.75,
                    coeff84[1] + coeff84[3]*1.75)
CoefTable$SD2.00 = c(coeff50[1] + coeff50[3]*2.0,
                  coeff84[1] + coeff84[3]*2.00)
#<<fishtrap, echo=FALSE>>= 
D50 = 0.055
Tstar = 0.045
S = 0.02
Tcrit = 9.81*1650*Tstar*D50
dcrit = Tcrit/(9810*S)
epsilon100 = exp(coeff50[1] + coeff50[3]*1.0)*100^coeff50[2]
D50low = 0.055*(1-epsilon100)
D50up = 0.055*(1+epsilon100)
Tcritlow = 9.81*1650*Tstar*D50low
dcritlow = Tcritlow/(9810*S)
Tcritup = 9.81*1650*Tstar*D50up
dcritup = Tcritup/(9810*S)

Q = 2.5
h = 0.23*Q^0.4
offset = h - dcrit
hfull = offset + 2*dcrit
Qfull = (hfull/0.23)^(1/0.4)
hlow = offset + 2*dcritlow
Qlow = (hlow/0.23)^(1/0.4)
hup = offset + 2*dcritup
Qup = (hup/0.23)^(1/0.4)


# Specify Fishtrap Creek peak flows for 1972 - 2011
Q = c(14, 3.77, 8.72, 12, 4.11, 6.03, 4.33, 3.31, 3.78, 9.64,
9.55, 9.27, 4.13, 11.7, 3.6, 5.85, 3.98, 3.88, 10.5, 6.01, 9.67,
10.8, 15.3, 10.1, 13.5, 5.05, 3.13, 8.68, 4.2, 6.3, 10.1, 8.8, 8.02, 
8.08, 8.02, 7.9, 9.73)

# Set up x axis tick positions and labels
Ttick = c(1.001,1.01,1.1,1.5,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,25,30,35,40,45,50,60,70,80,90,100)

# Fit a line by method of moments, along with 95% confidence intervals
KTtick = -(sqrt(6)/pi)*(0.5772 + log(log(Ttick/(Ttick-1))))
QTtick = mean(Q) + KTtick*sd(Q) 

RPlow = as.numeric(approx(x = QTtick, y = Ttick, xout = Qlow)[2])
RPup = as.numeric(approx(x = QTtick, y = Ttick, xout = Qup)[2])

epsilon500 = exp(coeff50[1] + coeff50[3]*1.0)*500^coeff50[2]
D50low2 = 0.055*(1-epsilon500)
D50up2 = 0.055*(1+epsilon500)
Tcritlow2 = 9.81*1650*Tstar*D50low2
dcritlow2 = Tcritlow2/(9810*S)
Tcritup2 = 9.81*1650*Tstar*D50up2
dcritup2 = Tcritup2/(9810*S)

hlow2 = offset + 2*dcritlow2
Qlow2 = (hlow2/0.23)^(1/0.4)
hup2 = offset + 2*dcritup2
Qup2 = (hup2/0.23)^(1/0.4)

RPlow2 = as.numeric(approx(x = QTtick, y = Ttick, xout = Qlow2)[2])
RPup2 = as.numeric(approx(x = QTtick, y = Ttick, xout = Qup2)[2])