alphaG <- function(vario.param, alphapval = 0.05, D = 1)
  {
    ax <- vario.param$range
    modcov <- vario.param$typmod

    if (modcov==3)
      {
        A <- pi*ax*ax/5
      }
    if (modcov==4)
      {
        A <- 2*pi*ax*ax
      }
    if (modcov==7)
      {
        A <- pi*ax*ax
      }
    if (modcov!=3 && modcov!=4 && modcov!=7)
      stop("this function is not implemented for this covariance model")
    N <- D/A
    alpha <- 1 - (1-alphapval)^(1/N)
  }

############################################################################

as.vario.param <- function(model,range,sill,nugget=0)
{
v <- c(model,range,sill,nugget)
  
if (v[1]=="linear") typmod <- 1
if (v[1]=="circular") typmod <- 2
if (v[1]=="spherical") typmod <- 3
if (v[1]=="exponential") typmod <- 4
if (v[1]=="pentaspherical") typmod <- 5
if (v[1]=="power") typmod <- 6
if (v[1]=="gaussian") typmod <- 7
if (v[1]=="bessel") typmod <- 8
if (v[1]=="hole bessel") typmod <- 9
if (v[1]=="pepitic") typmod <- 10

param <- list(typmod=typmod,range=as.numeric(v[2]),sill=as.numeric(v[3]),nugget=as.numeric(v[4]))
class(param) <- "VarioParam"
return(param)
}

############################################################################

buildcovgrad <- function(x0, y0, x, y, path=myfortranpath, typmod = 4, ax = 0.1, pal = 1, pep=0)
{
##                  typmod =  1 for bounded linear
##                  typmod =  2 for circular
##                  typmod =  3 for spherical
##                  typmod =  4 for exponential
##                  typmod =  5 for pentaspherical
##                  typmod =  6 for power function
##                  typmod =  7 for gaussian
##                  typmod =  8 for bessel function
##                  typmod =  9 for 2d hole effect with bessel function
##                  typmod = 10 pepitic

  np  <- length(x)
  dx0 <- as.vector(outer(x0, x, FUN = "-"))
  dy0 <- as.vector(outer(y0, y, FUN = "-"))
  dd <- sqrt(dx0^2 + dy0^2)
  tt <- c(ax, pep, pal, coef=1)
  cc <- rep(0,  np )
  gs <- rep(1,  np )

  pathname <- paste(path,"zac.so",sep="")
  dyn.load(pathname)
  cc <- 0 - .Fortran("calgderiv",
                     gs = double( np ),
                     as.integer(1),
                     as.double(dd),
                     as.integer( np ),
                     as.integer(typmod),
                     as.double(tt))$gs
  c1 <- (cc * dx0)/dd
  c2 <- (cc * dy0)/dd
  cd <- cbind(c1, c2)
  return(cd)
}

############################################################################

buildcovhess <- function(x0, y0, x, y, path=myfortranpath, typmod = 4, ax = 0.1, pal = 1, pep=0)
{
##                  typmod =  1 for bounded linear
##                  typmod =  2 for circular
##                  typmod =  3 for spherical
##                  typmod =  4 for exponential
##                  typmod =  5 for pentaspherical
##                  typmod =  6 for power function
##                  typmod =  7 for gaussian
##                  typmod =  8 for bessel function
##                  typmod =  9 for 2d hole effect with bessel function
##                  typmod = 10 pepitic

  np <- length(x)
  dx0 <- as.vector(outer(x, x0, FUN = "-"))
  dy0 <- as.vector(outer(y, y0, FUN = "-"))
  dd <- sqrt(dx0^2 + dy0^2)
  dx0 <- dx0/dd
  dy0 <- dy0/dd	
  tt <- c(ax, pep, pal, coef=1)		
  cch <- rep(0, np)
  gsh <- rep(1, np)

  pathname <- paste(path,"zac.so",sep="")
  dyn.load(pathname)
  cch <- 0 - .Fortran("calgd2f",
                      gsh = double(np),
                      as.integer(1),
                      as.double(dd),
                      as.integer(np),
                      as.integer(typmod),
                      as.double(tt))$gsh
  ccg <- rep(0, np)
  gsg <- rep(1, np)

  pathname <- paste(path,"zac.so",sep="")
  dyn.load(pathname)
  ccg <- 0 - .Fortran("calgderiv",
                      gsg = double(np),
                      as.integer(1),
                      as.double(dd),
                      as.integer(np),
                      as.integer(typmod),
                      as.double(tt))$gsg

  c1 <- cch * dx0 * dx0 + (ccg * dy0 * dy0)/dd
  c2 <- (cch - ccg/dd) * dx0 * dy0
  c3 <- (cch - ccg/dd) * dy0 * dx0
  c4 <- cch * dy0 * dy0 + (ccg * dx0 * dx0)/dd
  cd <- cbind(c1, c2, c3, c4)

  return(cd)
}

############################################################################

buildcovmat <- function(x, y, path=myfortranpath, typmod = 4, ax = 0.1, pal = 1, pep=0)
{
##                  typmod =  1 for bounded linear
##                  typmod =  2 for circular
##                  typmod =  3 for spherical
##                  typmod =  4 for exponential
##                  typmod =  5 for pentaspherical
##                  typmod =  6 for power function
##                  typmod =  7 for gaussian
##                  typmod =  8 for bessel function
##                  typmod =  9 for 2d hole effect with bessel function
##                  typmod = 10 pepitic

  np <- length(x)
  dd <- sqrt(outer(x, x, FUN = "-")^2 + outer(y, y, FUN = "-")^2)
  tt <- c(ax, pep, pal, coef=1)	
  cc <- matrix(1, np, np)
  gs <- rep(1, np)

  pathname <- paste(path,"zac.so",sep="")
  dyn.load(pathname)
  for(i in 1:np)
    {
      cc[i,  ] <- 1 - .Fortran("calgs",
                               gs = double(np),
                               as.integer(1),
                               as.double(dd[i,  ]),
                               as.integer(np),
                               as.integer(typmod),
                               as.double(tt))$gs
    }
  return(cc)
}

############################################################################

