function     [grd,gammah,S2,cnt]=variogram(x,y,var,res,range,anisotropy_para)     
% function [grd,gammah,S2,cnt]=variogram(x,y,var,res,range,anisotropy_para)
% compute normalized semi-variogram model parameter
%% INPUT
%    var=f(x,y)		 (x,y) coordinates of variable var
%    res   = lag resolution
%    range = range limit for semivariogram computation
%    anisotropy_para = [angle ratio]
%            angle = orientation angle from horizontal axis in degree
%            ratio = aspect ratio of the scale in vertical (Y) to that of horizontal (X)
%% OUTPUT
%    grd =  lag grid
%    gammah = gamma(h) 	or gamma(grd)   semivariogram
%    S2  =  variance at each lag
%    cnt =  no. of pairs at each lag
%

%%  Kriging Software Package  version 1.0 beta,   October 29, 1998
%%  Copyright (c) 1998, property of Dezhang Chu and Woods Hole Oceanographic
%%  Institution.  All Rights Reserved.

ni=length(var);				% number of data points of observations
indx_nan=find(isnan(var));
var(indx_nan)=[];
x(indx_nan)=[];
y(indx_nan)=[];

nnan=length(indx_nan); 
n=ni-nnan;				% usable number of data points      
disp(sprintf('usable data points = %g',n))

if n-nnan < 1 
   disp('No usable data, press any key to cont.');
   pause
   grd=nan;
   gammah=nan;
   S2=nan;
   cnt=nan;
   return
end

% dealing with anisotropy
ang=anisotropy_para(1);
aspect=anisotropy_para(2);
if abs(ang) > eps & abs(1-aspect) > 1e-5
%% observations
    cz=[x(:) y(:)];
    rx=1;
    ry=aspect;
    [cz,rot]=trans(cz,ang,rx,ry);
    x=cz(:,1);y=cz(:,2);
end

m=floor(range/res);
x=reshape(x,n,1);
y=reshape(y,n,1);
var=reshape(var,n,1);

cnt=zeros(m,1);
gamma_sum=zeros(m,1);
var_sum=zeros(m,1);
var2_sum=zeros(m,1);
gammah=zeros(m,1);
for i=1:n
   if rem(i,200) == 0 disp(i);end
   d=sqrt( (x(i)-x).^2+(y(i)-y).^2);
   [dsort,indx_sort]=sort(d);
   dsort=dsort(2:n);								% excluding
   indx_sort=indx_sort(2:n);
   var_sort=var(indx_sort);       
   dindx=floor(dsort/res)+1;
   k_acc=0;								% accumulated k index
   for k=1:m
     indxk=find(dindx(k_acc+1:n-1) == k) + k_acc ;
     nk=length(indxk);
     if nk > 0
        cnt(k)=cnt(k)+nk;
        k_acc=k_acc+nk;
        var_dif= var(i) - var_sort(indxk);	
        var_sum(k)=var_sum(k)+sum(var_sort(indxk));
        var2_sum(k)=var2_sum(k)+var_sort(indxk)'*var_sort(indxk);
        gamma_sum(k)=gamma_sum(k)+ var_dif'*var_dif;
      end
   end
end
indx=find( cnt ~= 0);
if length(indx) > 0
   cnt_nz=cnt(indx);
   var_mean=var_sum(indx)./cnt_nz;
   c0=var2_sum(indx)./cnt_nz-var_mean.^2;
   gammah(indx)=0.5*gamma_sum(indx)./(cnt_nz.*c0);
   S2(indx)=c0./gammah(indx);
else
   gammah=nan.*ones(m,1);
   S2=nan*ones(m,1);
end
grd=[0:res:res*m]';	% grid for plotting
gammah=[0;gammah];
cnt=[mean(cnt);cnt(:)];
