###                                   Appendix S2                                  ###
###                                                                                ###
### Excerpt of code written in Program R to explore the effectiveness of localized ###
### management on the spread of chronic wasting disease in a geographically closed ###
### population of white-tailed deer.                                               ###

######################################################################################
################################## Workspace Set-Up ##################################
######################################################################################

# Clear worspace environment and import required packages
{
  rm(list=ls())
  
  require(CircStats)
  require(pdist)
  require(matrixStats)
  require(RColorBrewer)
  require(plotrix)
  library(doSNOW)
}

# Record scenario-specific variables
{
  ## At a density of 06, disp.dist=500, mate.dist=350, esta.dist=500, herd.cull=0.065
  ## At a density of 12, disp.dist=250, mate.dist=250, esta.dist=500, herd.cull=0.175
  ## At a density of 18, disp.dist=125, mate.dist=200, esta.dist=500, herd.cull=0.435
}

######################################################################################
################################## Model Definition ##################################
######################################################################################

# Define Monte Carlo Simulation
CWD_MonCar<-function(cull=F, area.size=20000, cell.size=250, deer.dens=6, 
                     half.life=0.5, k=1, eta=0.0001, inf=3, B.d=0.0055, B.i=0.00120, 
                     test.prop=0.20, hunt.antl=0.44, hunt.less=0.21, buff=500, 
                     fail1=0.25, fail2=0.10, disp.dist=250, mate.dist=250, 
                     esta.dist=500, cull.size=sqrt(3), cull.rate=0.8, herd.cull=0.175, 
                     imagery=F, trace=F) {
  
  # Define simulation parameter values
  {
    cell.numb=(area.size+buff+buff)/cell.size
    
    Cell.N=matrix(0,cell.numb,cell.numb)
    Cell.Z=matrix(0,cell.numb,cell.numb)
    Cell.V=matrix(0,cell.numb,cell.numb)
    
    ro=c(rep(NA,15+1),.010,.014,.019,.026,.035,.046,.060,.075,.093,.113,
         .135,.159,.185,.214,.247,.287,.340,.419,.567,1)
    om=c(rep(NA,23+35+1),.009,.013,.018,.025,.034,.045,.058,.073,.091,.110,
         .131,.154,.179,.205,.234,.267,.306,.356,.432,.577,1)
    mu.weekly=c(rep(.009,1+35+44),.010,.013,.015,.018,.021,.024,.028,.032,.037,.042,.047,
                .053,.059,.065,.071,.078,.085,.093,.101,.109,.117,.125,.134,
                .144,.154,.164,.176,.189,.203,.221,.244,.276,.320,.393,.544,1)
    mu=1-(1-c(0.20,0.10,0.05))^(1/52)
    mu=matrix(c((mu.weekly-min(mu.weekly))/max(mu.weekly-min(mu.weekly))*(1-mu[1])+mu[1],
                (mu.weekly-min(mu.weekly))/max(mu.weekly-min(mu.weekly))*(1-mu[2])+mu[2],
                (mu.weekly-min(mu.weekly))/max(mu.weekly-min(mu.weekly))*(1-mu[3])+mu[3]),
              nrow=3,byrow=T)
    for (i in 4:33) mu=rbind(mu,mu[3,])
    gm=1-log(2)/(52*half.life)
    
    norm.step=c(3.616552,256.011520)
    norm.turn=c(5.937467,0.2119131)
    disp.step=c(3.697977,1000)
    disp.turn=c(0.1157447,0.005709542)
    
    cull.numb=0; useable=T;
    det.date=NA; det.year=NA; det.prev=NA; x.min=NA; x.max=NA; y.min=NA; y.max=NA; d.min=NA; d.max=NA
  }
  
  # Set up data record
  {
    dat=data.frame(S=rep(NA,52*20),E=NA,I=NA,C=NA,
                   Norm=NA,Disp=NA,Mati=NA,Kill=NA,Hunt=NA,Sede=NA,Rear=NA,
                   GroupF=NA,GroupM=NA,Breed=NA,Prions=NA)
    t=1
  }
  
  # Generate starting population
  {
    ## Initialize female portion of the population
    {
      tmp.numb=round((area.size/1000)^2*deer.dens/6/2)
      deer.X=rep(NA,tmp.numb); deer.Y=rep(NA,tmp.numb); tmp=1
      deer.X[1]=sample((buff+300):(area.size+buff-300),1)
      deer.Y[1]=sample((buff+300):(area.size+buff-300),1)
      while (tmp<length(deer.X)) {
        tmp.X=sample((buff+300):(area.size+buff-300),1)
        tmp.Y=sample((buff+300):(area.size+buff-300),1)
        if (min(as.matrix(pdist(cbind(tmp.X,tmp.Y),
                                cbind(deer.X[!is.na(deer.X)],deer.Y[!is.na(deer.Y)]))))>300) {
          tmp=tmp+1; deer.X[tmp]=tmp.X; deer.Y[tmp]=tmp.Y
        }
      }
      
      tmp.deer=sample(5:9,tmp.numb,T)
      deer=data.frame(ID=1:sum(tmp.deer),Sex=0,Age=3,Cell=1,State="N",X=0,Y=0,
                      HR.Cen.X=0,HR.Cen.Y=0,HR.Dist.X=0,HR.Dist.Y=0,Group=0,GroupSex=0,GroupSize=0,
                      LeadX=0,LeadY=0,LeadA=0,OldGroup=0,Mother=0,Week=-1,Stage=1,Mort=0,Rate=0,Follow=0)
      levels(deer$State)=c("N","D","M","E","K","H","S","R")
      ### Normal, Dispersal, Mating, Exploratory, Killed, Hunted, Stationary, Rearing
      
      deer$Group[deer$Sex==0]=rep(1:tmp.numb,tmp.deer)
      deer$HR.Cen.X[deer$Sex==0]=rep(deer.X,tmp.deer)
      deer$HR.Cen.Y[deer$Sex==0]=rep(deer.Y,tmp.deer)
      deer$HR.Dist.X[deer$Sex==0]=rep(sample(100:300,tmp.numb,T),tmp.deer)
      deer$HR.Dist.Y[deer$Sex==0]=rep(sample(100:300,tmp.numb,T),tmp.deer)
    }
    
    ## Add male portion of the population
    {
      tmp.numb=round((area.size/1000)^2*deer.dens/10/2)
      deer.X=rep(NA,tmp.numb); deer.Y=rep(NA,tmp.numb); tmp=1
      deer.X[1]=sample((buff+450):(area.size+buff-450),1)
      deer.Y[1]=sample((buff+450):(area.size+buff-450),1)
      while (tmp<length(deer.X)) {
        tmp.X=sample((buff+450):(area.size+buff-450),1)
        tmp.Y=sample((buff+450):(area.size+buff-450),1)
        if (min(as.matrix(pdist(cbind(tmp.X,tmp.Y),
                                cbind(deer.X[!is.na(deer.X)],deer.Y[!is.na(deer.Y)]))))>450) {
          tmp=tmp+1; deer.X[tmp]=tmp.X; deer.Y[tmp]=tmp.Y
        }
      }
      
      tmp.deer=sample(5:15,tmp.numb,T)
      deer=rbind(deer,data.frame(ID=max(deer$ID)+1:sum(tmp.deer),Sex=1,Age=3,Cell=1,State="N",X=0,Y=0,
                                 HR.Cen.X=0,HR.Cen.Y=0,HR.Dist.X=0,HR.Dist.Y=0,Group=0,GroupSex=1,GroupSize=0,
                                 LeadX=0,LeadY=0,LeadA=0,OldGroup=0,Mother=0,Week=-1,Stage=1,Mort=0,Rate=0,Follow=0))
      deer$Group[deer$Sex==1]=rep(max(deer$Group)+1:tmp.numb,tmp.deer)
      deer$HR.Cen.X[deer$Sex==1]=rep(deer.X,tmp.deer)
      deer$HR.Cen.Y[deer$Sex==1]=rep(deer.Y,tmp.deer)
      deer$HR.Dist.X[deer$Sex==1]=rep(sample(150:450,tmp.numb,T),tmp.deer)
      deer$HR.Dist.Y[deer$Sex==1]=rep(sample(150:450,tmp.numb,T),tmp.deer)
    }
    
    ## Initialize all group characteristics to identical starting points
    {
      deer$X=deer$HR.Cen.X
      deer$Y=deer$HR.Cen.Y
      deer$LeadX=deer$X
      deer$LeadY=deer$Y
    }
    
    ## Define core area boundaries
    {
      core.area=c(min(deer$X),max(deer$X),min(deer$Y),max(deer$Y))
    }
  }
  
  # Visualize home range positioning
  {
    if (imagery) {
      par(mfrow=c(2,1),mar=c(1,1,1,1))
      plot(c(NA,NA),xlim=c(0,area.size+2*buff),ylim=c(0,area.size+2*buff),xlab="",ylab="",xaxt="n")
      draw.ellipse(deer$X[deer$Sex==0],deer$Y[deer$Sex==0],
                   deer$HR.Dist.X[deer$Sex==0],deer$HR.Dist.Y[deer$Sex==0])
      plot(c(NA,NA),xlim=c(0,area.size+2*buff),ylim=c(0,area.size+2*buff),xlab="",ylab="")
      draw.ellipse(deer$X[deer$Sex==1],deer$Y[deer$Sex==1],
                   deer$HR.Dist.X[deer$Sex==1],deer$HR.Dist.Y[deer$Sex==1])
    }
  }
  
  # Set up graphs
  {
    if (imagery) {
      par(mfrow = c(2, 2),         # 2x2 layout
          oma = c(2, 2.5, 0, 0.1), # one rows of text at the outer left and bottom margin and half at the right and top
          mar = c(1, 1, 2, 1),     # space for one row of text at ticks and to separate plots
          mgp = c(2, 1, 0),        # axis label at 2 rows distance, tick labels at 1 row
          xpd = NA)                # allow content to protrude into outer margin (and beyond)
      
      barplot(c(length(deer$ID[deer$Stage==1]),length(deer$ID[deer$Stage<37 & deer$Stage>1]),
                length(deer$ID[deer$Stage<81 & deer$Stage>36]),length(deer$ID[deer$Stage>=81])),
              names=c("S","E","I","C"))
      image(Cell.N,col=terrain.colors(100),xaxt="n",yaxt="n")
      image(Cell.Z,col=topo.colors(100),xaxt="n",yaxt="n")
      image(Cell.V,col=heat.colors(100),xaxt="n",yaxt="n")
      title(xlab = paste("STATUS AT YEAR 0 WEEK 0"), outer = TRUE, line = 1)
    }
  }
  
  # Run simulations
  for (year in 1:28) {
    
    ## PRE-RUT (SEP 1 - OCT 31)
    for (week in 1:9) {
      
      ### Movement
      {
        #### Update group size
        {
          deer=deer[order(deer$Group),]
          deer$GroupSize=rep(table(deer$Group),table(deer$Group))
        }
        
        #### Normal State
        {
          ##### Generate lead step lengths and turn angles
          {
            tmp.deer=unique(subset(deer,State=="N" | State=="R")[,8:17])
            tmp.numb=nrow(tmp.deer)
            
            tmp.step=rweibull(tmp.numb,norm.step[1],norm.step[2])
            tmp.turn=rwrpcauchy(tmp.numb,norm.turn[1],norm.turn[2])+tmp.deer$LeadA
            tmp.back=rbinom(tmp.numb,1,1-exp(-(((tmp.deer$LeadX-tmp.deer$HR.Cen.X)^2/(2*tmp.deer$HR.Dist.X^2))+
                                                 ((tmp.deer$LeadY-tmp.deer$HR.Cen.Y)^2/(2*tmp.deer$HR.Dist.Y^2)))))
            tmp.turn[tmp.back==1]=atan2(tmp.deer$HR.Cen.Y[tmp.back==1]-tmp.deer$LeadY[tmp.back==1],
                                        tmp.deer$HR.Cen.X[tmp.back==1]-tmp.deer$LeadX[tmp.back==1])
          }
          
          ##### Adjust lead position accordingly
          {
            deer$LeadX[deer$State=="N" | deer$State=="R"]=
            deer$LeadX[deer$State=="N"| deer$State=="R"]+
            rep(cos(tmp.turn)*tmp.step,tmp.deer$GroupSize)
            deer$LeadY[deer$State=="N" | deer$State=="R"]=
              deer$LeadY[deer$State=="N"| deer$State=="R"]+
              rep(sin(tmp.turn)*tmp.step,tmp.deer$GroupSize)
            deer$LeadA[deer$State=="N" | deer$State=="R"]=rep(tmp.turn,tmp.deer$GroupSize)
          }
          
          ##### Adjust follower position accordingly
          {
            tmp.dist=sqrt((deer$LeadY[deer$State=="N"]-deer$Y[deer$State=="N"])^2+
                            (deer$LeadX[deer$State=="N"]-deer$X[deer$State=="N"])^2)-
            rexp(length(deer$Y[deer$State=="N"]),1/225)
            tmp.turn=atan2(deer$LeadY[deer$State=="N"]-deer$Y[deer$State=="N"],
                           deer$LeadX[deer$State=="N"]-deer$X[deer$State=="N"])
            deer$X[deer$State=="N"]=deer$X[deer$State=="N"]+cos(tmp.turn)*tmp.dist
            deer$Y[deer$State=="N"]=deer$Y[deer$State=="N"]+sin(tmp.turn)*tmp.dist
          }
        }
        
        #### Dispersal State
        {
          
          tmp.deer=subset(deer,State=="D")
          tmp.numb=sum(deer$State=="D")
          
          if (tmp.numb>0) {
            tmp.step=rweibull(tmp.numb,disp.step[1],disp.step[2])
            tmp.turn=rwrpcauchy(tmp.numb,disp.turn[1],disp.turn[2])+tmp.deer$LeadA
            
            tmp.deer$LeadA=tmp.turn  
            tmp.deer$X=tmp.deer$X+cos(tmp.deer$LeadA)*tmp.step
            tmp.deer$Y=tmp.deer$Y+sin(tmp.deer$LeadA)*tmp.step
            
            for (i in 1:tmp.numb) {
              if (tmp.deer$Sex[i]==0)
                tmp.comp=subset(deer,State=="N" & Age>1 & GroupSex==0 & GroupSize<12 & Group!=tmp.deer$OldGroup[i])
              if (tmp.deer$Sex[i]==1)
                tmp.comp=subset(deer,State=="N" & Age>1 & GroupSex==1 & GroupSize<20 & Group!=tmp.deer$OldGroup[i])
              if (min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.comp$X,tmp.comp$Y))))<disp.dist) {
                tmp.find=which.min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.comp$X,tmp.comp$Y))))
                tmp.deer$State[i]="N"
                tmp.find=tmp.comp[tmp.find,]
                tmp.deer[i,8:17]=tmp.find[8:17] # New deer inherits group characteristics of nearest deer
                tmp.deer$OldGroup[i]=0
              } else if (min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.comp$X,tmp.comp$Y))))>esta.dist) {
                if (tmp.deer$X[i]>core.area[1] & tmp.deer$X[i]<core.area[2] & tmp.deer$Y[i]>core.area[1] & tmp.deer$Y[i]<core.area[2]) {
                  tmp.deer$State[i]="N"
                  tmp.deer$GroupSex[i]=tmp.deer$Sex[i]
                  tmp.deer$HR.Cen.X[i]=tmp.deer$X[i]
                  tmp.deer$HR.Cen.Y[i]=tmp.deer$Y[i]
                  if (tmp.deer$Sex[i]==0) {
                    tmp.deer$HR.Dist.X[i]=sample(100:300,1,T)
                    tmp.deer$HR.Dist.Y[i]=sample(100:300,1,T)
                  } else {
                    tmp.deer$HR.Dist.X[i]=sample(150:450,1,T)
                    tmp.deer$HR.Dist.Y[i]=sample(150:450,1,T)
                  }
                  tmp.deer$Group[i]=max(deer$Group,tmp.deer$Group)+1
                  tmp.deer$OldGroup[i]=0
                }
              }
            }
            deer[deer$State=="D",]=tmp.deer
          }
        }
        
        #### Rearing State
        {
          if (sum(deer$State=="R")>0) {
            tmp.X=rep(NA,length(deer$Y[deer$State=="R"]))
            tmp.Y=rep(NA,length(deer$Y[deer$State=="R"]))
            for (i in 1:sum(deer$State=="R")) {
              tmp.X[i]=deer$X[deer$ID==deer$Mother[deer$State=="R"][i]]
              tmp.Y[i]=deer$Y[deer$ID==deer$Mother[deer$State=="R"][i]]
            }
            tmp.dist=sqrt((tmp.Y-deer$Y[deer$State=="R"])^2+
                            (tmp.X-deer$X[deer$State=="R"])^2)-rexp(length(deer$Y[deer$State=="R"]),0.1)
            tmp.turn=atan2(tmp.Y-deer$Y[deer$State=="R"],
                           tmp.X-deer$X[deer$State=="R"])
            deer$X[deer$State=="R"]=deer$X[deer$State=="R"]+cos(tmp.turn)*tmp.dist
            deer$Y[deer$State=="R"]=deer$Y[deer$State=="R"]+sin(tmp.turn)*tmp.dist
          }
        }
        
        #### Turn Around Deer at Study Area Boundaries
        {
          deer$LeadA[deer$State=="D" & (deer$X<=0 | deer$X>=area.size+2*buff | deer$Y<=0 | deer$Y>=area.size+2*buff)]=
          deer$LeadA[deer$State=="D" & (deer$X<=0 | deer$X>=area.size+2*buff | deer$Y<=0 | deer$Y>=area.size+2*buff)]+pi/2
          deer$X[deer$X<=0]=1
          deer$X[deer$X>=area.size+2*buff]=area.size+2*buff-1
          deer$Y[deer$Y<=0]=1
          deer$Y[deer$Y>=area.size+2*buff]=area.size+2*buff-1
        }
      }
      
      ### Behavioral State Transition
      {
        deer$State[deer$State=="R" & deer$Age>26/52]="N"
        
        deer$State[deer$State=="N" & deer$Sex==0 & deer$Age>1][rbinom(
          length(deer$State[deer$State=="N" & deer$Sex==0 & deer$Age>1]),1,1-.99^(1/52))==1]="D"
        deer$State[deer$State=="N" & deer$Sex==1 & deer$Age>1][rbinom(
          length(deer$State[deer$State=="N" & deer$Sex==1 & deer$Age>1]),1,1-.90^(1/52))==1]="D"
        deer$OldGroup[deer$State=="D" & deer$OldGroup==0]=deer$Group[deer$State=="D" & deer$OldGroup==0]
        deer$Group[deer$State=="D"]=0
        
        deer=deer[order(deer$Group),]
      }
      
      ### Update grid-based locations and associated grid matrices
      {
        deer$Cell=ceiling(deer$X/cell.size)+cell.numb*(ceiling(deer$Y/cell.size)-1)
        
        Cell.N[,]=0
        Cell.N[sort(unique(deer$Cell))]=Cell.N[sort(unique(deer$Cell))]+table(deer$Cell)
        
        Cell.Z[,]=0
        Cell.Z[sort(unique(deer$Cell[deer$Stage>36]))]=
          Cell.Z[sort(unique(deer$Cell[deer$Stage>36]))]+table(deer$Cell[deer$Stage>36])
        
        Cell.Rate=(k*log(1+(B.d*Cell.Z+B.i*Cell.V)/k))/(1-eta+eta*Cell.N)
      }
      
      ### Transmission
      {
        deer$Rate=1-exp(-Cell.Rate[deer$Cell])
        deer$Stage[deer$Stage==1][rbinom(sum(deer$Stage==1),1,deer$Rate[deer$Stage==1])==1]=2
      }
      
      ### Disease Stage Transition
      {
        tmp.deer=subset(deer,deer$Stage>=17 & deer$Stage<37)
        deer$Stage[deer$Stage>=17 & deer$Stage<37][rbinom(nrow(tmp.deer),1,ro[tmp.deer$Stage])==1]=37
        
        tmp.deer=subset(deer,deer$Stage>=60 & deer$Stage<81)
        deer$Stage[deer$Stage>=60 & deer$Stage<81][rbinom(nrow(tmp.deer),1,om[tmp.deer$Stage])==1]=81
      }
      
      ### Host Survival
      {
        deer$Mort=mu[cbind(deer$Age+1,deer$Stage)]
        deer$State[rbinom(length(deer$Mort),1,1-deer$Mort)==0]="K"
        
        #### Dependent fawns
        {
          tmp.deer=subset(deer,Age<9/52)
          deer$State[deer$Age<9/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="K" & deer$Sex==0])]="K"
          tmp.deer=subset(deer,Age>=9/52 & Age<26/52)
          deer$State[deer$Age>=9/52 & deer$Age<26/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="K" & deer$Sex==0])]="N"
        }
      }
      
      ### Shedding of Prions
      {
        Cell.V[sort(unique(deer$Cell[deer$Stage>36 & deer$State=="K"]))]=
        Cell.V[sort(unique(deer$Cell[deer$Stage>36 & deer$State=="K"]))]+
        100*table(deer$Cell[deer$Stage>36 & deer$State=="K"])
        
        deer=subset(deer,State!="K")
        
        Cell.V[sort(unique(deer$Cell[deer$Stage>36]))]=
          Cell.V[sort(unique(deer$Cell[deer$Stage>36]))]+2*table(deer$Cell[deer$Stage>36])
      }
      
      ### Prion Survival
      {
        Cell.V[,]=rbinom(cell.numb*cell.numb,Cell.V,gm)
      }
      
      ### Age, Week, and Stage Increment
      {
        deer$Age=deer$Age+1/52
        deer$Week=deer$Week-1
        deer$Stage[deer$Stage>1]=deer$Stage[deer$Stage>1]+1
      }
      
      ### Update records
      {      
        dat[t,]=c(sum(deer$Stage==1),sum(deer$Stage>=2 & deer$Stage<=36),
                  sum(deer$Stage>=37 & deer$Stage<=80),sum(deer$Stage>=81),
                  sum(deer$State=="N"),sum(deer$State=="D"),sum(deer$State=="M"),
                  sum(deer$State=="K"),sum(deer$State=="H"),
                  sum(deer$State=="S"),sum(deer$State=="R"),
                  mean(table(deer$Group[deer$Group>0 & deer$Sex==0])),
                  mean(table(deer$Group[deer$Group>0 & deer$Sex==1])),0,sum(Cell.V))
        t=t+1
      }
    }
    
    ## Update graphs and check conditions
    {
      if (imagery) {
        barplot(c(length(deer$ID[deer$Stage==1]),length(deer$ID[deer$Stage<37 & deer$Stage>1]),
                  length(deer$ID[deer$Stage<81 & deer$Stage>36]),length(deer$ID[deer$Stage>=81])),
                names=c("S","E","I","C"))
        image(Cell.N,col=terrain.colors(100),xaxt="n",yaxt="n")
        image(Cell.Z,col=topo.colors(100),xaxt="n",yaxt="n")
        image(Cell.V,col=heat.colors(100),xaxt="n",yaxt="n")
        if (cull) {
          title(xlab = paste("STATUS AT YEAR, ",year,", WEEK ",week,": DETECTED",sep=""), outer=T, line=0)
        } else {
          title(xlab = paste("STATUS AT YEAR, ",year,", WEEK ",week,": UNDETECTED",sep=""), outer=T, line=0)
        }        
        title(xlab = paste("PREVALENCE = ",round(sum(Cell.Z)/sum(Cell.N),5)," WITH A GROWTH OF ",
                           round(((sum(Cell.Z)/sum(Cell.N))/
                                    ((dat$I[t-52]+dat$C[t-52])/
                                       (dat$S[t-52]+dat$E[t-52]+dat$I[t-52]+dat$C[t-52]))),3),
                           sep=""), outer=TRUE, line=1)
      }
      if (!cull & ((sum(Cell.Z)/sum(Cell.N)>.4) | (sum(Cell.Z+Cell.V)==0 & year>inf))) {
        useable=F
        break()
      }
    }
    
    ## Revert all dispersal deer to normal states for beginning of RUT
    {
      tmp.deer=subset(deer,State=="D")
      tmp.numb=sum(deer$State=="D")
      
      if (tmp.numb>0) {
        for (i in 1:tmp.numb) {
          tmp.find=which.min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(deer$X[deer$State=="N" & deer$Age>1 & deer$GroupSex==tmp.deer$Sex[i]],
                                                                                      deer$Y[deer$State=="N" & deer$Age>1 & deer$GroupSex==tmp.deer$Sex[i]]))))
          tmp.deer$State[i]="N"
          tmp.find=deer[deer$State=="N" & deer$Age>1 & deer$Sex==tmp.deer$Sex[i],][tmp.find,]
          tmp.deer[i,8:17]=tmp.find[8:17] # New deer inherits group characteristics of nearest deer
          tmp.deer$OldGroup=0
        }
      }
      
      deer[deer$State=="D",]=tmp.deer
      ###### deer=subset(deer,State!="D") #######
    }
    
    ## RUT (NOV 1 - DEC 31)
    for (week in 10:17) {
      
      ### Movement
      {
        #### Update group size
        {
          deer=deer[order(deer$Group),]
          deer$GroupSize=rep(table(deer$Group),table(deer$Group))
        }
        
        #### Normal State
        {
          ##### Females maintaining group cohesion
          {
            ###### Generate lead step lengths and turn angles
            {
              tmp.deer=unique(subset(deer,(State=="N" & Sex==0) | Age<1)[,8:17])
              tmp.numb=nrow(tmp.deer)
              
              tmp.step=rweibull(tmp.numb,norm.step[1],norm.step[2])
              tmp.turn=rwrpcauchy(tmp.numb,norm.turn[1],norm.turn[2])+tmp.deer$LeadA
              tmp.back=rbinom(tmp.numb,1,1-exp(-(((tmp.deer$LeadX-tmp.deer$HR.Cen.X)^2/(2*tmp.deer$HR.Dist.X^2))+
                                                   ((tmp.deer$LeadY-tmp.deer$HR.Cen.Y)^2/(2*tmp.deer$HR.Dist.Y^2)))))
              tmp.turn[tmp.back==1]=atan2(tmp.deer$HR.Cen.Y[tmp.back==1]-tmp.deer$LeadY[tmp.back==1],
                                          tmp.deer$HR.Cen.X[tmp.back==1]-tmp.deer$LeadX[tmp.back==1])
            }
            
            ###### Adjust lead position accordingly
            {
              deer$LeadX[(deer$State=="N" & deer$Sex==0) | deer$Age<1]=
              deer$LeadX[(deer$State=="N" & deer$Sex==0) | deer$Age<1]+
              rep(cos(tmp.turn)*tmp.step,tmp.deer$GroupSize)
              deer$LeadY[(deer$State=="N" & deer$Sex==0) | deer$Age<1]=
                deer$LeadY[(deer$State=="N" & deer$Sex==0) | deer$Age<1]+
                rep(sin(tmp.turn)*tmp.step,tmp.deer$GroupSize)
              deer$LeadA[(deer$State=="N" & deer$Sex==0) | deer$Age<1]=rep(tmp.turn,tmp.deer$GroupSize)
            }
            
            ###### Adjust follower position accordingly
            {
              tmp.dist=sqrt((deer$LeadY[deer$State=="N" & deer$Sex==0]-deer$Y[deer$State=="N" & deer$Sex==0])^2+
                              (deer$LeadX[deer$State=="N" & deer$Sex==0]-deer$X[deer$State=="N" & deer$Sex==0])^2)-
              rexp(length(deer$Y[deer$State=="N" & deer$Sex==0]),0.00406)
              tmp.turn=atan2(deer$LeadY[deer$State=="N" & deer$Sex==0]-deer$Y[deer$State=="N" & deer$Sex==0],
                             deer$LeadX[deer$State=="N" & deer$Sex==0]-deer$X[deer$State=="N" & deer$Sex==0])
              deer$X[deer$State=="N" & deer$Sex==0]=deer$X[deer$State=="N" & deer$Sex==0]+cos(tmp.turn)*tmp.dist
              deer$Y[deer$State=="N" & deer$Sex==0]=deer$Y[deer$State=="N" & deer$Sex==0]+sin(tmp.turn)*tmp.dist
            }
          }
          
          ##### Males abondoning group dynamics
          {
            ###### Generate step lengths and turn angles
            {
              tmp.deer=subset(deer,State=="N" & Sex==1 & Age>1)[,6:17]
              tmp.numb=nrow(tmp.deer)
              
              tmp.step=rweibull(tmp.numb,norm.step[1],norm.step[2])
              tmp.turn=rwrpcauchy(tmp.numb,norm.turn[1],norm.turn[2])+tmp.deer$LeadA
              tmp.back=rbinom(tmp.numb,1,1-exp(-(((tmp.deer$X-tmp.deer$HR.Cen.X)^2/(2*tmp.deer$HR.Dist.X^2))+
                                                   ((tmp.deer$Y-tmp.deer$HR.Cen.Y)^2/(2*tmp.deer$HR.Dist.Y^2)))))
              tmp.turn[tmp.back==1]=atan2(tmp.deer$HR.Cen.Y[tmp.back==1]-tmp.deer$Y[tmp.back==1],
                                          tmp.deer$HR.Cen.X[tmp.back==1]-tmp.deer$X[tmp.back==1])
            }
            
            ###### Adjust position accordingly
            {
              deer$X[deer$State=="N" & deer$Sex==1 & deer$Age>1]=
              deer$X[deer$State=="N" & deer$Sex==1 & deer$Age>1]+cos(tmp.turn)*tmp.step
              deer$Y[deer$State=="N" & deer$Sex==1 & deer$Age>1]=
                deer$Y[deer$State=="N" & deer$Sex==1 & deer$Age>1]+sin(tmp.turn)*tmp.step
              deer$LeadA[deer$State=="N" & deer$Sex==1 & deer$Age>1]=tmp.turn
            }
          }
        }
        
        #### Mating State
        {
          if (sum(deer$State=="M")>0) {
            for (i in 1:sum(deer$State=="M")) {
              deer$LeadX[deer$State=="M"][i]=deer$X[deer$ID==deer$Follow[deer$State=="M"][i]]
              deer$LeadY[deer$State=="M"][i]=deer$Y[deer$ID==deer$Follow[deer$State=="M"][i]]
            }
            tmp.dist=sqrt((deer$LeadY[deer$State=="M"]-deer$Y[deer$State=="M"])^2+
                            (deer$LeadX[deer$State=="M"]-deer$X[deer$State=="M"])^2)-
              rexp(length(deer$Y[deer$State=="M"]),0.1)
            tmp.turn=atan2(deer$LeadY[deer$State=="M"]-deer$Y[deer$State=="M"],
                           deer$LeadX[deer$State=="M"]-deer$X[deer$State=="M"])
            deer$X[deer$State=="M"]=deer$X[deer$State=="M"]+cos(tmp.turn)*tmp.dist
            deer$Y[deer$State=="M"]=deer$Y[deer$State=="M"]+sin(tmp.turn)*tmp.dist
          }
        }
        
        #### Rearing State
        {
          if (sum(deer$State=="R")>0) {
            tmp.X=rep(NA,length(deer$Y[deer$State=="R"]))
            tmp.Y=rep(NA,length(deer$Y[deer$State=="R"]))
            for (i in 1:sum(deer$State=="R")) {
              tmp.X[i]=deer$X[deer$ID==deer$Mother[deer$State=="R"][i]]
              tmp.Y[i]=deer$Y[deer$ID==deer$Mother[deer$State=="R"][i]]
            }
            tmp.dist=sqrt((tmp.Y-deer$Y[deer$State=="R"])^2+
                            (tmp.X-deer$X[deer$State=="R"])^2)-rexp(length(deer$Y[deer$State=="R"]),0.1)
            tmp.turn=atan2(tmp.Y-deer$Y[deer$State=="R"],
                           tmp.X-deer$X[deer$State=="R"])
            deer$X[deer$State=="R"]=deer$X[deer$State=="R"]+cos(tmp.turn)*tmp.dist
            deer$Y[deer$State=="R"]=deer$Y[deer$State=="R"]+sin(tmp.turn)*tmp.dist
          }
        }
        
        #### Turn Around Deer at Study Area Boundaries
        {
          deer$LeadA[deer$State=="D" & (deer$X<=0 | deer$X>=area.size+2*buff | deer$Y<=0 | deer$Y>=area.size+2*buff)]=
          deer$LeadA[deer$State=="D" & (deer$X<=0 | deer$X>=area.size+2*buff | deer$Y<=0 | deer$Y>=area.size+2*buff)]+pi/2
          deer$X[deer$X<=0]=1
          deer$X[deer$X>=area.size+2*buff]=area.size+2*buff-1
          deer$Y[deer$Y<=0]=1
          deer$Y[deer$Y>=area.size+2*buff]=area.size+2*buff-1
        }
      }
      
      ### Behavioral State Transition
      {
        deer$State[deer$State=="R" & deer$Age>26/52]="N"
        
        #### Mate Finding
        {
          tmp.free=subset(deer,Sex==0 & Age>1 & Week<0)     # Assumes female fawns do not breed
          if (sum(deer$Sex==0 & deer$Age>1 & deer$Week<0)>0) {
            tmp.deer=subset(deer,Sex==1 & Age>1 & State=="N") # Assumes male fawns do not breed
            tmp.dist=as.matrix(pdist(cbind(tmp.deer$X,tmp.deer$Y),cbind(tmp.free$X,tmp.free$Y)))
            
            if (min(tmp.dist)<mate.dist) {
              tmp.deer=tmp.deer[which(rowMins(tmp.dist)<mate.dist),]
              tmp.dist=as.matrix(pdist(cbind(tmp.deer$X,tmp.deer$Y),cbind(tmp.free$X,tmp.free$Y)))
              for (i in 1:nrow(tmp.deer)) {
                tmp.free=subset(deer,Sex==0 & Age>0 & Week<0)
                tmp.dist=as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.free$X,tmp.free$Y)))
                if (min(tmp.dist)<mate.dist) {
                  tmp.find=tmp.free$ID[which.min(tmp.dist)]
                  deer$Week[deer$ID==tmp.find]=1
                  deer$Follow[deer$ID==tmp.find]=tmp.deer$ID[i]
                  
                  tmp.deer$State[i]="M"
                  tmp.deer$Follow[i]=tmp.find
                  tmp.deer$Week[i]=deer$Week[deer$ID==tmp.find]
                  
                  deer[deer$ID==tmp.deer$ID[i],]=tmp.deer[i,]
                }
              }
            }
          }
        }
        
        #### Mating
        if (sum(deer$Week[deer$Sex==1]==0)>0) {
          
          tmp.deer=subset(deer,Sex==1 & Week==0)
          for (i in 1:nrow(tmp.deer)) {
            deer$Follow[deer$ID==tmp.deer$Follow[i]]=0
            deer$State[deer$ID==tmp.deer$Follow[i]]="N"
            deer$Week[deer$ID==tmp.deer$Follow[i]]=sample(27:31,1)
          }
          deer$Follow[deer$Sex==1 & deer$Week==0]=0
          deer$State[deer$Sex==1 & deer$Week==0]="N"
        }
        
        deer=deer[order(deer$Group),]
      }
      
      ### Update grid-based locations and associated grid matrices
      {
        deer$Cell=ceiling(deer$X/cell.size)+cell.numb*(ceiling(deer$Y/cell.size)-1)
        
        Cell.N[,]=0
        Cell.N[sort(unique(deer$Cell))]=Cell.N[sort(unique(deer$Cell))]+table(deer$Cell)
        
        Cell.Z[,]=0
        Cell.Z[sort(unique(deer$Cell[deer$Stage>36]))]=
          Cell.Z[sort(unique(deer$Cell[deer$Stage>36]))]+table(deer$Cell[deer$Stage>36])
        
        Cell.Rate=(k*log(1+(B.d*Cell.Z+B.i*Cell.V)/k))/(1-eta+eta*Cell.N)
      }
      
      ### Transmission
      {
        deer$Rate=1-exp(-Cell.Rate[deer$Cell])
        deer$Stage[deer$Stage==1][rbinom(sum(deer$Stage==1),1,deer$Rate[deer$Stage==1])==1]=2
      }
      
      ### Disease Stage Transition
      {
        tmp.deer=subset(deer,deer$Stage>=17 & deer$Stage<37)
        deer$Stage[deer$Stage>=17 & deer$Stage<37][rbinom(nrow(tmp.deer),1,ro[tmp.deer$Stage])==1]=37
        
        tmp.deer=subset(deer,deer$Stage>=60 & deer$Stage<81)
        deer$Stage[deer$Stage>=60 & deer$Stage<81][rbinom(nrow(tmp.deer),1,om[tmp.deer$Stage])==1]=81
      }
      
      ### Host Survival
      {
        #### Hunter Harvest
        {
          deer$State[deer$Sex==1 & deer$Age>1][
            rbinom(length(deer$State[deer$Sex==1 & deer$Age>1]),1,(1-hunt.antl)^(1/8))==0]="H"
          deer$State[deer$Sex==0 | deer$Age<1][
            rbinom(length(deer$State[deer$Sex==0 | deer$Age<1]),1,(1-hunt.less)^(1/8))==0]="H"
        }
        
        #### Sample for CWD
        {
          tmp=rbinom(length(deer$X[deer$State=="H" & deer$Stage>=37]),1,test.prop)
          if (sum(tmp)>0) {
            cull=T
            det.date=c(det.date,(year-1)*52+week)
            det.year=c(det.year,year)
            det.prev=c(det.prev,sum(deer$Stage>=37)/nrow(deer))
            x.min=c(x.min,deer$X[deer$State=="H" & deer$Stage>=37][which(tmp==1)[1]]-cull.size*1000)
            x.max=c(x.max,deer$X[deer$State=="H" & deer$Stage>=37][which(tmp==1)[1]]+cull.size*1000)
            y.min=c(y.min,deer$Y[deer$State=="H" & deer$Stage>=37][which(tmp==1)[1]]-cull.size*1000)
            y.max=c(y.max,deer$Y[deer$State=="H" & deer$Stage>=37][which(tmp==1)[1]]+cull.size*1000)
            d.min=c(d.min,week)
            d.max=c(d.max,week+8)
          }
        }
        
        #### Culling
        if (cull & year>inf & cull.size<7.5) {
          harv.numb=sum(deer$State=="H")
          for (i in 2:length(d.min)) {
            if (week>=d.min[i] & week<=d.max[i] & year==det.year[i]) {
              deer$State[deer$X<x.max[i] & deer$X>x.min[i] & 
                           deer$Y<y.max[i] & deer$Y>y.min[i]][rbinom(
                             length(deer$State[deer$X<x.max[i] & deer$X>x.min[i] & 
                                                 deer$Y<y.max[i] & deer$Y>y.min[i]]),1,
                             1-(1-cull.rate)^(1/8))==1]="H"
              cull.numb=cull.numb+sum(deer$State=="H")-harv.numb
            }
          }
        }
        
        if (cull & year>inf & cull.size>7.5) {
          harv.numb=sum(deer$State=="H")
          if (week>=d.min[2] & week<=d.max[2] & year==det.year[2]) {
            deer$State[rbinom(length(deer$State),1,
                              1-(1-herd.cull)^(1/8))==1]="H"
            cull.numb=cull.numb+sum(deer$State=="H")-harv.numb
          }
        }
        
        #### Dependent fawns (harvest mortality)
        {
          tmp.deer=subset(deer,Age<9/52)
          deer$State[deer$Age<9/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="H" & deer$Sex==0])]="K"
          tmp.deer=subset(deer,Age>=9/52 & Age<26/52)
          deer$State[deer$Age>=9/52 & deer$Age<26/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="H" & deer$Sex==0])]="N"
        }
        
        #### Following mates (harvest mortality)
        {
          tmp.deer=subset(deer,Sex==0)
          deer$Week[deer$Sex==0][which(tmp.deer$Follow %in% deer$ID[deer$State=="H" & deer$Sex==1])]=-1
          
          tmp.deer=subset(deer,Sex==1 & State=="M")
          deer$Week[deer$Sex==1 & deer$State=="M"][which(tmp.deer$Follow %in% deer$ID[deer$State=="H" & deer$Sex==0])]=-1
          deer$State[deer$Sex==1 & deer$State=="M"][which(tmp.deer$Follow %in% deer$ID[deer$State=="H" & deer$Sex==0])]="N"
        }
        
        #### Remove hunted individuals
        deer=subset(deer,State!="H")
        
        #### Regular mortality
        deer$Mort=mu[cbind(deer$Age+1,deer$Stage)]
        deer$State[rbinom(length(deer$Mort),1,1-deer$Mort)==0]="K"
        
        #### Dependent fawns (natural mortality)
        {
          tmp.deer=subset(deer,Age<9/52)
          deer$State[deer$Age<9/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="K" & deer$Sex==0])]="K"
          tmp.deer=subset(deer,Age>=9/52 & Age<26/52)
          deer$State[deer$Age>=9/52 & deer$Age<26/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="K" & deer$Sex==0])]="N"
        }
        
        #### Following mates (natural mortality)
        {
          tmp.deer=subset(deer,Sex==0)
          deer$Week[deer$Sex==0][which(tmp.deer$Follow %in% deer$ID[deer$State=="K" & deer$Sex==1])]=-1
          
          tmp.deer=subset(deer,Sex==1 & State=="M")
          deer$Week[deer$Sex==1 & deer$State=="M"][which(tmp.deer$Follow %in% deer$ID[deer$State=="K" & deer$Sex==0])]=-1
          deer$State[deer$Sex==1 & deer$State=="M"][which(tmp.deer$Follow %in% deer$ID[deer$State=="K" & deer$Sex==0])]="N"
        }
      }
      
      ### Shedding of Prions
      {
        Cell.V[sort(unique(deer$Cell[deer$Stage>36 & deer$State=="K"]))]=
        Cell.V[sort(unique(deer$Cell[deer$Stage>36 & deer$State=="K"]))]+
        100*table(deer$Cell[deer$Stage>36 & deer$State=="K"])
        
        deer=subset(deer,State!="K")
        
        Cell.V[sort(unique(deer$Cell[deer$Stage>36]))]=
          Cell.V[sort(unique(deer$Cell[deer$Stage>36]))]+2*table(deer$Cell[deer$Stage>36])
      }
      
      ### Prion Survival
      {
        Cell.V[,]=rbinom(cell.numb*cell.numb,Cell.V,gm)
      }
      
      ### Age, Day, and Stage Increment
      {
        deer$Age=deer$Age+1/52
        deer$Week=deer$Week-1
        deer$Stage[deer$Stage>1]=deer$Stage[deer$Stage>1]+1
      }
      
      ### Update records
      {
        dat[t,]=c(sum(deer$Stage==1),sum(deer$Stage>=2 & deer$Stage<=36),
                  sum(deer$Stage>=37 & deer$Stage<=80),sum(deer$Stage>=81),
                  sum(deer$State=="N"),sum(deer$State=="D"),sum(deer$State=="M"),
                  sum(deer$State=="K"),sum(deer$State=="H"),
                  sum(deer$State=="S"),sum(deer$State=="R"),
                  mean(table(deer$Group[deer$Group>0 & deer$Sex==0])),
                  mean(table(deer$Group[deer$Group>0 & deer$Sex==1])),0,sum(Cell.V))
        t=t+1
      }
    }
    
    ## Update graphs and check conditions
    {
      if (imagery) {
        barplot(c(length(deer$ID[deer$Stage==1]),length(deer$ID[deer$Stage<37 & deer$Stage>1]),
                  length(deer$ID[deer$Stage<81 & deer$Stage>36]),length(deer$ID[deer$Stage>=81])),
                names=c("S","E","I","C"))
        image(Cell.N,col=terrain.colors(100),xaxt="n",yaxt="n")
        image(Cell.Z,col=topo.colors(100),xaxt="n",yaxt="n")
        image(Cell.V,col=heat.colors(100),xaxt="n",yaxt="n")
        if (cull) {
          title(xlab = paste("STATUS AT YEAR, ",year,", WEEK ",week,": DETECTED",sep=""), outer=T, line=0)
        } else {
          title(xlab = paste("STATUS AT YEAR, ",year,", WEEK ",week,": UNDETECTED",sep=""), outer=T, line=0)
        }        
        title(xlab = paste("PREVALENCE = ",round(sum(Cell.Z)/sum(Cell.N),5)," WITH A GROWTH OF ",
                           round(((sum(Cell.Z)/sum(Cell.N))/
                                    ((dat$I[t-52]+dat$C[t-52])/
                                       (dat$S[t-52]+dat$E[t-52]+dat$I[t-52]+dat$C[t-52]))),3),
                           sep=""), outer=TRUE, line=1)
      }
      if (!cull & ((sum(Cell.Z)/sum(Cell.N)>.4) | (sum(Cell.Z+Cell.V)==0 & year>inf))) {
        useable=F
        break()
      }
    }
    
    ## Remove any following behavior and re-initiate group dynamics at the end of RUT
    {
      deer$Follow=0
      deer$Week[deer$Weel<1]=-1
      deer$Week[deer$Day>=31]=31
      deer$State[deer$State=="M"]="N"
      
      for (i in unique(deer$Group)) {
        deer$LeadX[deer$Group==i]=mean(deer$LeadX[deer$Group==i])
        deer$LeadY[deer$Group==i]=mean(deer$LeadY[deer$Group==i])
        deer$LeadA[deer$Group==i]=runif(1,-pi,pi)
      }
    }
    
    ## Record percentage of non-fawn females breeding
    {
      dat$Breed[t-1]=sum(deer$Week[deer$Sex==0 & deer$Age>1]>0)/
      sum(deer$Week[deer$Sex==0 & deer$Age>1]<999)
    }
    
    ## GESTATION (JAN 1 - APR 30)
    for (week in 18:35) {
      
      ### Movement
      {
        #### Update group size
        {
          deer=deer[order(deer$Group),]
          deer$GroupSize=rep(table(deer$Group),table(deer$Group))
        }
        
        #### Normal State
        {
          ##### Generate lead step lengths and turn angles
          {
            tmp.deer=unique(subset(deer,State=="N" | State=="R")[,8:17])
            tmp.numb=nrow(tmp.deer)
            
            tmp.step=rweibull(tmp.numb,norm.step[1],norm.step[2])
            tmp.turn=rwrpcauchy(tmp.numb,norm.turn[1],norm.turn[2])+tmp.deer$LeadA
            tmp.back=rbinom(tmp.numb,1,1-exp(-(((tmp.deer$LeadX-tmp.deer$HR.Cen.X)^2/(2*tmp.deer$HR.Dist.X^2))+
                                                 ((tmp.deer$LeadY-tmp.deer$HR.Cen.Y)^2/(2*tmp.deer$HR.Dist.Y^2)))))
            tmp.turn[tmp.back==1]=atan2(tmp.deer$HR.Cen.Y[tmp.back==1]-tmp.deer$LeadY[tmp.back==1],
                                        tmp.deer$HR.Cen.X[tmp.back==1]-tmp.deer$LeadX[tmp.back==1])
          }
          
          ##### Adjust lead position accordingly
          {
            deer$LeadX[deer$State=="N" | deer$State=="R"]=
            deer$LeadX[deer$State=="N"| deer$State=="R"]+
            rep(cos(tmp.turn)*tmp.step,tmp.deer$GroupSize)
            deer$LeadY[deer$State=="N" | deer$State=="R"]=
              deer$LeadY[deer$State=="N"| deer$State=="R"]+
              rep(sin(tmp.turn)*tmp.step,tmp.deer$GroupSize)
            deer$LeadA[deer$State=="N" | deer$State=="R"]=rep(tmp.turn,tmp.deer$GroupSize)
          }
          
          ##### Adjust follower position accordingly
          {
            tmp.dist=sqrt((deer$LeadY[deer$State=="N"]-deer$Y[deer$State=="N"])^2+
                            (deer$LeadX[deer$State=="N"]-deer$X[deer$State=="N"])^2)-
            rexp(length(deer$Y[deer$State=="N"]),1/225)
            tmp.turn=atan2(deer$LeadY[deer$State=="N"]-deer$Y[deer$State=="N"],
                           deer$LeadX[deer$State=="N"]-deer$X[deer$State=="N"])
            deer$X[deer$State=="N"]=deer$X[deer$State=="N"]+cos(tmp.turn)*tmp.dist
            deer$Y[deer$State=="N"]=deer$Y[deer$State=="N"]+sin(tmp.turn)*tmp.dist
          }
        }
        
        #### Dispersal State
        {
          tmp.deer=subset(deer,State=="D")
          tmp.numb=sum(deer$State=="D")
          
          if (tmp.numb>0) {
            tmp.step=rweibull(tmp.numb,disp.step[1],disp.step[2])
            tmp.turn=rwrpcauchy(tmp.numb,disp.turn[1],disp.turn[2])+tmp.deer$LeadA
            
            tmp.deer$LeadA=tmp.turn  
            tmp.deer$X=tmp.deer$X+cos(tmp.deer$LeadA)*tmp.step
            tmp.deer$Y=tmp.deer$Y+sin(tmp.deer$LeadA)*tmp.step
            
            for (i in 1:tmp.numb) {
              if (tmp.deer$Sex[i]==0)
                tmp.comp=subset(deer,State=="N" & Age>1 & GroupSex==0 & GroupSize<12 & Group!=tmp.deer$OldGroup[i])
              if (tmp.deer$Sex[i]==1)
                tmp.comp=subset(deer,State=="N" & Age>1 & GroupSex==1 & GroupSize<20 & Group!=tmp.deer$OldGroup[i])
              if (min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.comp$X,tmp.comp$Y))))<disp.dist) {
                tmp.find=which.min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.comp$X,tmp.comp$Y))))
                tmp.deer$State[i]="N"
                tmp.find=tmp.comp[tmp.find,]
                tmp.deer[i,8:17]=tmp.find[8:17] # New deer inherits group characteristics of nearest deer
                tmp.deer$OldGroup[i]=0
              } else if (min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.comp$X,tmp.comp$Y))))>esta.dist) {
                if (tmp.deer$X[i]>core.area[1] & tmp.deer$X[i]<core.area[2] & tmp.deer$Y[i]>core.area[1] & tmp.deer$Y[i]<core.area[2]) {
                  tmp.deer$State[i]="N"
                  tmp.deer$GroupSex[i]=tmp.deer$Sex[i]
                  tmp.deer$HR.Cen.X[i]=tmp.deer$X[i]
                  tmp.deer$HR.Cen.Y[i]=tmp.deer$Y[i]
                  if (tmp.deer$Sex[i]==0) {
                    tmp.deer$HR.Dist.X[i]=sample(100:300,1,T)
                    tmp.deer$HR.Dist.Y[i]=sample(100:300,1,T)
                  } else {
                    tmp.deer$HR.Dist.X[i]=sample(150:450,1,T)
                    tmp.deer$HR.Dist.Y[i]=sample(150:450,1,T)
                  }
                  tmp.deer$Group[i]=max(deer$Group,tmp.deer$Group)+1
                  tmp.deer$OldGroup[i]=0
                }
              }
            }
            deer[deer$State=="D",]=tmp.deer
          }
        }
        
        #### Rearing State
        {
          if (sum(deer$State=="R")>0) {
            tmp.X=rep(NA,length(deer$Y[deer$State=="R"]))
            tmp.Y=rep(NA,length(deer$Y[deer$State=="R"]))
            for (i in 1:sum(deer$State=="R")) {
              tmp.X[i]=deer$X[deer$ID==deer$Mother[deer$State=="R"][i]]
              tmp.Y[i]=deer$Y[deer$ID==deer$Mother[deer$State=="R"][i]]
            }
            tmp.dist=sqrt((tmp.Y-deer$Y[deer$State=="R"])^2+
                            (tmp.X-deer$X[deer$State=="R"])^2)-rexp(length(deer$Y[deer$State=="R"]),0.1)
            tmp.turn=atan2(tmp.Y-deer$Y[deer$State=="R"],
                           tmp.X-deer$X[deer$State=="R"])
            deer$X[deer$State=="R"]=deer$X[deer$State=="R"]+cos(tmp.turn)*tmp.dist
            deer$Y[deer$State=="R"]=deer$Y[deer$State=="R"]+sin(tmp.turn)*tmp.dist
          }
        }
        
        #### Turn Around Deer at Study Area Boundaries
        {
          deer$LeadA[deer$State=="D" & (deer$X<=0 | deer$X>=area.size+2*buff | deer$Y<=0 | deer$Y>=area.size+2*buff)]=
          deer$LeadA[deer$State=="D" & (deer$X<=0 | deer$X>=area.size+2*buff | deer$Y<=0 | deer$Y>=area.size+2*buff)]+pi/2
          deer$X[deer$X<=0]=1
          deer$X[deer$X>=area.size+2*buff]=area.size+2*buff-1
          deer$Y[deer$Y<=0]=1
          deer$Y[deer$Y>=area.size+2*buff]=area.size+2*buff-1
        }
      }
      
      ### Behavioral State Transition
      {
        deer$State[deer$State=="R" & deer$Age>26/52]="N"
        
        deer$State[deer$State=="N" & deer$Sex==0 & deer$Age>1][rbinom(
          length(deer$State[deer$State=="N" & deer$Sex==0 & deer$Age>1]),1,1-.99^(1/52))==1]="D"
        deer$State[deer$State=="N" & deer$Sex==1 & deer$Age>1][rbinom(
          length(deer$State[deer$State=="N" & deer$Sex==1 & deer$Age>1]),1,1-.90^(1/52))==1]="D"
        deer$OldGroup[deer$State=="D" & deer$OldGroup==0]=deer$Group[deer$State=="D" & deer$OldGroup==0]
        deer$Group[deer$State=="D"]=0
        
        deer=deer[order(deer$Group),]
      }
      
      ### Update grid-based locations and associated grid matrices
      {
        deer$Cell=ceiling(deer$X/cell.size)+cell.numb*(ceiling(deer$Y/cell.size)-1)
        
        Cell.N[,]=0
        Cell.N[sort(unique(deer$Cell))]=Cell.N[sort(unique(deer$Cell))]+table(deer$Cell)
        
        Cell.Z[,]=0
        Cell.Z[sort(unique(deer$Cell[deer$Stage>36]))]=
          Cell.Z[sort(unique(deer$Cell[deer$Stage>36]))]+table(deer$Cell[deer$Stage>36])
        
        Cell.Rate=(k*log(1+(B.d*Cell.Z+B.i*Cell.V)/k))/(1-eta+eta*Cell.N)
      }
      
      ### Transmission
      {
        deer$Rate=1-exp(-Cell.Rate[deer$Cell])
        deer$Stage[deer$Stage==1][rbinom(sum(deer$Stage==1),1,deer$Rate[deer$Stage==1])==1]=2
      }
      
      ### Disease Stage Transition
      {
        tmp.deer=subset(deer,deer$Stage>=17 & deer$Stage<37)
        deer$Stage[deer$Stage>=17 & deer$Stage<37][rbinom(nrow(tmp.deer),1,ro[tmp.deer$Stage])==1]=37
        
        tmp.deer=subset(deer,deer$Stage>=60 & deer$Stage<81)
        deer$Stage[deer$Stage>=60 & deer$Stage<81][rbinom(nrow(tmp.deer),1,om[tmp.deer$Stage])==1]=81
      }
      
      ### Host Survival
      {
        #### Culling
        if (cull & year>inf & cull.size<7.5) {
          harv.numb=sum(deer$State=="H")
          for (i in 2:length(d.min)) {
            if (week>=d.min[i] & week<=d.max[i] & year==det.year[i]) {
              deer$State[deer$X<x.max[i] & deer$X>x.min[i] & 
                           deer$Y<y.max[i] & deer$Y>y.min[i]][rbinom(
                             length(deer$State[deer$X<x.max[i] & deer$X>x.min[i] & 
                                                 deer$Y<y.max[i] & deer$Y>y.min[i]]),1,
                             1-(1-cull.rate)^(1/8))==1]="H"
              cull.numb=cull.numb+sum(deer$State=="H")-harv.numb
            }
          }
        }
        
        if (cull & year>inf & cull.size>7.5) {
          harv.numb=sum(deer$State=="H")
          if (week>=d.min[2] & week<=d.max[2] & year==det.year[2]) {
            deer$State[rbinom(length(deer$State),1,
                              1-(1-herd.cull)^(1/8))==1]="H"
            cull.numb=cull.numb+sum(deer$State=="H")-harv.numb
          }
        }
        
        #### Dependent fawns (harvest mortality)
        {
          tmp.deer=subset(deer,Age<9/52)
          deer$State[deer$Age<9/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="H" & deer$Sex==0])]="K"
          tmp.deer=subset(deer,Age>=9/52 & Age<26/52)
          deer$State[deer$Age>=9/52 & deer$Age<26/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="H" & deer$Sex==0])]="N"
        }
        
        #### Removed hunted individuals
        deer=subset(deer,State!="H")
        
        #### Regular mortality
        deer$Mort=mu[cbind(deer$Age+1,deer$Stage)]
        deer$State[rbinom(length(deer$Mort),1,1-deer$Mort)==0]="K"
        
        #### Dependent fawns (natural mortality)
        {
          tmp.deer=subset(deer,Age<9/52)
          deer$State[deer$Age<9/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="K" & deer$Sex==0])]="K"
          tmp.deer=subset(deer,Age>=9/52 & Age<26/52)
          deer$State[deer$Age>=9/52 & deer$Age<26/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="K" & deer$Sex==0])]="N"
        }
      }
      
      ### Shedding of Prions
      {
        Cell.V[sort(unique(deer$Cell[deer$Stage>36 & deer$State=="K"]))]=
        Cell.V[sort(unique(deer$Cell[deer$Stage>36 & deer$State=="K"]))]+
        100*table(deer$Cell[deer$Stage>36 & deer$State=="K"])
        
        deer=subset(deer,State!="K")
        
        Cell.V[sort(unique(deer$Cell[deer$Stage>36]))]=
          Cell.V[sort(unique(deer$Cell[deer$Stage>36]))]+2*table(deer$Cell[deer$Stage>36])
      }
      
      ### Prion Survival
      {
        Cell.V[,]=rbinom(cell.numb*cell.numb,Cell.V,gm)
      }
      
      ### Age, Day, and Stage Increment
      {
        deer$Age=deer$Age+1/52
        deer$Week=deer$Week-1
        deer$Stage[deer$Stage>1]=deer$Stage[deer$Stage>1]+1
      }
      
      ### Update records
      {      
        dat[t,]=c(sum(deer$Stage==1),sum(deer$Stage>=2 & deer$Stage<=36),
                  sum(deer$Stage>=37 & deer$Stage<=80),sum(deer$Stage>=81),
                  sum(deer$State=="N"),sum(deer$State=="D"),sum(deer$State=="M"),
                  sum(deer$State=="K"),sum(deer$State=="H"),
                  sum(deer$State=="S"),sum(deer$State=="R"),
                  mean(table(deer$Group[deer$Group>0 & deer$Sex==0])),
                  mean(table(deer$Group[deer$Group>0 & deer$Sex==1])),0,sum(Cell.V))
        t=t+1
      }
    }
    
    ## Update graphs and check conditions
    {
      if (imagery) {
        barplot(c(length(deer$ID[deer$Stage==1]),length(deer$ID[deer$Stage<37 & deer$Stage>1]),
                  length(deer$ID[deer$Stage<81 & deer$Stage>36]),length(deer$ID[deer$Stage>=81])),
                names=c("S","E","I","C"))
        image(Cell.N,col=terrain.colors(100),xaxt="n",yaxt="n")
        image(Cell.Z,col=topo.colors(100),xaxt="n",yaxt="n")
        image(Cell.V,col=heat.colors(100),xaxt="n",yaxt="n")
        if (cull) {
          title(xlab = paste("STATUS AT YEAR, ",year,", WEEK ",week,": DETECTED",sep=""), outer=T, line=0)
        } else {
          title(xlab = paste("STATUS AT YEAR, ",year,", WEEK ",week,": UNDETECTED",sep=""), outer=T, line=0)
        }        
        title(xlab = paste("PREVALENCE = ",round(sum(Cell.Z)/sum(Cell.N),5)," WITH A GROWTH OF ",
                           round(((sum(Cell.Z)/sum(Cell.N))/
                                    ((dat$I[t-52]+dat$C[t-52])/
                                       (dat$S[t-52]+dat$E[t-52]+dat$I[t-52]+dat$C[t-52]))),3),
                           sep=""), outer=TRUE, line=1)
      }
      if (!cull & ((sum(Cell.Z)/sum(Cell.N)>.4) | (sum(Cell.Z+Cell.V)==0 & year>inf))) {
        useable=F
        break()
      }
    }
    
    ## FAWNING (MAY 1 - AUG 31)
    for (week in 36:52) {
      
      ### Movement
      {
        #### Update group size
        {
          deer=deer[order(deer$Group),]
          deer$GroupSize=rep(table(deer$Group),table(deer$Group))
        }
        
        #### Normal State
        {
          ##### Generate lead step lengths and turn angles
          {
            tmp.deer=unique(subset(deer,State=="N" | State=="S" | State=="R")[,8:17])
            tmp.numb=nrow(tmp.deer)
            
            tmp.step=rweibull(tmp.numb,norm.step[1],norm.step[2])
            tmp.turn=rwrpcauchy(tmp.numb,norm.turn[1],norm.turn[2])+tmp.deer$LeadA
            tmp.back=rbinom(tmp.numb,1,1-exp(-(((tmp.deer$LeadX-tmp.deer$HR.Cen.X)^2/(2*tmp.deer$HR.Dist.X^2))+
                                                 ((tmp.deer$LeadY-tmp.deer$HR.Cen.Y)^2/(2*tmp.deer$HR.Dist.Y^2)))))
            tmp.turn[tmp.back==1]=atan2(tmp.deer$HR.Cen.Y[tmp.back==1]-tmp.deer$LeadY[tmp.back==1],
                                        tmp.deer$HR.Cen.X[tmp.back==1]-tmp.deer$LeadX[tmp.back==1])
          }
          
          ##### Adjust lead position accordingly
          {
            deer$LeadX[deer$State=="N" | deer$State=="S" | deer$State=="R"]=
            deer$LeadX[deer$State=="N"| deer$State=="S" | deer$State=="R"]+
            rep(cos(tmp.turn)*tmp.step,tmp.deer$GroupSize)
            deer$LeadY[deer$State=="N" | deer$State=="S" | deer$State=="R"]=
              deer$LeadY[deer$State=="N"| deer$State=="S" | deer$State=="R"]+
              rep(sin(tmp.turn)*tmp.step,tmp.deer$GroupSize)
            deer$LeadA[deer$State=="N" | deer$State=="S" | deer$State=="R"]=rep(tmp.turn,tmp.deer$GroupSize)
          }
          
          ##### Adjust follower position accordingly
          {
            tmp.dist=sqrt((deer$LeadY[deer$State=="N"]-deer$Y[deer$State=="N"])^2+
                            (deer$LeadX[deer$State=="N"]-deer$X[deer$State=="N"])^2)-
            rexp(length(deer$Y[deer$State=="N"]),1/225)
            tmp.turn=atan2(deer$LeadY[deer$State=="N"]-deer$Y[deer$State=="N"],
                           deer$LeadX[deer$State=="N"]-deer$X[deer$State=="N"])
            deer$X[deer$State=="N"]=deer$X[deer$State=="N"]+cos(tmp.turn)*tmp.dist
            deer$Y[deer$State=="N"]=deer$Y[deer$State=="N"]+sin(tmp.turn)*tmp.dist
          }
        }
        
        #### Dispersal State
        {
          tmp.deer=subset(deer,State=="D")
          tmp.numb=sum(deer$State=="D")
          
          if (tmp.numb>0) {
            tmp.step=rweibull(tmp.numb,disp.step[1],disp.step[2])
            tmp.turn=rwrpcauchy(tmp.numb,disp.turn[1],disp.turn[2])+tmp.deer$LeadA
            
            tmp.deer$LeadA=tmp.turn  
            tmp.deer$X=tmp.deer$X+cos(tmp.deer$LeadA)*tmp.step
            tmp.deer$Y=tmp.deer$Y+sin(tmp.deer$LeadA)*tmp.step
            
            for (i in 1:tmp.numb) {
              if (tmp.deer$Sex[i]==0)
                tmp.comp=subset(deer,State=="N" & Age>1 & GroupSex==0 & GroupSize<12 & Group!=tmp.deer$OldGroup[i])
              if (tmp.deer$Sex[i]==1)
                tmp.comp=subset(deer,State=="N" & Age>1 & GroupSex==1 & GroupSize<20 & Group!=tmp.deer$OldGroup[i])
              if (min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.comp$X,tmp.comp$Y))))<disp.dist) {
                tmp.find=which.min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.comp$X,tmp.comp$Y))))
                tmp.deer$State[i]="N"
                tmp.find=tmp.comp[tmp.find,]
                tmp.deer[i,8:17]=tmp.find[8:17] # New deer inherits group characteristics of nearest deer
                tmp.deer$OldGroup[i]=0
              } else if (min(as.matrix(pdist(cbind(tmp.deer$X[i],tmp.deer$Y[i]),cbind(tmp.comp$X,tmp.comp$Y))))>esta.dist) {
                if (tmp.deer$X[i]>core.area[1] & tmp.deer$X[i]<core.area[2] & tmp.deer$Y[i]>core.area[1] & tmp.deer$Y[i]<core.area[2]) {
                  tmp.deer$State[i]="N"
                  tmp.deer$GroupSex[i]=tmp.deer$Sex[i]
                  tmp.deer$HR.Cen.X[i]=tmp.deer$X[i]
                  tmp.deer$HR.Cen.Y[i]=tmp.deer$Y[i]
                  if (tmp.deer$Sex[i]==0) {
                    tmp.deer$HR.Dist.X[i]=sample(100:300,1,T)
                    tmp.deer$HR.Dist.Y[i]=sample(100:300,1,T)
                  } else {
                    tmp.deer$HR.Dist.X[i]=sample(150:450,1,T)
                    tmp.deer$HR.Dist.Y[i]=sample(150:450,1,T)
                  }
                  tmp.deer$Group[i]=max(deer$Group,tmp.deer$Group)+1
                  tmp.deer$OldGroup[i]=0
                }
              }
            }
            deer[deer$State=="D",]=tmp.deer
          }
        }
        
        #### Rearing State
        {
          if (sum(deer$State=="R")>0) {
            tmp.X=rep(NA,length(deer$Y[deer$State=="R"]))
            tmp.Y=rep(NA,length(deer$Y[deer$State=="R"]))
            for (i in 1:sum(deer$State=="R")) {
              tmp.X[i]=deer$X[deer$ID==deer$Mother[deer$State=="R"][i]]
              tmp.Y[i]=deer$Y[deer$ID==deer$Mother[deer$State=="R"][i]]
            }
            tmp.dist=sqrt((tmp.Y-deer$Y[deer$State=="R"])^2+
                            (tmp.X-deer$X[deer$State=="R"])^2)-rexp(length(deer$Y[deer$State=="R"]),0.1)
            tmp.turn=atan2(tmp.Y-deer$Y[deer$State=="R"],
                           tmp.X-deer$X[deer$State=="R"])
            deer$X[deer$State=="R"]=deer$X[deer$State=="R"]+cos(tmp.turn)*tmp.dist
            deer$Y[deer$State=="R"]=deer$Y[deer$State=="R"]+sin(tmp.turn)*tmp.dist
          }
        }
        
        #### Turn Around Deer at Study Area Boundaries
        {
          deer$LeadA[deer$State=="D" & (deer$X<=0 | deer$X>=area.size+2*buff | deer$Y<=0 | deer$Y>=area.size+2*buff)]=
          deer$LeadA[deer$State=="D" & (deer$X<=0 | deer$X>=area.size+2*buff | deer$Y<=0 | deer$Y>=area.size+2*buff)]+pi/2
          deer$X[deer$X<=0]=1
          deer$X[deer$X>=area.size+2*buff]=area.size+2*buff-1
          deer$Y[deer$Y<=0]=1
          deer$Y[deer$Y>=area.size+2*buff]=area.size+2*buff-1
        }
      }
      
      ### Behavioral State Transition
      { 
        #### Births
        {
          if (sum(deer$Week[deer$Sex==0 & deer$Age>0 & deer$State=="N"]==0)>0) {
            tmp.deer=subset(deer,Sex==0 & Age>0 & State=="N" & Week==0)
            tmp.deer=tmp.deer[order(tmp.deer$Age),]
            tmp.size=colSums(rmultinom(length(tmp.deer$Age[tmp.deer$Age<2.5]),1,
                                       prob=c(fail1,(1-fail1)/20*c(5,14,1)))*0:3)
            tmp.size=c(tmp.size,colSums(rmultinom(length(tmp.deer$Age[tmp.deer$Age>2.5]),
                                                  1,prob=c(fail2,(1-fail2)/20*c(5,14,1)))*0:3))
            if (sum(tmp.size)>0) {
              tmp.fawn=data.frame(ID=max(deer$ID)+1:sum(tmp.size),Sex=sample(c(0,1),sum(tmp.size),T),Age=0,Cell=1,State="S",
                                  X=rep(tmp.deer$X,tmp.size),Y=rep(tmp.deer$Y,tmp.size),
                                  HR.Cen.X=rep(tmp.deer$HR.Cen.X,tmp.size),HR.Cen.Y=rep(tmp.deer$HR.Cen.Y,tmp.size),
                                  HR.Dist.X=rep(tmp.deer$HR.Dist.X,tmp.size),HR.Dist.Y=rep(tmp.deer$HR.Dist.Y,tmp.size),
                                  Group=rep(tmp.deer$Group,tmp.size),GroupSex=0,GroupSize=rep(tmp.deer$GroupSize,tmp.size),
                                  LeadX=rep(tmp.deer$LeadX,tmp.size),LeadY=rep(tmp.deer$LeadY,tmp.size),
                                  LeadA=rep(tmp.deer$LeadA,tmp.size),OldGroup=0,Mother=rep(tmp.deer$ID,tmp.size),
                                  Week=3,Stage=1,Mort=0,Rate=0,Follow=0)
              deer$State[deer$Sex==0 & deer$Age>0 & deer$State=="N" & deer$Week==0]="S"
              deer$Week[deer$Sex==0 & deer$Age>0 & deer$State=="N" & deer$Week==0]=3
              deer[nrow(deer)+1:sum(tmp.size),]=tmp.fawn
            }
          }
        }
        
        #### Yearling Dispersal
        {
          deer$State[deer$State=="N" & deer$Sex==0 & abs(deer$Age-1)<.01][rbinom(
            length(deer$State[deer$State=="N" & deer$Sex==0 & abs(deer$Age-1)<.01]),1,0.2)==1]="D"
          deer$State[deer$State=="N" & deer$Sex==1 & abs(deer$Age-1)<.01]="D"
        }
        
        #### Stationary Cesation (in case of fawn death)
        {
          tmp.deer=subset(deer,State=="S" & Age>1)
          if (sum(deer$ID[deer$State=="S" & deer$Age>1])>0) {
            for (i in 1:nrow(tmp.deer)) {
              if (sum(deer$Mother==tmp.deer$ID[i])==0) {
                tmp.deer$State[i]="N"
                tmp.deer$Week[i]=-1
              }
            }
            deer[deer$State=="S" & deer$Age>1,]=tmp.deer
          }
        }
        
        #### Standard Transitions
        {
          deer$State[deer$State=="R" & deer$Age>26/52]="N"
          
          deer$State[deer$State=="S" & deer$Week==0 & deer$Age>.5]="N"
          deer$State[deer$State=="S" & deer$Week==0 & deer$Age<.5]="R"
          
          deer$State[deer$State=="N" & deer$Sex==0 & deer$Age>1][rbinom(
            length(deer$State[deer$State=="N" & deer$Sex==0 & deer$Age>1]),1,1-.99^(1/52))==1]="D"
          deer$State[deer$State=="N" & deer$Sex==1 & deer$Age>1][rbinom(
            length(deer$State[deer$State=="N" & deer$Sex==1 & deer$Age>1]),1,1-.90^(1/52))==1]="D"
          deer$OldGroup[deer$State=="D" & deer$OldGroup==0]=deer$Group[deer$State=="D" & deer$OldGroup==0]
          deer$Group[deer$State=="D"]=0
          
          deer=deer[order(deer$Group),]
        }
      }
      
      ### Update grid-based locations and associated grid matrices
      {
        deer$Cell=ceiling(deer$X/cell.size)+cell.numb*(ceiling(deer$Y/cell.size)-1)
        
        Cell.N[,]=0
        Cell.N[sort(unique(deer$Cell))]=Cell.N[sort(unique(deer$Cell))]+table(deer$Cell)
        
        Cell.Z[,]=0
        Cell.Z[sort(unique(deer$Cell[deer$Stage>36]))]=
          Cell.Z[sort(unique(deer$Cell[deer$Stage>36]))]+table(deer$Cell[deer$Stage>36])
        
        Cell.Rate=(k*log(1+(B.d*Cell.Z+B.i*Cell.V)/k))/(1-eta+eta*Cell.N)
      }
      
      ### Transmission
      {
        deer$Rate=1-exp(-Cell.Rate[deer$Cell])
        deer$Stage[deer$Stage==1][rbinom(sum(deer$Stage==1),1,deer$Rate[deer$Stage==1])==1]=2
      }
      
      ### Disease Stage Transition
      {
        tmp.deer=subset(deer,deer$Stage>=17 & deer$Stage<37)
        deer$Stage[deer$Stage>=17 & deer$Stage<37][rbinom(nrow(tmp.deer),1,ro[tmp.deer$Stage])==1]=37
        
        tmp.deer=subset(deer,deer$Stage>=60 & deer$Stage<81)
        deer$Stage[deer$Stage>=60 & deer$Stage<81][rbinom(nrow(tmp.deer),1,om[tmp.deer$Stage])==1]=81
      }
      
      ### Host Survival
      {
        deer$Mort=mu[cbind(deer$Age+1,deer$Stage)]
        deer$State[rbinom(length(deer$Mort),1,1-deer$Mort)==0]="K"
        
        #### Dependent fawns
        {
          tmp.deer=subset(deer,Age<9/52)
          deer$State[deer$Age<9/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="K" & deer$Sex==0])]="K"
          tmp.deer=subset(deer,Age>=9/52 & Age<26/52)
          deer$State[deer$Age>=9/52 & deer$Age<26/52][which(tmp.deer$Mother %in% deer$ID[deer$State=="K" & deer$Sex==0])]="N"
        }
      }
      
      ### Shedding of Prions
      {
        Cell.V[sort(unique(deer$Cell[deer$Stage>36 & deer$State=="K"]))]=
        Cell.V[sort(unique(deer$Cell[deer$Stage>36 & deer$State=="K"]))]+
        100*table(deer$Cell[deer$Stage>36 & deer$State=="K"])
        
        deer=subset(deer,State!="K")
        
        Cell.V[sort(unique(deer$Cell[deer$Stage>36]))]=
          Cell.V[sort(unique(deer$Cell[deer$Stage>36]))]+2*table(deer$Cell[deer$Stage>36])
      }
      
      ### Prion Survival
      {
        Cell.V[,]=rbinom(cell.numb*cell.numb,Cell.V,gm)
      }
      
      ### Age, Day, and Stage Increment
      {
        deer$Age=deer$Age+1/52
        deer$Week=deer$Week-1
        deer$Stage[deer$Stage>1]=deer$Stage[deer$Stage>1]+1
      }
      
      ### Update records
      {
        dat[t,]=c(sum(deer$Stage==1),sum(deer$Stage>=2 & deer$Stage<=36),
                  sum(deer$Stage>=37 & deer$Stage<=80),sum(deer$Stage>=81),
                  sum(deer$State=="N"),sum(deer$State=="D"),sum(deer$State=="M"),
                  sum(deer$State=="K"),sum(deer$State=="H"),
                  sum(deer$State=="S"),sum(deer$State=="R"),
                  mean(table(deer$Group[deer$Group>0 & deer$Sex==0])),
                  mean(table(deer$Group[deer$Group>0 & deer$Sex==1])),0,sum(Cell.V))
        t=t+1
      }
    }
    
    ## Introduce infected individual at study area boundary
    if (year==inf) {
      tmp.sick=data.frame(ID=max(deer$ID)+1,Sex=1,Age=2,Cell=1,State="D",
                          X=NA,Y=NA,HR.Cen.X=NA,HR.Cen.Y=NA,HR.Dist.X=200,HR.Dist.Y=200,
                          Group=0,GroupSex=1,GroupSize=NA,LeadX=NA,LeadY=NA,LeadA=runif(1,-pi,pi),
                          OldGroup=0,Mother=0,Week=0,Stage=37,Mort=0,Rate=0,Follow=0)
      if (sample(c(0,1),1)==0) {
        tmp.sick$X=sample(1:(area.size-1),1)
        tmp.sick$Y=sample(c(1,area.size-1),1)
      } else {
        tmp.sick$X=sample(c(1,area.size-1),1)
        tmp.sick$Y=sample(1:(area.size-1),1)
      }
      tmp.sick$HR.Cen.X=tmp.sick$X+sample(c(-1,1),1)
      tmp.sick$HR.Cen.Y=tmp.sick$Y+sample(c(-1,1),1)
      tmp.sick$LeadX=tmp.sick$X
      tmp.sick$LeadY=tmp.sick$Y
      deer[nrow(deer)+1,]=tmp.sick
    }
    
    ## Update graphs and check conditions
    {
      if (imagery) {
        barplot(c(length(deer$ID[deer$Stage==1]),length(deer$ID[deer$Stage<37 & deer$Stage>1]),
                  length(deer$ID[deer$Stage<81 & deer$Stage>36]),length(deer$ID[deer$Stage>=81])),
                names=c("S","E","I","C"))
        image(Cell.N,col=terrain.colors(100),xaxt="n",yaxt="n")
        image(Cell.Z,col=topo.colors(100),xaxt="n",yaxt="n")
        image(Cell.V,col=heat.colors(100),xaxt="n",yaxt="n")
        if (cull) {
          title(xlab = paste("STATUS AT YEAR, ",year,", WEEK ",week,": DETECTED",sep=""), outer=T, line=0)
        } else {
          title(xlab = paste("STATUS AT YEAR, ",year,", WEEK ",week,": UNDETECTED",sep=""), outer=T, line=0)
        }        
        title(xlab = paste("PREVALENCE = ",round(sum(Cell.Z)/sum(Cell.N),5)," WITH A GROWTH OF ",
                           round(((sum(Cell.Z)/sum(Cell.N))/
                                    ((dat$I[t-52]+dat$C[t-52])/
                                       (dat$S[t-52]+dat$E[t-52]+dat$I[t-52]+dat$C[t-52]))),3),
                           sep=""), outer=TRUE, line=1)
      }
      if (!cull & ((sum(Cell.Z)/sum(Cell.N)>.4) | (sum(Cell.Z+Cell.V)==0 & year>inf))) {
        useable=F
        break()
      }
    }
    
    ## Break loop if 11 years have passed since detection
    {
      if ((useable & cull) | (useable & cull.size==0)) {
        if (year>det.year[2]+11) {
          break()
        }
      }
    }
  }
  
  # Adjust data frame and produce final graph of results
  {
    dat$Z=dat$I+dat$C
    dat$N=dat$S+dat$E+dat$I+dat$C
    if (imagery) {
      par(mfrow=c(1,1))
      plot(1:(t-1),dat$N[1:(t-1)],ylim=c(0,max(dat$N[1:(t-1)])),
           type="l",lwd=2,xlab="Week",ylab="Abundance")
      points(1:(t-1),dat$Z[1:(t-1)],type="l",lwd=2,col="red")
    }
  }
  
  # Return projected values
  if (year==28 | !useable) {
    return(c(0,rep(NA,16)))
  }
  if (useable & (cull | cull.size==0)) {
    final.prev=dat$Z[det.date[2]-1+520]/dat$N[det.date[2]-1+520]
    start.prev=dat$Z[det.date[2]-1]/dat$N[det.date[2]-1]
    final.size=dat$N[det.date[2]-1+520]
    start.size=dat$N[det.date[2]-1]
    
    print(paste("Population growth rate =",round((final.size/start.size)^(1/10),3)))
    print(paste("Prevalence growth rate =",round((final.prev/start.prev)^(1/10),3)))
    print(paste("Final prevalence =",round(final.prev,5)))
    
    if (trace) 
      return(list(var=c(useable,cull,deer.dens,half.life,test.prop,cull.size,cull.rate,k,eta,
                        final.prev,start.prev,final.size,start.size,
                        det.date[2],cull.numb,length(det.date)-1,dat$Prions[det.date[2]-1]),
                  dat=dat))
    if (!trace)
      return(c(useable,cull,deer.dens,half.life,test.prop,cull.size,cull.rate,k,eta,
               final.prev,start.prev,final.size,start.size,
               det.date[2],cull.numb,length(det.date)-1,dat$Prions[det.date[2]-1]))
  }
}

