% bumps.m generates a field with bumps centered at Poisson points

m = 30;                                      % avg number of points
s = 0;
N = -1;
while (s<m)
  s = s - log(rand);
  N = N + 1;                                  % N ~ Poisson(m)
end

a = rand(N,1);                                % Poisson process points in R^2
b = rand(N,1);
c = randn(N,1);                               % bump amplitude

gridsize = [30 30];                           % number of grid points in x, y
gridrange= [0 1 0 1];                         % default [xmin xmax ymin ymax]
[x,y]    = easygrid(gridrange, gridsize);     % matrices of x and y values
Z        = zeros(gridsize);

for i=1:gridsize(1),
  for j=1:gridsize(2),
    s = 0;
    for k=1:N,
      d2 = (x(i,j)-a(k))^2 + (y(i,j)-b(k))^2;
       s  = s + c(k)*exp(-100*d2);               % Gaussian bumps
%       s  = s + abs(c(k))*max(1-10*sqrt(d2),0); % positive, linear to zero
%       s  = s + c(k)*max(1-10*sqrt(d2),0);      % linear down to zero
    end
    Z(i,j) = s;
  end
end
 
clf
mesh(x,y,Z)
rotate3d on

pause

pcolor(x,y,Z)
shading flat