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)