yrfun<-function(Fmult,M,sa,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,M,sa,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)
}
project<-function(years,Nhist,Fmort,sa,M,wa,alpha,K,CV.N=0,CV.pred=0,CV.R=0){
  N0<-Nhist*exp(rnorm(1)*CV.N)    # Rather silly - everything has the same CV etc
  Ytraj<-c()
  Straj<-c()
  for(y in years){
    Fimplem<-Fmort*exp(rnorm(1)*CV.pred)  # Implementation error
    Z<-Fimplem*sa+M
    C<-(Fimplem*sa/Z)*(1-exp(-Z))*N0
    S<-sum(wa*pa*N0)
    R<-(alpha*S/(1+S/K))*exp(rnorm(1)*CV.R) # Recr with 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,N=N0))
}