buildcovvec<- function(x0, y0, x, y, path=myfortranpath, typmod = 4, ax = 0.1, pal = 1, pep = 0)
{
##                  typmod =  1 for bounded linear
##                  typmod =  2 for circular
##                  typmod =  3 for spherical
##                  typmod =  4 for exponential
##                  typmod =  5 for pentaspherical
##                  typmod =  6 for power function
##                  typmod =  7 for gaussian
##                  typmod =  8 for bessel function
##                  typmod =  9 for 2d hole effect with bessel function
##                  typmod = 10 pepitic

  np <- length(x)
  dd <- as.vector(sqrt(outer(x0, x, FUN = "-")^2 + outer(y0, y, FUN = "-")^2))
  tt <- c(ax, pep, pal, coef=1)
  cc <- rep(0, np)
  gs <- rep(1, np)

  pathname <- paste(path,"zac.so",sep="")
  dyn.load(pathname)
  cc <- 1 - .Fortran("calgs",
                     gs = double(np),
                     as.integer(1),
                     as.double(dd),
                     as.integer(np),
                     as.integer(typmod),
                     as.double(tt))$gs
  return(cc)
}

############################################################################

cc <- function(X2, talpha, path=myfortranpath)
{
  nx <- nrow(X2)
  ny <- ncol(X2)
  bc <- as.integer(na2zero(X2) > talpha)
  bc <- as.vector(bc)
  lab <- rep(0, nx * ny)
  kvol <- 0
  ncc <- 0
  ncmax1 <- 0
  ncmax2 <- 0
  ncmax3 <- 0

  pathname <- paste(path,"zac.so",sep="")
  dyn.load(pathname)
  cc <- .Fortran("cc",
                 as.integer(nx),
                 as.integer(ny),
                 as.integer(bc),
                 as.integer(lab),
                 as.integer(kvol),
                 as.integer(ncc),
                 as.integer(ncmax1),
                 as.integer(ncmax2),
                 as.integer(ncmax3))
  
  invisible(return(list(ncc=cc[[6]], lab=cc[[4]])))
}

############################################################################

elim <- function(reszac, only.sign = T)
{
  x <- reszac$x
  y <- reszac$y
  z <- reszac$z
  xygrid <- reszac$xygrid
  cclab <- reszac$cclab
  p.value <- reszac$p.value
  alphapval <- reszac$alphapval
  
  xgrid <- xygrid$x
  ygrid <- xygrid$y
  np <- length(x)
  xx <- x
  yy <- y
  zz <- z

  ncc <- length(p.value)

  PL <- NULL

  if (only.sign == T)
    {
      for(nc in 1:ncc)
        {
          if (p.value[nc] <= alphapval)
            {
              coord <- which(cclab==nc,arr.ind=T)
              coordi <- coord[,1]
              coordj <- coord[,2]

              if (all(coordi==coordi[1])==F && all(coordj==coordj[1])==F)
                {
                  p <- chull(coordi,coordj)
                  pl <- list(x = xgrid[coordi[p]],y = ygrid[coordj[p]])
                  if (length(p.value) > 1)
                    { 
                      M<-list(pl)
                      PL<-c(PL,M)
                    }
                  if (length(p.value) == 1) PL<- list(pl)
                }
              else
               {
	          pl <- list(x=coordi,y=coordj)	
                  if (length(p.value) > 1)
                    { 
                      M<-list(pl)
                      PL<-c(PL,M)
                    }
                  if (length(p.value) == 1) PL<- list(pl)		
               }		
            }
        }
    }
  
  else
    {
      for(nc in 1:ncc)
        {
          coord <- which(cclab==nc,arr.ind=T)
          coordi <- coord[,1]
          coordj <- coord[,2]  
          
          if (all(coordi==coordi[1])==F && all(coordj==coordj[1])==F)
            {
              p <- chull(coordi,coordj)
              pl <- list(x = xgrid[coordi[p]],y = ygrid[coordj[p]]) 
              if (length(p.value) > 1)
                { 
                  M<-list(pl)
                  PL<-c(PL,M)
                }
              if (length(p.value) == 1) PL<- list(pl)
            }
         }
    }

  if (length(PL) == 1) ncc <- 1
  else
    {
      ncc <- length(PL)
    }

  for (i in 1:np)
    {
      for(nc in ncc)
        {
         pip <- inout(pts=rbind(c(x[i],y[i]),c(x[i],y[i])),poly=cbind(PL[[nc]]$x,PL[[nc]]$y),bound=NULL,quiet=TRUE)[1]
          if ( pip == T )
            {
              xx[i] <- NA
              yy[i] <- NA 	
              zz[i] <- NA
            }
        }
    }
  x <- xx[!(is.na(xx))]
  y <- yy[!(is.na(yy))]
  z <- zz[!(is.na(zz))]

  npelim <- sum(is.na(xx))
  np <- length(x)

  X <- NULL
  Y <- NULL
  Z <- NULL
  XX <- NULL
  YY <- NULL
  ZZ <- NULL
  xelim <- NULL
  yelim <- NULL
  zelim <- NULL
  ncelim <- 0

  for (i in 1:(np-1))
    { 
      for (j in (i+1):np)
        {
          RES <- NULL
          for(nc in ncc)
            {
              l <- length(PL[[nc]]$x)

	      if (l > 1)
                {
		for(m in 1:(l-1))
		      { 
		      res <- intersection(x[i],y[i],x[j],y[j],PL[[nc]]$x[m],PL[[nc]]$y[m],PL[[nc]]$x[m+1],PL[[nc]]$y[m+1])
		      RES <- c(RES,res)
		      }
		}	      
             }
          if (sum(RES)==0)
            {
              X <- c(X,x[i])
              Y <- c(Y,y[i])
              Z <- c(Z,z[i])
              XX <- c(XX,x[j])
              YY <- c(YY,y[j])
              ZZ <- c(ZZ,z[j])
            }
          else
            {
              ncelim <- ncelim+1
              xelim <- c(xelim,x[i],x[j])
              yelim <- c(yelim,y[i],y[j])
              zelim <- c(zelim,z[i],z[j])
            }
        }
    } 
  invisible(return(list(X=X,Y=Y,XX=XX,YY=YY,Z=Z,ZZ=ZZ,xelim=xelim, yelim=yelim, zelim=zelim, PL=PL)))
}

############################################################################

