# Leslie Matrix from pgs. 51-52
t=50  # number of years
A=matrix(0,4,4)
A[1,] = c(0,  0,  2,  4)
A[2,1] = 0.5
A[3,2] = 1
A[4,3] = 0.5
A

# Build a matrix for storing population output from 'for loop'
# 'n' will store the number of individuals as follows:
#   - each row represents a year
#   - each column represents an age (stage) class
#   - value in row 2, column 3 = the number in yr 2 for age class 3
n=matrix(0,t,4)
n[1,]=c(100, 0, 0, 0)

# Project population forward and store output
for (yr in 2:t) {
    n[yr,] <- A%*%n[(yr-1),] # %*% = matrix multiply
}

# 'N' stores the population totals for each year 
N=apply(n,1,sum) # calculate the row sums of 'n'

# 'Lambda' will hold the annual growth rates
Lambda=N[2:t]/N[1:(t-1)]

# Store age (stage) distributions 
stage.structure=n/N   
# you can use the next line to see the rounded values if you like
#round(stage.structure,3)

# Make some plots
oldpar <- par(no.readonly = TRUE) # Store default graphics settings
par(mfrow=c(3, 1)) # Set graphics window to 3 rows of plots in 1 column
matplot(1:t,n,'o',xlab="Year",ylab="Population Size")
matplot(1:t,log(n),'o',xlab="Year",ylab="ln(Population Size)")
plot(1:(t-1),Lambda,'o', xlab="Year")

windows() # open new graphics window (IT WILL BE ON TOP OF THE OTHER GRAPH)
par(mfrow=c(2, 1)) # Set graphics window to 2 rows of plots in 1 column
matplot(1:10,n[1:10,],"o",xlab="Year",ylab="Population Size")
matplot(1:t,stage.structure,"o",xlab="Year",ylab="Proportion in Age Class")

par(oldpar) # default (original) graphics parameter settings restored

#========================================#
# New functions used this week include:  #
#   '%*%'                                #
#   'apply'                              #
#   'round'                              #
#   'matplot'                            #
#   'windows'                            #
#   'par'                                #
#========================================#