selF1<-0
selF2<-0.25
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:1000*100
Rhat<-alpha*Srange/(1+Srange/K)
par(mfcol=c(1,2))
plot(Srange,Rhat,type='l',xlab="S ('000 t)",ylab="R (millions)",ylim=c(0,15000))

sr.crash<-srfun(Fcrash,sa,M,wa,pa)
lines(Srange,Srange/sr.crash)
text(150000,8500,"Beverton-Holt curve")
text(150000,8000,expression(R=alpha*S/(1+S/K)))
text(15000,6500,"Slope for\n Fcrash is")
text(15500,5500,expression(alpha))


F.traj<-seq(0,1.5,0.01)
sr.traj<-sapply(F.traj,srfun,sa,M,wa,pa)
plot(F.traj,sr.traj,type='l',ylab="S/R",xlab="F")
lines(c(Fcrash,Fcrash,0),c(0,sr.crash,sr.crash))
text(.1,sr.crash+0.75,expression(1/alpha))
text(Fcrash+.2,1,"Fcrash")