eval.alpha <- function(x, y, xygrid, vario.param, alphalim = c(0.00005,0.05), alphapval = 0.05, ni = 100,allgrid=F)
{
  nx <- length(xygrid$x)
  ny <- length(xygrid$y) 
  
  X2<-array(0,dim=c(nx,ny,ni))
  Ldet<-array(0,dim=c(nx,ny,ni))
  Z <- array(0,dim=c(length(x),ni))

  if (allgrid==T) xygrid$mask <- matrix(T,nrow=nx,ncol=ny)
  
  for(i in 1:ni)
    {
      sim <- sim.data(x = x, y = y, xygrid = xygrid, vario.param = vario.param)
      X2[,,i] <- sim$X2
      Ldet[,,i] <- sim$Ldet
      Z[,i] <- sim$z
    }
  
  nsim <- alphapval*ni
  nsim <- round(nsim,0)

  alpha0 <- alphalim[1]
  alpha1 <- alphalim[2]

  n0 <- nb.zac(X2, Ldet, xygrid = xygrid, alpha = alpha0, alphapval = alphapval)
  n1 <- nb.zac(X2, Ldet, xygrid = xygrid, alpha = alpha1, alphapval = alphapval)

  if ((n0 < nsim) == T && (n1 > nsim) == T)
    { 
      repeat
        {
          alpha <- (alpha0 + alpha1)/2
          n <- nb.zac(X2, Ldet, xygrid = xygrid, alpha = alpha, alphapval = alphapval)
          if (n > nsim)
            { 
              alpha1 <- alpha
            }
          if (n < nsim)
            { 
              alpha0 <- alpha
            }
          if (n == nsim)
            {
              ALPHA <- NULL
              N <- NULL 
              alpha5 <- alpha
              alpha1 <- alpha
              for(i in 1:4) 
                {
                  alpha <- (alpha0 + alpha1)/2
                  n <- nb.zac(X2, Ldet, xygrid = xygrid, alpha = alpha, alphapval = alphapval)
                  N <- c(N,n)
                  ALPHA <- c(ALPHA,alpha)
                  
                  if (n < nsim)
                    { 
                      alpha0 <- alpha
                    }
                  if (n == nsim)
                    { 
                      alpha1 <- alpha
                    }
                }
              l <- length(ALPHA)
              if (sum(N != nsim) == l)
                alpha <- alpha5
              else
                {
                  alpha <- min(ALPHA[N==nsim])
                }
              break
            }
        }
    }

  if (n0 >= nsim) stop("alpha0 should be smaller")
  
  if (n1 <= nsim) stop("alpha1 should be greater")
 
  invisible(return(alpha))
}

############################################################################

intersection <- function(x1,y1,x2,y2,x3,y3,x4,y4)
{
 res <- F 

  x21 <- (x2-x1)
  y21 <- (y2-y1)

  x43 <- (x4-x3)
  y43 <- (y4-y3) 

  D <- y21*x43-y43*x21

  if (D != 0) 
    {
      x31 <- (x3-x1)
      y31 <- (y3-y1)
      T1 <- (y31*x43-y43*x31)/D
      T2 <- (y31*x21-y21*x31)/D
      if ( (T1 < 0 || T1 > 1) || (T2 < 0 || T2 > 1)) {}
      else
        {
          res <- T
        } 
    }

  invisible(return(res)) 
}

############################################################################

make.grid <- function(nx,ny,poly=NULL)
  {
    if ((nx < 2) || (ny < 2)) stop("the grid must be at least of size 2x2")

    if (is.null(poly))
      {
        poly <- cbind(x[chull(x,y)],y[chull(x,y)])
      }
    
    xrang <- range(poly[, 1], na.rm = TRUE)
    yrang <- range(poly[, 2], na.rm = TRUE)
    xmin <- xrang[1]
    xmax <- xrang[2]
    ymin <- yrang[1]
    ymax <- yrang[2]
    
    xinc <- (xmax-xmin)/nx
    yinc <- (ymax-ymin)/ny
    
    xc <- xmin-xinc/2
    yc <- ymin-yinc/2
    xgrid <- rep(0,nx)
    ygrid <- rep(0,ny)
    xgrid[1] <- xc + xinc
    ygrid[1] <- yc + yinc

    for (i in 2:nx)
      {
        xgrid[i] <- xgrid[i-1]+xinc
      }
    for (i in 2:ny)
      {
        ygrid[i] <- ygrid[i-1]+yinc
      }

    yy <- matrix(xgrid,nx,ny)
    xx <- t(yy)
    yy <- matrix(ygrid,nx,ny)

    X <- as.vector(xx)
    Y <- as.vector(yy)

    poly <- rbind(poly,poly[1,])
    pts <- inpip(pts=cbind(X,Y),poly)

    X[pts] <- TRUE
    X[X!=TRUE] <- FALSE
    mask <- matrix(X,ncol=ny,nrow=nx,byrow=T)

    invisible(return(list(x=xgrid,y=ygrid,xx=xx,yy=yy,pts=pts,xinc=xinc,yinc=yinc,mask=matrix(as.logical(mask),nx,ny))))
  }

############################################################################

na2zero <- function(M)
{
	n <- nrow(M)
	m <- ncol(M)
	for(i in 1:n) {
		for(j in 1:m) {
			if(is.na(M[i, j]))
				M[i, j] <- 0
		}
	}
	invisible(return(M))
}

############################################################################

nb.zac <- function(X2, Ldet, xygrid, alpha = 0.01, alphapval = 0.05)
{
  nx <- nrow(X2)
  ny <- ncol(X2)

  bx <- xygrid$xinc
  by <- xygrid$yinc

  ni <-  dim(X2)[[3]]
  
  Nsim <- 0
  for(j in 1:ni)
    {     
      N <- 0
      T2 <- na2zero(X2[,,j])
      ldet <- Ldet[,,j]
      
      cc <- cc(T2, qchisq(1-alpha, 2))
      ncc <- cc$ncc
      if (ncc != 0)
        {
          cclab <- matrix(cc$lab, ncol = ny)
          tmax <- rep(0, ncc)
          detlambda <- rep(0, ncc)
          smaxpix <- rep(0, ncc)
          smax <- rep(0, ncc)
          prob <- rep(0, ncc)
          for(nc in 1:ncc)
            {
              smaxpix[nc] <- sum(cclab == nc)
              smax[nc] <- bx * by * (smaxpix[nc])
              tmax[nc] <- max(T2[cclab == nc])
              detlambda[nc] <- ldet[(cclab == nc) & (T2 >= tmax[nc] - 1e-07)]
              prob[nc] <- exp(( - qchisq(1-alpha, 2) * smax[nc] * sqrt(detlambda[nc]))/6.2832)
              
              if(prob[nc] <= alphapval && sum(is.na(prob[nc] <= alphapval)) == 0)
                {
                  N <- N + 1
                }
            }
          if(N > 0)
            {
              Nsim <- Nsim + 1
            }
        }
    }

  invisible(return(Nsim))
}

############################################################################

