FlowFishFishing: Simulations of coastal fisheries

#' 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))
}