my.empvgram<-function(ys,us,vs,nbins,half=1){
  #ys is vector of responses
  #us is vector of u-indices (horizontal coordinates)
  #vs is vector of v-indices (vertical coordinates)
  #half denotes whether or not to restrict estimates
  #to one half the maximum distance between points
  #(default=true)
  #--------------------------------------------------
  #distances
  uu<-outer(us,us,FUN="-")
  uu<-uu[upper.tri(uu)]
  vv<-outer(vs,vs,FUN="-")
  vv<-vv[upper.tri(vv)]
  ds<-sqrt(uu^2+vv^2)
  #--------------------------------------------------
  #squared value differences
  yy<-outer(ys,ys,FUN="-")
  yy<-yy[upper.tri(yy)]
  #---------------------------------------------------
  #if half=1 only use data up to 1/2 maximum distance
  if(half==1){
  cut<-max(ds)/2
  yy<-yy[ds<=cut]
  ds<-ds[ds<=cut]}
  #----------------------------------------------------
  #create  bins
  ovals<-order(ds)
  ds<-ds[ovals]
  yy<-yy[ovals]
  bwid<-max(ds)/nbins
  bins<-bwid*(0:nbins)
  #-----------------------------------------------------
  #empirical variograms (Matheron and Cressie robust)
  bcnt<-0
  gs<-NULL
  rgs<-NULL
  ns<-NULL
  mids<-NULL
  repeat{
    bcnt<-bcnt+1
    tup<-bins[bcnt+1]
    tlow<-bins[bcnt]
    tys<-yy[(ds>tlow)&(ds<=tup)]
    tn<-length(tys)
    tmid<-(tlow+tup)/2
    #Matheron
    tg<-0.5*(1/tn)*sum(tys^2)
    #Cressie and Hawkins
    tgr<-0.5*(1/((0.457*tn+0.494)*tn^3))*(sum(abs(tys)^0.5))^4
    gs<-c(gs,tg)
    rgs<-c(rgs,tgr)
    ns<-c(ns,tn)
    mids<-c(mids,tmid)
    if(bcnt==nbins) break
   }
  #------------------------------------------------------
  res<-data.frame(n=ns,d=mids,g=gs,rgs=rgs)
  return(res)
  }