######################################################################################
############################## Simulation Initialization #############################
######################################################################################

# Set up data frame, initialize time stamp, and set number of iterations and cores
{
  res=data.frame(useable=0,cull=0,deer.dens=0,half.life=0,test.prop=0,cull.size=0,
                 cull.rate=0,k=0,eta=0,final.prev=NA,start.prev=NA,final.size=NA,
                 start.size=NA,det.date=0,cull.numb=0,cull.times=0,prions=NA)
  TimIni=Sys.time()
  iter.max=250
  core=4
}

# Initialize multi-core processing
cl <- makeCluster(core, type="SOCK")
registerDoSNOW(cl)

# Loop until 250 useable iterations are achieved
while (sum(res$useable)<iter.max) {
  
  ## Concurrently run multiple iterations
  tmp=foreach(i=1:core, .combine=rbind, .packages=c('CircStats','pdist','matrixStats')) %dopar% (
    CWD_MonCar(test.prop=0.02, area.size=20000, deer.dens=18, cull.size=sqrt(99), mate.dist=200, disp.dist=125,
               B.d=0.0055, B.i=0.00120, cull.rate=0.0, herd.cull=0.435, half.life=0.5, inf=3))
  
  ## Temporarily remove final growth rates
  res$PrevGrowth=NULL
  res$SizeGrowth=NULL
  
  ## Add results to data frame
  tmp=as.data.frame(tmp); names(tmp)=names(res); res=rbind(res,tmp)
  
  ## Calculate final growth rates
  res$PrevGrowth=(res$final.prev/res$start.prev)^(1/10)
  res$SizeGrowth=(res$final.size/res$start.size)^(1/10)
  
  ## Check inclusion conditions
  res$useable[res$PrevGrowth>=1.45]=0
  res$useable[res$SizeGrowth< 0.85]=0
  res$useable[res$final.prev> 0.50]=0
  
  ## Plot final and start prevalence and population size
  boxplot(res$final.prev,res$start.prev,outline=F,names=c("Final","Start"))
  
  ## Display progress
  print(paste("Completion: ",round(sum(res$useable)/iter.max*100,4),"% at ",Sys.time(),sep=""))
  
  ## Save partially completed file
  write.csv(res,file=paste("Dens",mean(res$deer.dens[res$useable==1]),"Test",mean(res$test.prop[res$useable==1])*100,
                           ".Size",mean(res$cull.size[res$useable==1])^2,".Rate",mean(res$cull.rate[res$useable==1])*100,
                           ".Life",mean(res$half.life[res$useable==1]),".Part.txt",sep=""), row.names=F)
}

# Evaluate model speed
TimDif=Sys.time()-TimIni
TimDif

# Terminate multi-core processing
stopCluster(cl)

# Export model results
write.csv(res,file=paste("Dens",mean(res$deer.dens[res$useable==1]),"Test",mean(res$test.prop[res$useable==1])*100,
                         ".Size",mean(res$cull.size[res$useable==1])^2,".Rate",mean(res$cull.rate[res$useable==1])*100,
                         ".Life",mean(res$half.life[res$useable==1]),".txt",sep=""), row.names=F)