#Build Commands Fig 3.1
 curve_number<-seq(1,100,1)
 Ia_black<-1.299918-0.217235*0+0.211964*0-0.007108*curve_number
 Ia_blue<-1.299918-0.217235*1+0.211964*1-0.007108*curve_number
 Ia_green<-1.299918-0.217235*0+0.211964*1-0.007108*curve_number
 Ia_red<-1.299918-0.217235*1+0.211964*0-0.007108*curve_number
 #
MCL_S<-seq(1,20000,5)
Cl_black=0.8406-0.2035*0+0.3349*0-0.00005*MCL_S
Cl_blue=0.8406-0.2035*1+0.3349*1-0.00005*MCL_S
Cl_green=0.8406-0.2035*0+0.3349*1-0.00005*MCL_S
Cl_red=0.8406-0.2035*1+0.3349*0-0.00005*MCL_S
 #
leinhard<-function(t,trms,n,beta){
term1<-beta/(gamma(n/beta))
term2<-(n/beta)^(n/beta)
term3<-(1/trms)
term4<-(t/trms)^(n-1)
term5<-((n/beta)*(t/trms)^beta)
term1*term2*term3*term4*exp(-term5)}
#
trms_black<-10^(1.25258-0.36105*0+0.10253*0)*(MCL_S)^0.40163
trms_blue<-10^(1.25258-0.36105*1+0.10253*1)*(MCL_S)^0.40163
trms_green<-10^(1.25258-0.36105*0+0.10253*1)*(MCL_S)^0.40163
trms_red<-10^(1.25258-0.36105*1+0.10253*0)*(MCL_S)^0.40163
Tp_black<-trms_black/((3.262480/(3.262480-1))^0.5)
Tp_blue<-trms_blue/((3.262480/(3.262480-1))^0.5)
Tp_green<-trms_green/((3.262480/(3.262480-1))^0.5)
Tp_red<-trms_red/((3.262480/(3.262480-1))^0.5)
Qp_black<-leinhard(Tp_black,trms_black,3.262480,2)
Qp_blue<-leinhard(Tp_blue,trms_blue,3.262480,2)
Qp_green<-leinhard(Tp_green,trms_green,3.262480,2)
Qp_red<-leinhard(Tp_red,trms_red,3.262480,2)
in2cfs<-(640*43560)/(60*12)
#
plot(MCL_S,Qp_black*in2cfs,log="xy",tck=1,type="l",lwd=4,xlim=c(1,20000),ylim=c(20,5000),xlab="MCL/S in miles",ylab="Qp in cfs-hr/in-sq.mi.")
lines(MCL_S,Qp_blue*in2cfs,col=rgb(0,0,1,1),lwd=4)
lines(MCL_S,Qp_green*in2cfs,col=rgb(0,1,0,1),lwd=4)
lines(MCL_S,Qp_red*in2cfs,col=rgb(1,0,0,1),lwd=4)
#
plot(MCL_S,Tp_black/60,log="xy",type="l",lwd=4,xlim=c(1,20000),ylim=c(0.2,20),tck=1,xlab="MCL/S in miles",ylab="Tp in hours")
lines(MCL_S,Tp_blue/60,col=rgb(0,0,1,1),lwd=4)
lines(MCL_S,Tp_green/60,col=rgb(0,1,0,1),lwd=4)
lines(MCL_S,Tp_red/60,col=rgb(1,0,0,1),lwd=4)
#
plot(MCL_S,Cl_black,tck=1,type="l",lwd=4,xlim=c(1,20000),ylim=c(0,1.2),xlab="MCL/S in miles",ylab="Cl in inches per hour")
lines(MCL_S,Cl_blue,lwd=4,col=rgb(0,0,1,1))
lines(MCL_S,Cl_green,lwd=4,col=rgb(0,1,0,1))
lines(MCL_S,Cl_red,lwd=4,col=rgb(1,0,0,1))
#
plot(curve_number,Ia_black,lwd=4,xlim=c(50,100),ylim=c(0.4,1.1),xlab="CN (curve number)",ylab="Ia in inches",type="l",tck=1)
lines(curve_number,Ia_blue,lwd=4,col=rgb(0,0,1,1))
lines(curve_number,Ia_green,lwd=4,col=rgb(0,1,0,1))
lines(curve_number,Ia_red,lwd=4,col=rgb(1,0,0,1))