plot.zac <- function(reszac, alpha0 = NULL, ...)
{
  X <- reszac$xygrid$x
  Y <- reszac$xygrid$y
  poly <- reszac$poly
  
  zac <- reszac$zac
  if (is.null(zac)) stop("There is no ZAC")

  if (is.null(alpha0))
    {
      image(X,Y,zac,zlim=c(0,2),...)
      polygon(poly,lwd=2)
    }
  else
    {
      zac.cclab <- reszac$cclab
      zac.pvalue <- reszac$p.value
      zac.ncc <- length(zac.pvalue)
      alphapval <- reszac$alphapval
      alpha <- reszac$alpha
      if (alpha0 <= alpha) stop("alpha0 must be greater than alpha")

      X2 <- reszac$X2
      x2 <- as.array(na2zero(X2))

      cc <- cc(x2, qchisq(1-alpha0,2))
      ncc <- cc[[1]]
      CClab <- matrix(cc[[2]],ncol=length(Y))

      l <- sum(CClab!=0)
      num <- NULL
      Xi <- matrix(0,ncol=ncc,nrow=l)
      Xj <- matrix(0,ncol=ncc,nrow=l)
  
      for (nc in 1: ncc)
        {
          N <- which(CClab==nc,arr.ind=T)
          Xi[,nc] <- c(N[,1],rep(0,(l-length(N[,1]))))
          Xj[,nc] <- c(N[,2],rep(0,(l-length(N[,2]))))
        }
  
      for (mc in 1:zac.ncc)
        {
          M <- which(zac.cclab==mc,arr.ind=T)
          Ci <- M[,1]
          Cj <- M[,2]

          for (nc in 1:ncc)
            {
              for(i in 1:length(unique(M[,1])))
                {
                  if (sum(Xi[,nc] == unique(Ci)[i]) != 0)
                    {
                      cj <- unique(Cj[Ci==unique(Ci)[i]])
                      xj <- unique(Xj[Xi[,nc]==unique(Ci)[i],nc])
                      
                      for(j in 1:length(cj))
                        {
                          for(k in 1:length(xj))
                            {
                              if ((cj[j]==xj[k]) && (zac.pvalue[mc] <= alphapval)) num <- c(num,nc)
                            }
                        }
                    }
                }
            }
        }
      num <- unique(num)
  
      Ci <- NULL
      Cj <- NULL

      l <- length(num)
      for(i in 1:l)
        {
          Ci <- c(Ci,Xi[Xi[,num[i]]!=0,num[i]])
          Cj <- c(Cj,Xj[Xj[,num[i]]!=0,num[i]])
    }

      zac0 <- zac
      M <- which(!(is.na(zac)),arr.ind=T)
      for(i in 1:length(M[,1]))
        { 
          zac0[M[i,1],M[i,2]] <- 0
        }

      for(i in 1:length(Ci))
        {
          zac0[Ci[i],Cj[i]] <- 2
        }
      image(X,Y,zac0,zlim=c(0,2),...)
      polygon(poly,lwd=2)
    }  
}

############################################################################

 rmultnorm <-function(n, mu, vmat, tol = 1e-07)
{
  p <- ncol(vmat)
  if(length(mu) != p)
    stop("mu vector is the wrong length")

  M<-summary(as.vector(vmat-t(vmat)))
  if(M[6] > tol)
  stop("vmat not symmetric")

 vs <- La.svd(vmat) #,method = c("dgesdd", "dgesvd"))
 vsqrt <- t(t(vs$vt) %*% (t(vs$u) * sqrt(vs$d)))
  ans <- matrix(rnorm(n * p), nrow = n) %*% vsqrt
  ans <- sweep(ans, 2, mu, "+")
  
  dimnames(ans) <- list(NULL, dimnames(vmat)[[2]])
  
  return(ans)
}

############################################################################

