KnoxM.test<-function(x,y,time,del1,del2,Nrep){
# -----------------------------------------------------------------
#  An example of execution
#
#     dat<-scan("d:/demo/McHardy.dat",list(x=0,y=0,time=0))
#     del1<-2; del2<-0; Nrep<-999
#     out<-KnoxM.test(dat$x,dat$y,dat$time,del1,del2,Nrep)
#
# ---------------------------------------------------------------------
#  ARGUEMENTS
#     x    : a vector of x-coordinates of case
#     y    : a vector of y-coordinates of case
#     time : a vector of observed times
#     del1 : a measure of closeness in space
#     del2 : a measure of closeness in time
#     Nrep : The number of Monte Carlo replications, e.g., 999, 9999
#
#  VALUES
#     Knox.T            : test statistic
#     Freq              : a vector of simulated test statistics under the null
#     Simulated.p.value : simulated p-value
#
#-----------------------------------------------------------------------
n<-length(x)
ktmon<-1:(Nrep+1)
sdis<-matrix(0,n,n)
tdis<-matrix(0,n,n)
as<-sdis
at<-tdis
for (i in 1:n){
   for (j in 1:n){
     sdis[i,j]<-sqrt( (x[i]-x[j])^2 + (y[i]-y[j])^2 )
     tdis[i,j]<- abs( time[i] - time[j])
   }
}
for(i in 1:n){
   for (j in 1:n){
      as[i,j] <- ifelse(sdis[i,j]<= del1, 1, 0)
      at[i,j] <- ifelse(tdis[i,j]<= del2, 1, 0)
   }
   as[i,i]<-0
   at[i,i]<-0
}
s1<-0
for(i in 1:n){
   for (j in 1:n){
   s1<-s1+as[i,j] * at[i,j]
    }
}
obst <- s1/2
#
# Monte Carlo Hypothesis testing 
#
for(k in 1:Nrep){
timeR<-sample(time)

for (i in 1:n){
   for (j in 1:n){
     tdis[i,j]<- abs( timeR[i] - timeR[j])
   }
}
for(i in 1:n){
   for (j in 1:n){
      at[i,j] <- ifelse(tdis[i,j]<= del2, 1, 0)
   }
   at[i,i]<-0
}
s1<-0
for(i in 1:n){
   for (j in 1:n){
   s1<-s1+as[i,j] * at[i,j]
    }
}
ktmon[k] <- s1/2
}
ktmon[Nrep+1]<-obst
r<-length(ktmon[ktmon>=obst])
p<-r/(Nrep+1)
list(Knox.T=obst , Freq=ktmon, Simulated.p.value=p)
}