age.sim.var.normK <- function(m1, mdev, nstart, Klim, floor, it.len) {
    # simulates trajectory of age or stage-structured population,
    #   where fecundities of 3rd and 4th ages/stages normally distributed
    #   and survival rates normally distributed
    # Function arguments:
    #	m1: Age or stage-structured population matrix
    #   mdev: standard deviations: [F2, F3&4, P134, P23, P24] 
    #	nstart:   Starting population, vector of age class abundances
    #   Klim: population ceiling, summed over all age classes
    #   floor: population floor; when sum(nt)< floor, simulation ends
    # 	it.len:   length of iteration (e.g., 20, for 20 yr)
      zero <- 0.0000000000000000
      m1.mean <- m1
      m1.var <- m1.mean
      m1.var[1,2] <- rnorm(1, m1[1,2], mdev[1])  # age2 fecundity
      m1.var[1,3] <- rnorm(1, m1[1,3], mdev[2])  # age3 fecundity 
      m1.var[1,4] <- rnorm(1, m1[1,4], mdev[2])  # age4 fecundity
      m1.var[2,1] <- rnorm(1, m1[2,1], mdev[3])  # age1 survival 
      m1.var[4,3] <- rnorm(1, m1[4,3], mdev[3])  # age3 survival
      m1.var[4,4] <- rnorm(1, m1[4,4], mdev[3])  # age4 survival
      m1.var[3,2] <- rnorm(1, m1[3,2], mdev[4])  # survival juv->disp
      m1.var[4,2] <- rnorm(1, m1[4,2], mdev[5])  # survival juv->adult
      nt <- nstart 
      nt1 <- nstart
      ntraj <- nstart
      dim <- length(nstart)
      ntraj <- matrix(nrow=it.len, ncol=dim)
      j <- 0	
      while (j < it.len) {
          j <- j + 1
      nt <- m1.var %*% nt1
      nt <- trunc(nt)
      if (sum(nt) > Klim) {
          nt <- nt * Klim/sum(nt)
          }
      if (sum(nt) < floor) { 
    #                           nt <- zero
          j <- it.len
          }
      nt1 <- nt
      ntraj[j,] <- nt
      }
      ntraj
}

sim.age.pva.normK <- function(m1, mdev, nstart, Klim, floor, it.len, sim.num) { 
    # Performs "population viability analysis" for age-structured population.
    #     Created for population with 4 age/stage classes; easily modified 
    #        to other matrix dimensions. 
    #     Fecundities and survival rates normally distributed.  
    # Output: An array of the population age distributions at the end of each
    #     simulation.
    # Each row i contains number in each age class at end of the ith
    #     simulation.
    # A row containing all zeros implies extinction for that simulation.
    # More sophisticated versions would treat all matrix elements as variable.
    #     That project is left as an exercise. 
    # Function arguments: 
    #	  m1: Matrix of age/stage-specific fecundities and survival rates, 
    #         must be square matrix.
    #	  nstart:   starting population age distribution, 
    #               must be same dimension as matrix. 
    #     it.len:   length of iteration (e.g., 20, for 20 yr)
    #     sim.num:  number of simulations to run
    # Calls function "age.sim.var.normK", which simulates pop growth age 
    #   structure with normally distributed Fx and Px 
        col <- length(nstart)
 	sim.i <- NULL
	sim.res <- matrix(ncol=col)
	i <- 0
        while (i < sim.num) {
	    i <- i + 1
	    sim.end <- nstart
            sim.i <- age.sim.var.normK(m1, mdev, nstart, Klim, floor, it.len)
            sim.end <- sim.i[it.len,]
	    sim.res <- rbind(sim.res, sim.end)
        }
        sim.res
}