sim.data <- function(x, y, xygrid, vario.param)
{
  nx <- length(xygrid$x)
  ny <- length(xygrid$y)
  
  typmod <- vario.param$typmod
  range <- vario.param$range
  sill <- vario.param$sill +  vario.param$nugget
  pep <- vario.param$nugget	

  np <- length(x)
  covmat <-  buildcovmat(x, y, typmod = typmod, ax = range, pal = vario.param$sill, pep = pep)
  invcovmat <- solve(covmat)
  z <- as.vector(rmultnorm(1, rep(0, np), covmat))
  
  S <- array(0, dim = c(nx, ny, 4))
  rho <- array(0, dim = c(nx, ny))
  s11 <- array(0, dim = c(nx, ny))
  s22 <- array(0, dim = c(nx, ny))
  ds1 <- array(0, dim = c(nx, ny, 2))
  ds2 <- array(0, dim = c(nx, ny, 2))
  dc12 <- array(0, dim = c(nx, ny, 2))
  drho <- array(0, dim = c(nx, ny, 2))
  L2.11 <- array(0, dim = c(nx, ny))
  L2.12 <- L2.11
  L2.22 <- L2.11
  L1.11 <- L2.11
  L1.12 <- L2.11
  L1.22 <- L2.11
  w1 <- array(0, dim = c(nx, ny))
  w2 <- array(0, dim = c(nx, ny))
  zk <- array(0, dim = c(nx, ny))
  X2 <- array(0, dim = c(nx, ny))
  U1 <- array(0, dim = c(nx, ny))
  U2 <- array(0, dim = c(nx, ny))
  L11 <- array(0, dim = c(nx, ny))
  L12 <- L11
  L22 <- L11
  Ldet <- array(0, dim = c(nx, ny))

  for(i in 1:nx)
    {
      for(j in 1:ny)
        {
          if (xygrid$mask[i,j]==T)
            {
              x0 <- xygrid$x[i]
              y0 <- xygrid$y[j]
  
              c0 <- buildcovvec(x0, y0, x, y, typmod = typmod, ax = range, pal = vario.param$sill, pep = pep)	
              dc0 <- buildcovgrad(x0, y0, x, y, typmod = typmod, ax = range, pal = vario.param$sill, pep = pep)	
              ddc0 <- buildcovhess(x0, y0, x, y, typmod = typmod, ax = range, pal = vario.param$sill, pep = pep)	
              G <- t(dc0) %*% invcovmat %*% dc0
              S[i, j,  ] <- as.vector(G)
              s11[i, j] <- S[i, j, 1]
              s22[i, j] <- S[i, j, 4]
              rho[i, j] <- S[i, j, 2]/sqrt(S[i, j, 1] * S[i, j, 4])	

              ds1[i, j, 1] <- t(ddc0[, 1]) %*% invcovmat %*% dc0[, 1]/s11[i, j]
              ds1[i, j, 2] <- t(ddc0[, 3]) %*% invcovmat %*% dc0[, 1]/s11[i, j]
              ds2[i, j, 1] <- t(ddc0[, 2]) %*% invcovmat %*% dc0[, 2]/s22[i, j]
              ds2[i, j, 2] <- t(ddc0[, 4]) %*% invcovmat %*% dc0[, 2]/s22[i, j]
              dc12[i, j, 1] <- t(ddc0[, 1]) %*% invcovmat %*% dc0[, 2] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 2]
              dc12[i, j, 2] <- t(ddc0[, 3]) %*% invcovmat %*% dc0[, 2] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 4]
              drho[i, j, 1] <- dc12[i, j, 1]/(sqrt(s11[i, j]) * sqrt(s22[i, j])) - rho[i, j] * (ds1[i, j, 1] + ds2[i, j, 1])
              drho[i, j, 2] <- dc12[i, j, 2]/(sqrt(s11[i, j]) * sqrt(s22[i, j])) - rho[i, j] * (ds1[i, j, 2] + ds2[i, j, 2])

              L2.11[i, j] <- (t(ddc0[, 2]) %*% invcovmat %*% ddc0[, 2])/s22[i, j] - (ds2[i, j, 1])^2
              L2.12[i, j] <- (t(ddc0[, 2]) %*% invcovmat %*% ddc0[, 4])/s22[i, j] - ds2[i, j, 1] %*% ds2[i, j, 2]
              L2.22[i, j] <- (t(ddc0[, 4]) %*% invcovmat %*% ddc0[, 4])/s22[i, j] - (ds2[i, j, 2])^2

              L1.11[i, j] <- rho[i, j] * (2 * ds1[i, j, 1] * (drho[i, j, 1] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 2]/(sqrt(s11[i, j]) * sqrt(s22[i, j]))) + (2 * ds2[i, j, 1] * t(ddc0[, 1]) %*% invcovmat %*% dc0[, 2])/(sqrt(s11[i, j]) * sqrt(s22[i, j])) - (2 * t(ddc0[, 1]) %*% invcovmat %*% ddc0[, 2])/(sqrt(s11[i, j]) * sqrt(s22[i, j])))
              L1.11[i, j] <- L1.11[i, j] + rho[i, j]^2 * (t(ddc0[, 2]) %*% invcovmat %*% ddc0[, 2]/s22[i, j] - ds2[i, j, 1]^2 - 2 * ds1[i, j, 1] * ds2[i, j, 1])
              L1.11[i, j] <- L1.11[i, j] + drho[i, j, 1]^2 - ds1[i, j, 1]^2 + t(ddc0[, 1]) %*% invcovmat %*% ddc0[, 1]/s11[i, j] - (2 * drho[i, j, 1] * t(ddc0[, 1]) %*% invcovmat %*% dc0[, 2])/sqrt(s11[i, j] * s22[i, j])
              L1.11[i, j] <- L1.11[i, j] * (1/(1 - rho[i, j]^2))
              L1.11[i, j] <- L1.11[i, j] - (rho[i, j]^2 * drho[i, j, 1]^2)/(1 - rho[i, j]^2)^2

              L1.12[i, j] <- rho[i, j] * (ds1[i, j, 1] * (drho[i, j, 2] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 4]/(sqrt(s11[i, j]) * sqrt(s22[i, j]))) + ds1[i, j, 2] * (drho[i, j, 1] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 2]/(sqrt(s11[i, j]) * sqrt(s22[i, j]))) + (ds2[i, j, 2] * t(ddc0[, 1]) %*% invcovmat %*% dc0[, 2])/(sqrt(s11[i, j]) * sqrt(s22[i, j])) + (ds2[i, j, 1] * t(ddc0[, 3]) %*% invcovmat %*% dc0[, 2])/(sqrt(s11[i, j]) * sqrt(s22[i, j])) - t(ddc0[, 1]) %*% invcovmat %*% ddc0[, 4]/(sqrt(s11[i, j]) * sqrt(s22[i, j])) - t(ddc0[, 3]) %*% invcovmat %*% ddc0[, 2]/(sqrt(s11[i, j]) * sqrt(s22[i, j])))
              L1.12[i, j] <- L1.12[i, j] + rho[i, j]^2 * (t(ddc0[, 2]) %*% invcovmat %*% ddc0[, 4]/s22[i, j] - ds2[i, j, 1] * ds2[i, j, 2] - ds1[i, j, 2] * ds2[i, j, 1] - ds1[i, j, 1] * ds2[i, j, 2])
              L1.12[i, j] <- L1.12[i, j] + drho[i, j, 1] * drho[i, j, 2] - ds1[i, j, 1] * ds1[i, j, 2] + t(ddc0[, 1]) %*% invcovmat %*% ddc0[, 3]/s11[i, j] - (drho[i, j, 2] * t(ddc0[, 1]) %*% invcovmat %*% dc0[, 2])/sqrt(s11[i, j] * s22[i, j]) - (drho[i, j, 1] * t(ddc0[, 3]) %*% invcovmat %*% dc0[, 2])/sqrt(s11[i, j] * s22[i, j])
              L1.12[i, j] <- L1.12[i, j] * (1/(1 - rho[i, j]^2))
              L1.12[i, j] <- L1.12[i, j] - (rho[i, j]^2 * drho[i, j, 1] * drho[i, j, 2])/(1 - rho[i, j]^2)^2
              L1.22[i, j] <- rho[i, j] * (2 * ds1[i, j, 2] * (drho[i, j, 2] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 4]/(sqrt(s11[i, j]) * sqrt(s22[i, j]))) + (2 * ds2[i, j, 2] * t(ddc0[, 3]) %*% invcovmat %*% dc0[, 2])/(sqrt(s11[i, j]) * sqrt(s22[i, j])) - (2 * t(ddc0[, 3]) %*% invcovmat %*% ddc0[, 4])/(sqrt(s11[i, j]) * sqrt(s22[i, j])))
              L1.22[i, j] <- L1.22[i, j] + rho[i, j]^2 * (t(ddc0[, 4]) %*% invcovmat %*% ddc0[, 4]/s22[i, j] - ds2[i, j, 2]^2 - 2 * ds1[i, j, 2] * ds2[i, j, 2])
              L1.22[i, j] <- L1.22[i, j] + drho[i, j, 2]^2 - ds1[i, j, 2]^2 + t(ddc0[, 3]) %*% invcovmat %*% ddc0[, 3]/s11[i, j] - (2 * drho[i, j, 2] * t(ddc0[, 3]) %*% invcovmat %*% dc0[, 2])/sqrt(s11[i, j] * s22[i, j])
              L1.22[i, j] <- L1.22[i, j] * (1/(1 - rho[i, j]^2))
              L1.22[i, j] <- L1.22[i, j] - (rho[i, j]^2 * drho[i, j, 2]^2)/(1 - rho[i, j]^2)^2
              w1[i, j] <- t(dc0[, 1]) %*% invcovmat %*% z
              w2[i, j] <- t(dc0[, 2]) %*% invcovmat %*% z
              zk[i, j] <- t(c0) %*% invcovmat %*% z
            }
        }
    }

  U1 <- w1/sqrt(S[,  , 1] * (1 - rho^2)) - (rho * w2)/sqrt(S[,  , 4] * (1 - rho^2))
  U2 <- w2/sqrt(S[,  , 4])
  X2 <- U1 * U1 + U2 * U2
  L11 <- (L1.11 * U1^2 + L2.11 * U2^2)/X2
  L12 <- (L1.12 * U1^2 + L2.12 * U2^2)/X2
  L22 <- (L1.22 * U1^2 + L2.22 * U2^2)/X2
  Ldet <- L11 * L22 - L12^2

  invisible(return(list(X2=X2,Ldet=Ldet,x=x,y=y,z=z)))
}

