% covariogram2.m generates field values, computes a sample covariogram, and compares to the actual covariogram

clear

meth = 4;

switch meth
% --------------------------- Set up x and y values on a grid
case 1,
gridsize = [20 20];                           % number of grid points in x, y
gridrange= [0 5 0 5];                       % default [xmin xmax ymin ymax]
[x,y]    = easygrid(gridrange, gridsize);     % matrices of x and y values
x        = reshape(x,prod(gridsize),1);
y        = reshape(y,prod(gridsize),1);

% --------------------------- Set up uniformly distributed x and y values
case 2,
numpoints= 400;
x        = 10*rand(numpoints,1);
y        = 10*rand(numpoints,1);

% --------------------------- Set up blocks of x and y values
case 3,
numblocks= 10;
perblock = 40;
x = [];
y = [];
for i=1:numblocks,
  a = 5*rand;                                   % lower left corner of block
  b = 5*rand;
  x = [x; a+rand(perblock,1)];                  % append new x values to old
  y = [y; b+rand(perblock,1)];
end

% --------------------------- Set up blocks of x and y values
otherwise,
perblock = 40;
x = [];
y = [];
[a, b] = easygrid([0 4 0 4],[3 3]);
for i=1:3,
 for j=1:3,
  x = [x; a(i,j)+rand(perblock,1)];             % append new x values to old
  y = [y; b(i,j)+rand(perblock,1)];
 end
end

end

% --------------------------- Generate field values
m        = 1000;                               % number of Fourier modes

A        = randn(1,m)*sqrt(2/m);                % amplitudes of Fourier modes
phi      = 2 * pi * rand(m,1);                  % random phase
norm     = 0.1*(100/0.1).^rand(m,1);            % almost scale invariant
norm     = 1 + 2*round(rand(m,1)) + rand(m,1);  % two wide bands
norm     = 1 + 2*round(rand(m,1)) + 0.1*rand(m,1);% two narrow bands
norm     = sqrt(-2*log(rand(m,1)));             % rougly exponential
norm     = 3 + rand(m,1);                       % one wide band
theta    = 2 * pi * rand(m,1);                  % arguments of wave numbers
K        = [norm.*cos(theta) norm.*sin(theta)]; % isotropic wave numbers

for i=1:length(x),
  Z(i,1) = A * cos(K*[x(i); y(i)] + phi);         % compute field value at x,y
end

figure(1)
clf

subplot(3,2,1)
hist(Z,20)
a = axis;
axis([-max(abs(Z)) max(abs(Z)) a(3) a(4)]);
title('Histogram of field values')

subplot(3,2,2)
plot(x,y,'*')
title('Locations of data points')
drawnow

n = prod(size(x));                      % the total number of entries in x
d = dist([x y], [x'; y']);              % Euclidean distances between points
D = reshape(d,n^2,1);                   % rearrange these into a vector

mZ= mean(Z);
e = (Z-mZ)*(Z-mZ)';                     % all possible products
E = reshape(e,n^2,1);                   % square and rearrange into a vector

subplot(3,1,2)
L = min(160000,length(D));               % limit the number of points to plot
plot(D(1:L), E(1:L), '.');              % plot products vs. point dist
xlabel('Distances in the plane')
ylabel('Covariance products')
drawnow

subplot(3,1,3)
a=[];
b=[];
a(1) = 0;                               % distance zero
b(1) = var(Z);                          % estimate variance at zero lag
md   = 0.9*max(D);                      % most of the maximum distance
nb   = 30;                              % maximum number of bins to use
j    = 2;                               % current bin
for i=1:nb,                             % form bins
  L = find(D > (i-1)*md/nb & D <= i*md/nb); 
                                        % indices with certain point distances
  if length(L) > 0,                     % this bin is nonempty
    a(j) = mean(D(L));                  % average these point distances 
    b(j) = trimmean(E(L),0);            % average the corresponding E values
    j    = j+1;
  end
end
plot(a,b);
hold on
xlabel('Distances in the plane')
ylabel('Covariance')
title('Actual (red), sample (blue), and normalized sample (green) covariogram')

plot(a,b/b(1),'g')

% --------------------------- Use wave number norms to estimate covariance
r = (0:120)*md/120;                     % distances at which to estimate covar
B = sum(besselj(0,norm*r))/m;           % estimate E J_0(|K|*r)

plot(r,B,'r')
grid on