function plothump(n,seed)
% create a sum of n elliptical Gaussian humps centered in [0,1]x[0,1]

% Paul Heckbert		9 Oct 2000

if nargin<1
    n = 3;
end
if nargin<2
    seed = 0;
end
amp_min = -1;	% amplitude range (negative to make depressions, not hills)
amp_max = -.2;
sigma_max = .4;	% standard deviation range
sigma_min = .02;

% initialize the random number generator (for repeatable results)
rand('state',seed);

% center of each Gaussian
cen = rand(n,2);

% amplitude of each Gaussian
amp = amp_min + (amp_max-amp_min)*rand(n,1);

% standard deviations of the major and minor axes of each Gaussian
% Gaussian i is circular iff sigma(i,1)==sigma(i,2)
sigma = sigma_min + (sigma_max-sigma_min)*rand(n,2);

% rotation angle (radians CCW) of each Gaussian
rot = pi*rand(n,1);

% Gaussian is amp*exp(-.5*(u^2/sigma1^2+v^2/sigma2^2))
% where [u v]' = [cos(rot) sin(rot); -sin(rot) cos(rot)] [x-cenx y-ceny]'

% evaluate function on a fine grid
xvec = (0:1/400:1)';
[xg,yg] = meshgrid(xvec,xvec);
zg = zeros(size(xg));
for i = 1:n
    x = xg-cen(i,1);
    y = yg-cen(i,2);
    c = cos(rot(i));
    s = sin(rot(i));
    ax = sigma(i,1)^-2;
    ay = sigma(i,2)^-2;
    A = -[c^2*ax+s^2*ay c*s*(ax-ay); c*s*(ax-ay) s^2*ax+c^2*ay]/2;
    zg = zg + amp(i)*exp( A(1,1)*x.^2 + 2*A(1,2)*x.*y + A(2,2)*y.^2 );
end

figure(1);
contour(xg,yg,zg,20);	% contour plot

% contour(xg,yg,zg,(-1:.05:.9));	% contour plot to see peak

% figure(2);
% contour(xg,yg,zg,-.3.^(0:50));	% contour plot to see tails

% figure(3);
% surf(xg,yg,-zg);		% filled surface