############################################################################

vario.zac <- function(reszac, dist.seq = seq(10, 200, by = 10), only.sign = T)
{
  x <- reszac$x
  y <- reszac$y
  z <- reszac$z
  cclab <- reszac$cclab
  p.value <- reszac$p.value
  alphapval <- reszac$alphapval
  
  if (is.null(cclab)) stop("There is no ZAC")

  if (only.sign==T && sum(p.value < alphapval)==0) stop("There is no significant ZAC")
  
  elim <- elim(reszac, only.sign = only.sign)

  X <- elim$X
  Y <- elim$Y
  Z <- elim$Z
  XX <- elim$XX
  YY <- elim$YY
  ZZ <- elim$ZZ
  
  xelim <- elim$xelim
  yelim <- elim$yelim

  np <- length(X)
  ncelim <- (length(x)*(length(x)-1)/2) - np

  h <- NULL
  gamma <- NULL
  for(i in 1:np)
    {
      dx <- (X[i]-XX[i])^2
      dy <- (Y[i]-YY[i])^2
      dd <- sqrt(dx + dy)
      h <- c(h,dd)
      g <- 0.5 * (Z[i]-ZZ[i])^2
      gamma <- c(gamma,g)
    }

  nlag <- length(dist.seq)
  g <- rep(0, nlag)
  d <- rep(0, nlag)
  n <- rep(0, nlag)
  
  b <- h < dist.seq[1]
  n[1] <- sum(b)
  d[1] <- sum(h[b])
  g[1] <- sum(gamma[b])

  if(n[1] > 0)
    {
      d[1] <- d[1]/n[1]
      g[1] <- g[1]/n[1]
    }
  for(ilag in 2:nlag)
    {
      b <- (h < dist.seq[ilag]) & (h > dist.seq[ilag - 1])
      n[ilag] <- sum(b)
      d[ilag] <- sum(h[b])
      g[ilag] <- sum(gamma[b])
      if(n[ilag] > 0)
        {
          d[ilag] <- d[ilag]/n[ilag]
          g[ilag] <- g[ilag]/n[ilag]
        }
    }

  vario <- list(n=n[n!=0], u=d[n!=0], v=g[n!=0], ncelim=ncelim, max.dist=max(dist.seq), trend = "cte", beta.ols = NULL, lambda=NULL)
  class(vario) <- "variogram"
  return(vario)
}

############################################################################

