# simpop.r
# statistical stock assessment
# "few" true ages, a=1,...,A and a plus group
#  Data:  indices for each age group
#         and catch (yield)
#         w: vector of length A+1 - weight at age 
#  Parameters:
#   Recr, Fmort,
#   selpat: vector of length A+1 - selection pattern
#   Ninit : vector of length A+1 - initial pop size
#   M,q
# 
#  Dimensions:  A=#true ages in data
#               numyears= # years
#               Apop=#ages in pop

R0<-100
M<-0.7
numyears<-10
CVrecr<-0.3
Apop<-3
Recr<-R0*exp(rnorm(numyears)*CVrecr)
Ninit<-Recr[1]*exp(-(0:Apop)*M)
Nmat<-Ninit
N0<-Ninit
Fmort<-rep(0.2,numyears)
A<-3
q<-c(0.1,rep(0.5,Apop))
qoldest<-q[length(q)]
w<-c(0.03,0.15,0.20,rep(0.25,Apop-A+1))
selpat<-c(0,rep(1,Apop))
Y<-c()
for(y in 1:(numyears-1)){
  Z<-Fmort[y]*selpat+M
  C<-((Fmort[y]*selpat)/Z)*(1-exp(-Z))*N0
  Y<-c(Y,sum(w*C))
  N1<-c(Recr[y+1],N0[1:(Apop-1)]*exp(-Fmort[y]*selpat[1:(Apop-1)]-M),
        N0[Apop]*exp(-Fmort[y]*selpat[Apop]-M)+
        N0[Apop+1]*exp(-Fmort[y]*selpat[Apop+1]-M))
  Nmat<-rbind(Nmat,N1)
  N0<-N1

}
Z<-Fmort[numyears]*selpat+M
C<-((Fmort[numyears]*selpat)/Z)*(1-exp(-Z))*N0
Y<-c(Y,sum(w*C))

dimnames(Nmat)<-list(Years=1:(numyears),Ages=c(1:Apop,"+"))
I<-t(t(Nmat)*q)
U<-qoldest*c(Nmat%*%w)
I<-cbind(I[1:numyears,1:3],U)
params.true<-log(c(Fmort,Recr,q))     # True values of all parameters