sim.data <- function(seed){

  if (!missing(seed)) set.seed(seed)

  q <- c(1.0,1.0,1.0)
  r <- 3.25
  theta <- c(0.20,0.25)
  sdpro <- c(0.35,0.40,0.20)
  sdobs <- c(1.0,0.9,0.75)
  ## Make z matrix for all years to make coding easy
  z.juv1 <- c(rep(c(1.4,0.7),each=10),rep(0.6,15))
  z.juv2 <- rep(0.4,35)
  z.adult <- c(rep(c(0.1,0.2),each=10),rep(0.4,15))
  z <- matrix(c(z.juv1,z.juv2,z.adult),nrow=3,ncol=35,byrow=TRUE)

  N.matrix <- data.frame(matrix(NA,ncol=40,nrow=3))
  names(N.matrix) <- paste('X',seq(1965,2004,1),sep='')

  N.matrix.1969 <- c(1.0,0.2,4.0)
  N.matrix$X1969 <- N.matrix.1969
  N.matrix[3,1:5] <- 4.0 # 4 million skates to start matrix calcs

  for (j in 6:40){
    N.matrix[1,j] <-  (N.matrix[3,j-5] * (r/2) + N.matrix[1,j-1] * (1 - theta[1])) * exp(-z[1,j-5]) # Deterministic 
    N.matrix[1,j] <- N.matrix[1,j]*exp(rnorm(1,0,sdpro[1]) - 0.5*sdpro[1]^2) #### ADD PROCESS ERROR
    
    N.matrix[2,j] <-  (N.matrix[1,j-1]  * theta[1] + N.matrix[2,j-1]  * (1 - theta[2]))* exp(-z[2,j-5])
    N.matrix[2,j] <- N.matrix[2,j]*exp(rnorm(1,0,sdpro[1]) - 0.5*sdpro[2]^2)
    
    N.matrix[3,j] <- (N.matrix[2,j-1]  * theta[2] + N.matrix[3,j-1])* exp(-z[3,j-5])
    N.matrix[3,j] <- N.matrix[3,j]*exp(rnorm(1,0,sdpro[3]) - 0.5*sdpro[3]^2)     
  }
  N.matrix <- N.matrix[,-seq(1,5,1)]
  N.obs <- N.matrix
  for (i in 1:nrow(N.obs)){
    for (j in 1:ncol(N.obs)){
      N.obs[i,j] <- N.obs[i,j]*exp(rnorm(1,0,sdobs[i]) - 0.5*sdobs[i]^2)
    }}
  final.dat <- data.frame(year=seq(1970,2004,1),t(N.obs))
  names(final.dat) <- c('year','juv1','juv2','adult')

### set up true.dat file
  true.dat <- data.frame(q1=q[1],q2=q[2],q3=q[3],Recruitment=r,theta.1=theta[1],theta.2=theta[2],sdpro1=sdpro[1],sdpro2=sdpro[2],sdpro3=sdpro[3],sdobs1=sdobs[1],sdobs2=sdobs[2],sdobs3=sdobs[3],z.juv1.1=unique(z.juv1)[1],z.juv1.2=unique(z.juv1)[2],z.juv1.3=unique(z.juv1)[3],z.juv2.1=unique(z.juv2),z.juv2.2=unique(z.juv2),z.juv2.3=unique(z.juv2),z.adult.1=unique(z.adult)[1],z.adult.2=unique(z.adult)[2],z.adult.3=unique(z.adult)[3],N.juv1.1969=N.matrix.1969[1],N.juv2.1969=N.matrix.1969[2],N.adult.1969=N.matrix.1969[3])

true.dat2 <- as.data.frame(t(true.dat))
names(true.dat2) <- 'value'
return(list(final.dat,true.dat2))
}
final.dat <- sim.data()[[1]]
true.dat <- sim.data()[[2]]
#===== data files
  dat.name <- 'skate.dat'
  bugdat.name <- 'skate_bugs.dat'
  true.dat.name <- 'true.dat'
  write.table(final.dat,row.names=F, col.names=T,  file=bugdat.name)
  write.table(true.dat,row.names=T, col.names=T,  file=true.dat.name)

  unlink(dat.name)
  zz <- file(description = dat.name, open = "w", blocking = TRUE, encoding = getOption("encoding"))
  cat( "# number of years of data", file = zz, sep = "\n")
  cat(nrow(final.dat), file = zz, sep = "\n")
  cat("# Recruitment lag", file = zz, sep = "\n")
  cat(5, file = zz, sep = "\n")
  cat('#Data', file = zz, sep = "\n")
  write.table(final.dat,row.names= FALSE, col.names=FALSE, append= TRUE, file=dat.name)