deteczac <- function(x, y, z, vario.param, nx=50, ny=50, poly=NULL, alpha = 0.01, alphapval = 0.05, kriging = 2)
{
  if (length(x) != length(y)) stop("x and y must have the same length")
  if (length(x) != length(z)) stop("z must have the same length as x and y")

  if(is.null(poly))
    {
      poly <- cbind(x[chull(x,y)],y[chull(x,y)])
    }
  
  xygrid <- make.grid(nx=nx,ny=ny,poly=poly)
  bx <- xygrid$xinc
  by <- xygrid$yinc
  
  typmod <- vario.param$typmod
  range <- vario.param$range
  sill <- vario.param$sill +  vario.param$nugget
  pep <- vario.param$nugget
  
  covmat <- buildcovmat(x, y, typmod = typmod, ax = range, pal = vario.param$sill, pep = pep)	
  invcovmat <- solve(covmat)
  
  np <-  length(x)
  un <- matrix(rep(1,np),ncol=1)
  den <- as.vector((t(un)%*%invcovmat%*%un))
  K <- invcovmat - invcovmat%*%un%*%t(un)%*%invcovmat/den

  z.old <- z
  if (kriging != 1) 
    {
#      z <- z/sqrt(sill)
      invcovmat <- K 
      kriging <- "ordinary"
    }
  if (kriging == 1)
    {
      z <- (z-mean(z))
#      z <- (z-mean(z))/sqrt(sill)
      kriging <- "simple"
    }

  S <- array(0, dim = c(nx, ny, 4))
  rho <- array(0, dim = c(nx, ny))
  s11 <- array(0, dim = c(nx, ny))
  s22 <- array(0, dim = c(nx, ny))
  w1 <- array(0, dim = c(nx, ny))
  w2 <- array(0, dim = c(nx, ny))
  zk <- array(0, dim = c(nx, ny))
  X2 <- array(0, dim = c(nx, ny))
  U1 <- array(0, dim = c(nx, ny))
  U2 <- array(0, dim = c(nx, ny))
  
  for(i in 1:nx)
    {
      for(j in 1:ny)
        {
          if (xygrid$mask[i,j]==T)
            {
              x0 <- xygrid$x[i]
              y0 <- xygrid$y[j]
          
              c0 <- buildcovvec(x0, y0, x, y, typmod = typmod, ax = range, pal = vario.param$sill, pep = pep)	
              dc0 <- buildcovgrad(x0, y0, x, y, typmod = typmod, ax = range, pal = vario.param$sill, pep = pep)	
              ddc0 <- buildcovhess(x0, y0, x, y, typmod = typmod, ax = range, pal = vario.param$sill, pep = pep)	
              G <- t(dc0) %*% invcovmat %*% dc0
              S[i, j,  ] <- as.vector(G)
              s11[i, j] <- S[i, j, 1]
              s22[i, j] <- S[i, j, 4]
              rho[i, j] <- S[i, j, 2]/sqrt(S[i, j, 1] * S[i, j, 4])
              
              w1[i, j] <- t(dc0[, 1]) %*% invcovmat %*% z
              w2[i, j] <- t(dc0[, 2]) %*% invcovmat %*% z
              if (kriging != 1)
                {
                  zk[i, j] <- t(c0) %*% solve(covmat) %*% z + ((1-t(c0)%*%solve(covmat)%*%un)/den)%*%t(un)%*%solve(covmat)%*%z
                }
              else
                {
                  zk[i, j] <- t(c0) %*% solve(covmat) %*% z
                }
            }
        }
    }
  U1 <- w1/sqrt(S[,  , 1] * (1 - rho^2)) - (rho * w2)/sqrt(S[,  , 4] * (1 - rho^2))
  U2 <- w2/sqrt(S[,  , 4])
  X2 <- U1 * U1 + U2 * U2
  
  Ci <- NULL
  Cj <- NULL
  cclab <- NULL
  zac <- NULL
  p.value <- NULL
  Xmax <- NULL
  Ymax <- NULL
  Tmax <- NULL
  Smax <- NULL
  detlambda <- NULL

  T2 <- na2zero(X2)

  cc <- cc(T2, qchisq(1-alpha, 2))
  ncc <- cc$ncc
  if (ncc != 0)
    {
      cclab <- matrix(cc$lab, ncol = ny)
      Tmax <- rep(0, ncc)
      detlambda <- rep(0, ncc)
      Smaxpix <- rep(0, ncc)
      Smax <- rep(0, ncc)
      p.value <- rep(0, ncc)
      for(nc in 1:ncc)
        {
          Smaxpix[nc] <- sum(cclab == nc)
          Smax[nc] <- bx * by * (Smaxpix[nc])
          Tmax[nc] <- max(T2[cclab == nc])
          coordmax <- which(T2==Tmax[nc],arr.ind=T)
          Ci <- c(Ci,coordmax[1])
          Cj <- c(Cj,coordmax[2])

          ds1 <- rep(0, 2)
          ds2 <- rep(0, 2)
          dc12 <- rep(0, 2)
          drho <- rep(0, 2)
          x0 <- xygrid$x[coordmax[1]]
          y0 <- xygrid$y[coordmax[2]]
          Xmax <- c(Xmax,x0)
          Ymax <- c(Ymax,y0)

          c0 <- buildcovvec(x0, y0, x, y, typmod = typmod, ax = range, pep = pep, pal = vario.param$sill)
          dc0 <- buildcovgrad(x0, y0, x, y, typmod =  typmod, ax = range, pep = pep, pal = vario.param$sill)
          ddc0 <- buildcovhess(x0, y0, x, y, typmod = typmod, ax = range, pep = pep, pal = vario.param$sill)
          G <- t(dc0) %*% invcovmat %*% dc0
          S <- as.vector(G)
          s11 <- S[1]
          s22 <- S[4]
          rho <- S[2]/sqrt(S[1] * S[4])
          
          ds1[1] <- t(ddc0[, 1]) %*% invcovmat %*% dc0[, 1]/s11
          ds1[2] <- t(ddc0[, 3]) %*% invcovmat %*% dc0[, 1]/s11
          ds2[1] <- t(ddc0[, 2]) %*% invcovmat %*% dc0[, 2]/s22
          ds2[2] <- t(ddc0[, 4]) %*% invcovmat %*% dc0[, 2]/s22
          dc12[1] <- t(ddc0[, 1]) %*% invcovmat %*% dc0[, 2] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 2]
          dc12[2] <- t(ddc0[, 3]) %*% invcovmat %*% dc0[, 2] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 4]
          drho[1] <- dc12[1]/(sqrt(s11) * sqrt(s22)) - rho * (ds1[1] + ds2[1])
          drho[2] <- dc12[2]/(sqrt(s11) * sqrt(s22)) - rho * (ds1[2] + ds2[2])	

          L2.11 <- (t(ddc0[, 2]) %*% invcovmat %*% ddc0[, 2])/s22 - (ds2[1])^2
          L2.12 <- (t(ddc0[, 2]) %*% invcovmat %*% ddc0[, 4])/s22 - ds2[1] %*% ds2[2]
          L2.22 <- (t(ddc0[, 4]) %*% invcovmat %*% ddc0[, 4])/s22 - (ds2[2])^2	
          
          L1.11 <- rho * (2 * ds1[1] * (drho[1] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 2]/(sqrt(s11) * sqrt(s22))) + (2 * ds2[1] * t(ddc0[, 1]) %*% invcovmat %*% dc0[, 2])/(sqrt(s11) * sqrt(s22)) - (2 * t(ddc0[, 1]) %*% invcovmat %*% ddc0[, 2])/(sqrt(s11) * sqrt(s22)))
          L1.11 <- L1.11 + rho^2 * (t(ddc0[, 2]) %*% invcovmat %*% ddc0[, 2]/s22 - ds2[1]^2 - 2 * ds1[1] * ds2[1])
          L1.11 <- L1.11 + drho[1]^2 - ds1[1]^2 + t(ddc0[, 1]) %*% invcovmat %*% ddc0[, 1]/s11 - (2 *  drho[1] * t(ddc0[, 1]) %*% invcovmat %*% dc0[,2])/sqrt(s11 * s22)
          L1.11 <- L1.11 * (1/(1 - rho^2))
          L1.11 <- L1.11 - (rho^2 * drho[1]^2)/(1 - rho^2)^2
          
          L1.12 <- rho * (ds1[1] * (drho[2] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 4]/(sqrt(s11) * sqrt(s22))) + ds1[2] * (drho[1] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 2]/(sqrt(s11) * sqrt(s22))) + (ds2[2] * t(ddc0[, 1]) %*% invcovmat %*% dc0[, 2])/(sqrt(s11) * sqrt(s22)) + (ds2[1] * t(ddc0[, 3]) %*% invcovmat %*% dc0[, 2])/(sqrt(s11) * sqrt(s22)) - t(ddc0[, 1]) %*% invcovmat %*% ddc0[, 4]/(sqrt(s11) * sqrt(s22)) - t(ddc0[, 3]) %*% invcovmat %*% ddc0[, 2]/(sqrt(s11) * sqrt(s22)))
          L1.12 <- L1.12 + rho^2 * (t(ddc0[, 2]) %*% invcovmat %*% ddc0[, 4]/s22 - ds2[1] * ds2[2] - ds1[2] * ds2[1] - ds1[1] * ds2[2])
          L1.12 <- L1.12 + drho[1] * drho[2] - ds1[1] * ds1[2] + t(ddc0[, 1]) %*% invcovmat %*% ddc0[,3]/s11 - (drho[2] * t(ddc0[, 1]) %*% invcovmat %*% dc0[, 2])/sqrt(s11 * s22) - (drho[1] * t(ddc0[, 3]) %*% invcovmat %*% dc0[,2])/sqrt(s11 * s22)
          L1.12 <- L1.12 * (1/(1 - rho^2))
          L1.12 <- L1.12 - (rho^2 * drho[1] * drho[2])/(1 - rho^2)^2

          L1.22 <- rho * (2 * ds1[2] * (drho[2] + t(dc0[, 1]) %*% invcovmat %*% ddc0[, 4]/(sqrt(s11) * sqrt(s22))) + (2 * ds2[2] * t(ddc0[, 3]) %*% invcovmat %*% dc0[, 2])/(sqrt(s11) * sqrt(s22)) - (2 * t(ddc0[, 3]) %*% invcovmat %*% ddc0[, 4])/(sqrt(s11) * sqrt(s22)))
          L1.22 <- L1.22 + rho^2 * (t(ddc0[, 4]) %*% invcovmat %*% ddc0[, 4]/s22 - ds2[2]^2 - 2 * ds1[2] * ds2[2])
          L1.22 <- L1.22 + drho[2]^2 - ds1[2]^2 + t(ddc0[,3]) %*% invcovmat %*% ddc0[, 3]/s11 - (2 * drho[2] * t(ddc0[, 3]) %*% invcovmat %*% dc0[,2])/sqrt(s11 * s22)
          L1.22 <- L1.22 * (1/(1 - rho^2))
          L1.22 <- L1.22 - (rho^2 * drho[2]^2)/(1 - rho^2)^2
              
          L11 <- (L1.11 * U1[coordmax[1], coordmax[2]]^2 + L2.11 * U2[coordmax[1], coordmax[2]]^2)/X2[coordmax[1], coordmax[2]]
          L12 <- (L1.12 * U1[coordmax[1], coordmax[2]]^2 + L2.12 * U2[coordmax[1], coordmax[2]]^2)/X2[coordmax[1], coordmax[2]]
          L22 <- (L1.22 * U1[coordmax[1], coordmax[2]]^2 + L2.22 * U2[coordmax[1], coordmax[2]]^2)/X2[coordmax[1], coordmax[2]]
          detlambda[nc] <- L11 * L22 - L12^2

          p.value[nc] <- exp((- qchisq(1-alpha, 2) * Smax[nc] * sqrt(detlambda[nc]))/6.2832)
        }

      N <- which(p.value < alphapval, arr.ind=T)
      signb <- rep(F, ncc)
      signb[N] <- T 
      zac <- X2
      
      for(i in 1:nx)
        {
          for(j in 1:ny)
            {
              if (!(is.na(zac[i,j])))
                {
                  if (cclab[i, j] == 0) zac[i, j] <- 0
                  else { zac[i, j] <- signb[cclab[i, j]] + 1 }
                }
            }
        }
    }

    invisible(return(list(zac=zac,X2=X2,zk=zk,detlambda=detlambda,p.value=p.value,cclab=cclab,xygrid=xygrid,Xmax=Xmax,Ymax=Ymax,Tmax=Tmax,Smax=Smax,x=x,y=y,z=z.old,vario.param=vario.param,alpha=alpha,alphapval=alphapval,poly=poly)))
}

