% covariogram2.m uses (x,y,Z) data to produce a sample covariogram

s = size(x);                            % number of rows and columns in x
n = prod(s);                            % the total number of entries in x

X = reshape(x,n,1);                     % put all x values into a single vector
Y = reshape(y,n,1);                     % put all y values into a single vector
F = reshape(Z,n,1);                     % put all field values into a vector
F = F - mean(F);                        % subtract the mean

d = dist([X Y], [X'; Y']);              % Euclidean distances between points
D = reshape(d,n^2,1);                   % rearrange these into a vector

e = F*F';                               % covariance data
E = reshape(e,n^2,1);                   % rearrange these into a vector, square

figure(2)                               % work in a second figure window
clf

if n < 1000,
subplot(2,1,1)
plot(D, E, '.');                        % plot data distances vs. point dist
xlabel('Distances in the plane')
ylabel('Covariance data')

subplot(2,1,2)
end

r = range(D);
mm = min(D);
for i=1:40,
  L = find(D >= mm+(i-1)*r/40 & D < mm+i*r/40); 
                                        % indices with certain point distances
  a(i) = mean(D(L));                    % average these point distances 
  b(i) = mean(E(L));                    % average E values
end
a=a(find(a*0==0));                      % leave out missing values (NaN)
b=b(find(b*0==0));
plot(a,b);
hold on
xlabel('Distances in the plane')
ylabel('Sample covariance at given distance')