%%%% XXX meshgrid!!

%% ----------------------------------------------------------------
clf;
rotate3d on;

%% ----------------------------------------------------------------
%n = 4;
n = 20;
%n = 100;

xlo = -2.0;
xhi =  2.0;
nx  =  n;
dx  = (xhi - xlo) / nx;

ylo = -2.0;
yhi =  2.0;
ny  =  n;
dy  = (yhi - ylo) / ny;

[x,y] = meshgrid(xlo:dx:xhi, ylo:dy:yhi);

%% ----------------------------------------------------------------
z = x.^2 + y.^2;

z1 = 5 * exp(-x.^2) .* exp(-y.^2);
z2 = -3 * exp(-(x+1).^2) .* exp(-(y-1).^2);
z2 = -3 * exp(-(x+1).^2) .* exp(-(y).^2);
z3 = -2 * exp(-(x+0.4).^2) .* exp(-(y+1.2).^2);
z4 = 1.5 * sin(x.^2) .* cos(y.^2);
z = z1 + z2 + z3 + z4;
%z = z1 + z2 + z4;
%z = z1 + z2;
%z = z1;

z = x.*exp(-x.^2-y.^2);

%% ----------------------------------------------------------------
hold on
contour(z);
figure;
rotate3d on;
surf(z);

%% ----------------------------------------------------------------
% The diff operators reduce one dimension by 1 but not the other.
% The gradient operator is dimension-preserving.  Only the latter is
% suitable input for a quiver plot.

[zx, zy] = gradient(z, dx, dy);

%contour(zx);
%contour(zy);

% Vector field with arrows.
figure; % Open up a second window
quiver(x, y, zx, zy);
