`SCR.fn` <-
function(){
### This function simulates spatial capture-recapture data and fits the model
# described in Royle and Gardner (Camera trapping book chapter) 
# which is a binomial observation model (not multinomial)
# as described in Royle et al. (J. Appl. Ecology) 

library(spdep)
library(R2WinBUGS)

#Set sigma, lam0, and beta parameters
sigma=3
beta=1
lam0=.6
J=10          # number of capture periods

#Create a 10x10 grid of traps

x=4:13
y=4:13
locs=expand.grid(x,y)
ntraps=dim(locs)[1]


#Set larger area boundaries (0,0) to (18,18) for "buffer" area
llx=0
upx=18

plot(c(llx,upx), c(llx, upx), typ='n')
points(locs)

#Simulate N home range centers
N=28
bearx=runif(N, 0,upx)
beary=runif(N, 0,upx)

points(bearx, beary, col='red')

#Make the Dij matrix of distances from home center to trap locations
D=matrix(0, nrow=ntraps, ncol=N)
for(i in 1:N){
       D[ ,i] = spDistsN1(as.matrix(locs), c(bearx[i], beary[i]))
}

#Make the K matrix
k=matrix(0, nrow=ntraps, ncol=N)
k=exp(-D^2/sigma)

#to make the encounter history for each simulated bear
#first get at probability of encounter

Eo=matrix(NA, nrow=ntraps, ncol=N)
pmean=matrix(NA, nrow=ntraps, ncol=N)
for(i in 1:ntraps){
for(j in 1:N){
Eo[i,j] = lam0*k[i,j]
pmean[i,j] = 1 -exp(-Eo[i,j])    # Pr(capture) = 1-exp(-lam0*k[i,j]) , which differs from Royle et al. (J. Appl. Ecol.)
                                 ### k[i,j] = exp(-d[i,j]^2 / sigma)
}
}

#Now, make enounter history for each bear over T time steps

nij=matrix(0, nrow=ntraps, ncol=N)
for(i in 1:ntraps){
for(j in 1:N){
nij[i,j]=rbinom(1, J, pmean[i,j])
}
}

nij=t(nij)

##########Set up the variables for WinBUGS

Madd=100
M=Madd+N
aug=matrix(0, nrow=Madd, ncol=ntraps)

Yarray=as.matrix(rbind(nij, aug))

trapmat=as.matrix(locs)

xl=llx
yl=llx
yu=upx
xu=upx

data1<-list(y=Yarray, M=M, trapmat=trapmat, J=J, ntraps=ntraps,  xu=xu, xl=xl, yu=yu, yl=yl) 

sink("SCRmodel.txt")
cat("
model {

sigma~dunif(0, 15)
lam0~dgamma(.1,.1)
psi~dunif(0, 1)

for (i in 1:M){
	
      z[i]~dbern(psi)
	SX[i]~dunif(xl, xu)
	SY[i]~dunif(yl, yu)

   for(j in 1:ntraps) {
		D2[i,j] <- pow(SX[i]-trapmat[j,1], 2) + pow(SY[i]-trapmat[j,2],2)
		Eo[i,j] <- lam0*exp(-D2[i,j]/sigma)
		pmean[i,j]<-1-(exp(-Eo[i,j]))
		tmp[i,j]<-pmean[i,j]*z[i]
		y[i,j]~dbin(tmp[i,j],J)
		}
	}
N<-sum(z[1:M])
}
",fill=TRUE)
sink()

params = list('lam0', 'sigma', 'N', 'psi')

sigma=5
psi=.6
SX=as.vector(runif(M,2, 10))
SY=as.vector(runif(M,2, 10))

z=as.vector(rbinom(M, 1, .5))


inits =  function() {list(z=z, psi=psi, lam0=.2, SX=SX, SY=SY, sigma=sigma) }

fit = bugs(data1, inits, params, model.file="SCRmodel.txt",
       debug=FALSE, n.chains=3, n.iter=1200, n.burnin=1000, n.thin=3)
 
fit


}