'Calculate the classical estimate of the emperical variogram

'make a list of fields that will be included
theTable = av.GetActiveDoc
theVtab = theTable.GetVtab
aFieldlist = List.Make
for each f in theVtab.GetFields
'only look for the numeric, visible fields
  if ( f.IsVisible ) then
    aFieldList = aFieldList.Add(f.GetAlias)
  end
end

'get the variable to use as the x coordinates
TheXName = MsgBox.ListAsString(aFieldList,
  "Choose the field to use as x coordinates","X")
  if (TheXName <> nil) then
    TheX = theVtab.FindField(TheXName)
  end

'get the variable to use as the y coordinates
TheYName = MsgBox.ListAsString(aFieldList,
  "Choose the field to use as y coordinates","Y")
  if (TheYName <> nil) then
    TheY = theVtab.FindField(TheYName)
  end

'get the variable to use as the data
TheDataName = MsgBox.ListAsString(aFieldList,
  "Choose the field to use as data","MP Data")
  if (TheDataName <> nil) then
    TheData = theVtab.FindField(TheDataName)
  end
  
'get the maximum distance
maxdist = MsgBox.Input("Give the maximum distance to use",
          "Maximum Distance","0").AsNumber


DistList = list.Make
NhList = List.Make
SSDList = List.Make

for each rec in theVtab
  if(theVtab.ReturnValue(theData,rec).IsNull) then
    continue
  else
    x1 = TheVTab.ReturnValue(TheX,rec)
    y1 = TheVTab.ReturnValue(TheY,rec)
    for each rec2 in theVtab
      if(theVtab.ReturnValue(theData,rec2).IsNull) then
        continue
      else
        x2 = TheVTab.ReturnValue(TheX,rec2)
        y2 = TheVTab.ReturnValue(TheY,rec2)
        if((y2<y1) or ((x2<=x1)and(y2=y1))) then
          continue
        else
          dist = (((x1-x2)^2) + ((y1-y2)^2)).sqrt
          if(dist>maxdist) then
            continue
          else
            ssd = ((theVtab.ReturnValue(theData,rec) - 
                    theVtab.ReturnValue(theData,rec2))^2)
            num = DistList.FindByValue(dist)
            if(num = -1)then
              DistList.Add(dist)
              NhList.Add(1)
              SSDList.Add(ssd)
            elseif(num > -1) then
              NhList.Set(num,(NhList.Get(num) + 1))
              SSDList.Set(num,(SSDList.Get(num) + ssd))
            end
          end
        end
      end   
    end
  end
end

newVtab = Vtab.MakeNew("em_vgm.dbf".AsFileName,dBASE)

'Make the Vtab editable
if (newVtab.IsEditable = False) then
  newVtab.StartEditingWithRecovery
end

vgm = Field.Make("VgrmEst",#FIELD_FLOAT,12,4)
dis = Field.Make("Distance",#FIELD_FLOAT,12,4)
nh = Field.Make("#Pairs",#FIELD_FLOAT,12,4)

newVtab.AddFields({dis,nh,vgm})

numrec = (DistList.Count - 1)
for each i in 0..numrec
  rec = newVtab.AddRecord
  vgmest = ((1/NhList.Get(i))*(SSDList.Get(i)))
  newVtab.SetValue(dis,rec,DistList.Get(i))
  newVtab.SetValue(nh,rec,NhList.Get(i))
  newVtab.SetValue(vgm,rec,vgmest)
end

newtable = Table.Make(newVtab)