intlik2 <- function(parm,y=y,X=traplocs,delta=.3,ssbuffer=2){

   Xl <- min(X[,1]) - ssbuffer
   Xu <- max(X[,1]) + ssbuffer
   Yu <- max(X[,2]) + ssbuffer
   Yl <- min(X[,2]) - ssbuffer

   xg <- seq(Xl+delta/2,Xu-delta/2,delta)
   yg <- seq(Yl+delta/2,Yu-delta/2,delta)
   npix.x <- length(xg)
   npix.y <- plength(yg)

   G <- cbind(rep(xg,npix.y),sort(rep(yg,npix.x)))
   nG <- nrow(G)
   D <- e2dist(X,G)
   area<- (delta*delta)*nrow(G)

   # extract the parameters from the input vector
   alpha0 <- parm[1]
   alpha1 <- exp(parm[2])
   n0 <- exp(parm[3])  # note parm[3] lives on the real line
   probcap <- plogis(alpha0)*exp(-alpha1*D*D)
   Pm <- matrix(NA,nrow=nrow(probcap),ncol=ncol(probcap))
   ymat <- rbind(y,rep(0,ncol(y))) ## tack on all-zero history
   lik.marg <- rep(NA,nrow(ymat))
   for(i in 1:nrow(ymat)){
     Pm[1:length(Pm)] <- dbinom(rep(ymat[i,],nG), K,
                   probcap[1:length(Pm)], log=TRUE))
     lik.cond <- exp(colSums(Pm))
     lik.marg[i] <- sum( lik.cond*(1/nG) )
   }
   nv <- c(rep(1,length(lik.marg)-1),n0)
## part1 here is the combinatorial term,
##  log(factorial(N)) = lgamma(N+1)
   part1 <- lgamma(nrow(y)+n0+1) - lgamma(n0+1)
   part2 <- sum(nv*log(lik.marg))
   return( -1*(part1+ part2) )
}