selF1<-0
selF2<-0.15
selF3<-0.35
Fcrash<-0.8
selF4<-1.1

Linf<-160
k<-0.1
beta<-3
cond<-0.02
ages<-1:14
t0<-0
la<-Linf*(1-exp(-k*(ages-t0)))            # Mean length at age
la<-round(la,2)
wa<-cond*la**beta/1000                    # Mean weight at age in kg
wa<-round(wa,2)

s50<-5                                    # Age at 50% selection
sa<-round(1/(1+exp(-1.1*(ages-s50))),2)   # Selection at age
p50<-5.5
pa<-round(1/(1+exp(-2*(ages-p50))),2)   # Proportion mature at age

M<-0.2

yrfun<-function(Fmult,sa,M,wa){
  Fmort<-Fmult*sa
  Z<-Fmort+M
  prop<-(Fmort/Z)*(1-exp(-Z))
  Ztemp<-c(0,Z[1:(length(Z)-1)])
  cumZ<-exp(-cumsum(Ztemp))
  C<-prop*cumZ
  Y<-sum(wa*C)
  return(Y)
}
srfun<-function(Fmult,sa,M,wa,pa){
  Fmort<-Fmult*sa
  Z<-Fmort+M
  Ztemp<-c(0,Z[1:(length(Z)-1)])
  cumZ<-exp(-cumsum(Ztemp))
  S<-sum(wa*pa*cumZ)
  return(S)
}

#alpha<-0.0005
alpha<-1/srfun(Fcrash,sa,M,wa,pa)
K<-20000

Srange<-0:2000*100
Rhat<-alpha*Srange/(1+Srange/K)
plot(Srange,Rhat,type='l',xlab="S ('000 t)",ylab="R (millions)",ylim=c(0,10000))

sr2<-srfun(selF2,sa,M,wa,pa)
sr3<-srfun(selF3,sa,M,wa,pa)
sr4<-srfun(selF4,sa,M,wa,pa)
#sr.crash<-srfun(Fcrash,sa,M,wa,pa)
#lines(Srange,Srange/sr.crash)
text(150000,6500,"Beverton-Holt curve")
text(150000,6000,expression(R=alpha*S/(1+S/K)))
text(35000,7000,"Replacement")
#text(15000,4500,"Slope for\n Fcrash is")
#text(15000,6000,expression(alpha))

S0<-100000
R0<-alpha*S0/(1+S0/K)
S1<-R0*sr4
R1<-alpha*S1/(1+S1/K)
S2<-R1*sr4
R2<-alpha*S2/(1+S2/K)
S3<-R2*sr4
R3<-alpha*S3/(1+S3/K)
S4<-R3*sr4
R4<-alpha*S4/(1+S4/K)
lines(c(S0,S0,S1,S1,S2,S2,S3,S3,S4),c(0,R0,R0,R1,R1,R2,R2,R3,R3))
lines(Srange,Srange/sr4)