############################################################################
##
##   Example
##
############################################################################

Poly <- matrix(c(-0.03373950,0.85876125,0.09335064,0.95309688,0.19502274,0.93594494,0.36532352,1.02742192,0.63983822,1.03599788,0.78217917,0.82159873,0.78726277,0.74441503,0.94993814,0.54716781,0.96010535,0.45283219,1.02619222,0.37850715,1.02619222,0.23271572,0.88385127,0.19555320,0.69067427,0.05833774,0.50003907,-0.01884595,0.41615958,-0.01884595,0.30940387,0.15839068,0.29669485,0.22128110,0.21535717,0.33276866,0.19502274,0.46426681,0.15181210,0.53287454,-0.03373950,0.79872949),ncol=2,byrow=T)

ExData <- data.frame(cbind(x=c(0.989387614,0.958398982,0.908775750,0.904678972,0.822680721,0.824101052,0.752239288,0.729133409,0.688450551,0.649470215,0.616788044,0.904405707,0.886905761,0.823256474,0.793600954,0.775143508,0.735815113,0.678437319,0.675585790,0.578629726,0.562208935,0.518658718,0.791976339,0.697609864,0.612873658,0.567973320,0.551470660,0.534848737,0.462404382,0.434985808,0.405912987,0.706177866,0.690205830,0.660276079,0.655820368,0.558055184,0.551608301,0.483603636,0.449180089,0.408128546,0.423961190,0.366070178,0.321788165,0.622850235,0.640157746,0.575379158,0.505157974,0.479397771,0.437952947,0.428227817,0.372117232,0.333894311,0.306205708,0.282100766,0.220077104,0.540939901,0.506546105,0.490978641,0.375033434,0.296583597,0.279687751,0.253559201,0.230994890,0.173548972,0.109420025,0.084788276,0.471446609,0.431603809,0.375117475,0.350659685,0.277097037,0.248813626,0.241692239,0.226259083,0.150498806,0.098012727,0.077832860,0.035371879,-0.002252952),y=c(0.27340028,0.39116496,0.41959470,0.53063274,0.56345333,0.62971853,0.70359586,0.80161311,0.86296083,0.93407874,1.01725310,0.25804114,0.31752633,0.40069766,0.46066455,0.53595709,0.60970899,0.67339696,0.71518960,0.79361786,0.89808914,0.96781694,0.18179489,0.31756438,0.51655359,0.60823373,0.66933675,0.78744125,0.86290821,0.89722246,0.98363218,0.11822387,0.21360591,0.26460425,0.34027458,0.40640354,0.50459790,0.57092884,0.60486767,0.69756296,0.79647329,0.84084031,0.93006004,0.06982767,0.19178458,0.23449142,0.28855526,0.34041481,0.44251153,0.52110956,0.57889074,0.68016198,0.74172667,0.83876599,0.89712398,0.03302984,0.10789601,0.21602212,0.34050364,0.47099897,0.54638264,0.62116702,0.71796638,0.75896351,0.84848549,0.90502517,0.01652016,0.06454687,0.11141228,0.24862178,0.28947084,0.35624128,0.41214631,0.50566493,0.58845720,0.68800935,0.68722286,0.80947797,0.84083877),z=c(332,337,135,315,336,259,252,231,252,255,226,307,329,321,132,143,269,253,261,224,240,261,352,322,365,316,234,263,234,239,241,332,307,314,320,322,353,342,269,253,251,253,271,335,311,287,326,326,260,312,449,242,262,238,256,345,327,315,311,205,142,351,214,219,238,205,259,257,280,381,286,228,226,287,161,256,257,265,268)))