supthin = function(z,lambda,f,b=mean(lambda)){
## z = data, lambda = conditional intensity at pts, f = function to compute lambda, 
## and b = resulting rate.
## First thin, then superpose
keepz = list()
for(i in 1:z$n){
if(runif(1) < b/lambda[i]){
keepz$t = c(keepz$t,z$t[i])
keepz$lon = c(keepz$lon,z$lon[i])
keepz$lat = c(keepz$lat,z$lat[i])
}
}
candn = rpois(1,b*X1*Y1*T)
candt = sort(runif(candn)*T)
candx = runif(candn)*X1
candy = runif(candn)*Y1
for(i in 1:candn){
v = f(candt[i],candx[i],candy[i],z)
if(v < b){
if(runif(1) < (b-v)/b){
keepz$t = c(keepz$t,candt[i])
keepz$lon = c(keepz$lon,candx[i])
keepz$lat = c(keepz$lat,candy[i])
}}
}
keepz$lon = keepz$lon[order(keepz$t)]
keepz$lat = keepz$lat[order(keepz$t)]
keepz$t = sort(keepz$t)
keepz$n = length(keepz$t)
keepz
}