TangoFPD.test<-function(oo, ee, d, lam, s, Nrep){
# -----------------------------------------------------------------
#  An example of execution
#
#     d<-0:9 ; Nrep<-999 ; lam<- 2*(1:20), s<-0:4
#     ob<- c(0, 1, 3, 2, 1, 8, 7, 4, 0, 2)
#     ex<- c(0.041,0.701,1.128,1.634,1.868,5.738,4.550,3.617,3.792,4.931)
#     out<-TangoFPD.test(ob,ex,d,lam,Nrep)
#
# ---------------------------------------------------------------------
#  ARGUEMENTS
#     oo      : a vector of observed number of cases sorted in the ascending 
#                 order of the distance to the point source
#     ee      : a vector of conditional expected number of cases sorted 
#                 in the ascending order of the distance to the point source
#     d       : a vector of distances to the point source
#     lam     : a vector of values of lambda
#     s       : a vector of values of parameter s
#     Nrep    : The number of Monte Carlo replications, e.g., 999, 9999
#
#  VALUES
#     Test.statistic : a vector of test statistic T_decline for each of lambda 
#     Max.statistic  : max(T_decline)
#     p.value        : simulated p-value
#     Freq           : a vector of simulated test statistics under the null
#
#-----------------------------------------------------------------------

p<-ee/sum(ee)
q<-oo/sum(oo)
k<-length(s)
len<-length(lam)
dim<-matrix(0,len,k)
nn<-sum(oo)
nr<-length(oo)
tval0<- -200
for (i in 1:len){
ram<-lam[i]
for (j in 1:k){
t<-s[j]
if (t<=0.00001){
   g0<- 1
   g1<- exp( -4 * (d/ram)^2 )
}
if (t>=0.00001){
   g0<- 1 - 4*d*(d-t)/t/t 
   g1<- exp( -4 * ((d-t)/ram)^2 )
}
g<-ifelse(d<=t, g0, g1)
h1<- sum( g * (q-p) ) * nn
h2<- ( sum( g^2 * p ) - ( sum( g * p ) )^2 ) * nn
tval<- h1 / sqrt( h2  )
dim[i,j]<-tval
tval0<-max(tval, tval0)
}
}
tan<-tval0
#
nloop<-Nrep+1
for (i in 2:nloop){
aa<-rmultinom(1,nn,p)
q<-aa/nn
tval2<- -200
for (i in 1:len){
ram<-lam[i]
for (j in 1:k){
t<-s[j]
if (t<= 0.00001){
   g0<- 1 
   g1<- exp( -4 * (d/ram)^2 )
}
if (t>=0.00001){
   g0<- 1 - 4*d*(d-t)/t/t 
   g1<- exp( -4 * ((d-t)/ram)^2 )
}
g<-ifelse(d<=t, g0, g1)
h1<- sum( g * (q-p) ) * nn
h2<- ( sum( g^2 * p ) - ( sum( g * p ) )^2 ) * nn
tval<- h1/ sqrt( h2  )
tval2<-max(tval,tval2)
}
}
tan1<-tval2
tan<-c(tan,tan1)
}
p1<- rank(-tan)[1]/(Nrep+1)
list(Max.statistic=tan[1], Test.statistics=dim, p.value=p1, Freq=tan)
}