Você encontra a explicação desse script clicando aqui.

Você pode fazer o download desse script e dos arquivos de parâmetros (simetrica49.csv e simetrica49_varpar.csv) clicando aqui.

def.par=par(no.readonly=TRUE)
varPares=function(x,y,n,dir,tol,largBanda){
  vd=as.vector(c(sin(dir),cos(dir)))
  altCone=largBanda/tan(tol)
  largBanda=largBanda^2
  library(pracma)
  k=0
  pares=matrix(c(rep(0,n*(n-1))),
  nrow=(n*(n-1)/2),ncol=2)
  if (tol < pi/2) {
    for (i in 1:(n-1)){
      for (j in (i+1):n){
        vp=as.vector(c((x[i]-x[j]),(y[i]-y[j])))
        dist=abs(dot(vd,vp))  #produto escalar
        lb=dot(vp,vp)-dist*dist
        if (dist > altCone){
          if (lb < largBanda){
            k=k+1;pares[k,1]=i;pares[k,2]=j}
        } else{
        lbCone=(dist*tan(tol))^2
        if (lb < lbCone){
          k=k+1;pares[k,1]=i;pares[k,2]=j}
        }
      }
    }
  } else {
    for (i in 1:(n-1)){
      for (j in (i+1):n){
        k=k+1;pares[k,1]=i;pares[k,2]=j}
    }  
  }
  npares=k
  outPares=pares[1:npares,]
  return(outPares)
}
print("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
setwd("C:\\GK21\\dados\\simetrica49")
dados=read.csv("simetrica49.csv",sep=";",header=T)
x=dados$X; y=dados$Y; z=dados$Zgauss
n=nrow(dados)
deg2rad=pi/180
#parametros do variograma
param=read.csv("simetrica49_varpar.csv",sep=";",header=T)
ndir=nrow(param)
dir=param$dir*deg2rad; tol=param$tol*deg2rad; largBanda=param$largBanda
passo=param$passo; tolPasso=param$tolPasso; npassos=param$npassos
for (k in 1:ndir){
  osPares=varPares(x,y,n,dir[k],tol[k],largBanda[k])
  print(osPares)
}

Você pode fazer o download desse script e dos arquivos de parâmetros (simetrica49.csv e simetrica49_varpar.csv) clicando aqui.