'Calculate the classical and robust estimates of the emperical 
'semivariogram.
'program written by Kari Jovaag (kari@iastate.edu) July 1997.

'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","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

'get the maximum number of lags to use
numlag = MsgBox.Input("Give the maximum number of lags to use",
         "Number of Lags","20").AsNumber

'Make a list for the lags
LagList = list.Make

'Calculate the cutoff points for the lags
lag = maxdist/numlag
for each i in 0..numlag
  LagList.Add(i*lag)
end

'Make lists for the avg distances, number of pairs at each distance,
'sum of the squared differences at each distance, and sum of the
'square root of the absolute value of the differences(SRob)
DistList = list.Make
NhList = List.Make
SSDList = List.Make
SRobList = List.Make

'Set initial values to 0
for each i in 0..(numlag-1)
  DistList.Add(0)
  NhList.Add(0)
  SSDList.Add(0)
  SRobList.Add(0)
end

'calculate the number of pairs, sum of squared differences,
'sum of the square root of the absolute value of the differences,
'and sum of the distances for each lag
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)
'make sure each pair is included only once
        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)
            srob = ((theVtab.ReturnValue(theData,rec) - 
                    theVtab.ReturnValue(theData,rec2)).abs.sqrt)
            for each i in 1..numlag
              low = LagList.Get(i-1) 
              hi = LagList.Get(i)
              if((dist >= low) and (dist < hi))then
                num = i-1
              else
                continue
              end
            end
            DistList.Set(num,(DistList.Get(num) + dist))
            NhList.Set(num,(NhList.Get(num) + 1))
            SSDList.Set(num,(SSDList.Get(num) + ssd))
            SRobList.Set(num,(SRobList.Get(num) + srob))
          end
        end
      end   
    end
  end
end

'make a new table for the results
newVtab = Vtab.MakeNew("em_vgm.dbf".AsFileName,dBASE)

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

'make fields for the 2 estimates, avg distance, and # of pairs
cla = Field.Make("ClasicalSVg",#FIELD_FLOAT,12,4)
rob = Field.Make("RobustSVg",#FIELD_FLOAT,12,4)
dis = Field.Make("Distance",#FIELD_FLOAT,12,4)
nh = Field.Make("#Pairs",#FIELD_FLOAT,12,4)
newVtab.AddFields({dis,nh,cla,rob})

'calculate the estimates and put the results in the table
for each i in 0..(numlag-1)
  nhval = NhList.Get(i)
  if(nhval = 0) then
    continue
  else
    rec = newVtab.AddRecord
    avgdis = DistList.Get(i)/nhval
    claest = ((1/nhval)*(SSDList.Get(i)))/2
    robest = ((((1/nhval)*(SRobList.Get(i)))^4)/
                                (0.457 + (0.494/nhval)))/2
    newVtab.SetValue(dis,rec,avgdis)
    newVtab.SetValue(nh,rec,nhval)
    newVtab.SetValue(cla,rec,claest)
    newVtab.SetValue(rob,rec,robest)
  end
end

'make and open a table document to display the results
newtable = Table.Make(newVtab)
newTable.GetWin.Open