#' Iterate the MLPA model
#'
#' @param params A list with all the biological and biophysical parameters
#' needed to run the model. Will normally be generated by makeMLPAparams()
#' @param num_times The number of time steps to iterate
#' @param ivec An integer vector with length num_times, with each element
#' specifying the index of the dispersal kernel (third dimension of params$kernel)
#' to use in each simulation year. Will normally be generated with sample();
#' included so that the same sequence of years can be anaylzed across model variants.
#' @param N0 Initial population abundance in each patch, in units of number of
#' individuals per unit of habitat. All age classes in all habitats are intialized
#' to the same abundance
#'
#' @return A list with two components. Nt is a three-dimensional array with the
#' abundance in each age class at each location at each time. The dimensions are
#' location, age, and time respecively. It includes time zero, so the third
#' dimension has length num_times+1.
#'
#' The second component is Wt, which is a two-dimensional array of spwaning
#' biomass, arranged by space and time.
MLPAiter <- function(params, num_times, ivec, N0=1) {
# Make the parameters locally available, and ensure that they don't persist
attach(params)
on.exit(detach(params))
# Get the number of sites
nsites <- dim(kernel)[1]
ws <- wx[-(1:(xM-1))]
# Set up the arrays to hold the results
Nt <- array(0, c(nsites, xmax, num_times+1))
Wt <- array(0, c(nsites, num_times+1))
Nt[,,1] <- N0
Wt[,1] <- drop(Nt[,-(1:(xM-1)) ,1]) %*% ws
for (ti in 1:num_times) {
# Choose the kernel
i <- ivec[ti]
CA <- drop(kernel[,,i])
St <- CA %*% Wt[,ti]
Rt <- alpha*St/(1 + beta*St)
Nt[,-1,ti+1] <- p*Nt[,-xmax,ti]
Nt[,1,ti+1] <- Rt
Wt[,ti+1] <- Nt[, -(1:(xM-1)), ti+1] %*% ws
}
return(list(Nt=Nt, Wt=Wt))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.