selF1<-0
selF2<-0.25
selF3<-0.35
selF4<-1.1                                # Chosen as Fcrash

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
recrCV<-0.4

yrfun<-function(Fmult){
  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){
  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)
}
years<-1:5
project<-function(Nhist,Fmort,sa,M,wa,alpha,K){
  N0<-Nhist
  Z<-Fmort*sa+M
  Ytraj<-c()
  Straj<-c()
  for(y in years){
    C<-(Fmort*sa/Z)*(1-exp(-Z))*N0
    S<-sum(wa*pa*N0)
    #R<-(alpha*S/(1+S/K))*exp(rnorm(1)*recrCV) # Recr with variation
    R<-(alpha*S/(1+S/K))                  # Recr without variation
    N1<-N0*exp(-Z)
    N0<-c(R,N1[1:(length(N1)-1)])
    Ytraj<-c(Ytraj,sum(wa*C))
    Straj<-c(Straj,S)
  }
  return(list(Y=Ytraj,S=Straj))
}
alpha<-1/srfun(selF4)               # Make this Fcrash
K<-20000

Frange<-(0:(selF4*100))/100         # Range for plotting
sr<-sapply(Frange,srfun)
yr<-sapply(Frange,yrfun)
Srange<-K*(alpha*sr-1)
Srange<-ifelse(Srange>0,Srange,0)
Rrange<-alpha*Srange/(1+Srange/K)
Yrange<-yr*Rrange
#par(mfrow=c(2,2))
#plot(Frange,yr,type='l',lwd=3)
#plot(Frange,sr,type='l',lwd=3)
#plot(Srange,Rrange,type='l',lwd=3)

Fhist<-selF3
Zhist<-Fhist*sa+M
srhist<-srfun(Fhist)
Shist<-K*(alpha*srhist-1)
Rhist<-alpha*Shist/(1+Shist/K)
Nhist<-Rhist*exp(cumsum(-c(0,Zhist[1:(length(Zhist)-1)])))
#barplot(Nhist,names.arg=ages)
#plot(Srange,Yrange,type='l',lwd=3)
#lines(Srange,0.10*Srange,lwd=3,col="blue")