MaxAge <- 50
N <- matrix(0,nrow=MaxAge,ncol=3)

# Set M and agerec (note that first age-class is the age-at-recruitmetn - hence the -1)
M <-0.18
AgeRec <- 5-1

S <- rep(exp(-M),3)
Sf <- rep(exp(-M)^0.694,3)
print(S)
print(Sf)

X <- matrix(c(0.1761,0.0000,0.0000,0.7052,0.2206,0.0000,0.1187,0.7794,1.0000),ncol=3,byrow=T)
print(X)

# classes 2 and 3 are mature
fec <- c(0,0.001165731,0.001930932)

N[1,] <- c(1,0,0)

for (Iage in 2:MaxAge)
{
 # growth at the of the year after mortality 
 N[Iage,] <- X%*%(S*N[Iage-1,])
}  

Nsum1 <- 0; Nsum2 <- 0
FecAge <- rep(0,MaxAge)
for (Iage in 1:MaxAge)
 { 
  Nsum1 <- Nsum1 + Iage*sum(N[Iage,]*Sf*fec)
  Nsum2 <- Nsum2 + sum(N[Iage,]*Sf*fec)
  FecAge[Iage] <- sum(N[Iage,]*Sf*fec)/sum(N[Iage,])
 }
  
# Check age structure
par(mfrow=c(2,1),oma=c(3,5,5,8),mar=c(5,4,1,1))
AgeS <- apply(N,1,sum)
plot(AgeRec+1:MaxAge,AgeS,xlab="Age",ylab="Relative abundance")
plot(AgeRec+1:MaxAge,FecAge,xlab="Age",ylab="Relative fecundity")

# find generation time
cat("Generation time =",Nsum1/Nsum2+AgeRec)