Mantel.test<-function(x,y,time,c1,c2,Nrep){
# -----------------------------------------------------------------
#  An example of execution
#
#     dat<-scan("d:/demo/McHardy.dat",list(x=0,y=0,time=0))
#     c1<-2; c2<-0; Nrep<-999
#     out<-KnoxM.test(dat$x,dat$y,dat$time,c1,c2,Nrep)
#
# ---------------------------------------------------------------------
#  ARGUEMENTS
#     x    : a vector of x-coordinates of case
#     y    : a vector of y-coordinates of case
#     time : a vector of observed times
#     c1   : a constant for Mantel's measure of closeness in space
#     c2   : a constant for Mantel's measure of closeness in time
#     Nrep : The number of Monte Carlo replications, e.g., 999, 9999
#
#  VALUES
#     Mantel.T          : test statistic
#     Expected          : expected value of test statistic
#     Simulated.p.value : simulated p-value
#     Freq              : a vector of simulated test statistics under the null
#
#-----------------------------------------------------------------------

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] <- 1/(sdis[i,j]+c1)
      at[i,j] <- 1/(tdis[i,j]+c2)
   }
   as[i,i]<-0
   at[i,i]<-0
}
s1<-0; s2<-0; s3<-0
for(i in 1:n){
   for (j in 1:n){
   s1<-s1+as[i,j] * at[i,j]
   s2<-s2+as[i,j]
   s3<-s3+at[i,j]  
   }
}
obst <- s1/2
Ms <- s2/2
Mt <- s3/2
ee  <- Ms*Mt/(n*(n-1)/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] <- 1/(tdis[i,j]+c2)
   }
   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(Mantel.T=obst, Expected=ee, Simulated.p.value=p,Freq=ktmon)
}