function z = hump1(x,y)
% a sum of n 2-D elliptical Gaussian humps centered in [0,1]x[0,1]
% evaluated at point (x,y)
% x and y scalars, returns a scalar

% Paul Heckbert		9 Oct 2000

n = 3;		% number of humps
seed = 4;	% for random number generator

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]'


% now evaluate the function
z = 0;
for i = 1:n
    c = cos(rot(i));
    s = sin(rot(i));
    uv = [c s; -s c]*[x-cen(i,1); y-cen(i,2)];	% translated & rotated (x,y)
    z = z + amp(i)*exp(-.5*sum((uv./sigma(i,:)').^2));
end