function [xx,yy,z] = planet_map(num_pts_x_and_y)
% By Neil E. Cotter, 3 July 2002.
% Create map of a planet.
% Input:
%  num_pts_x_and_y = number of points to generate in both
%     x and y directions.  Gives num_pts_x_and_y^2 points.
% Outputs:
%  x = x-coordinates (longitude)
%  y = y-coordinates (latitude)
%  z = z-coordinates (altitude)


% Create the x,y pts
spacing = 10/(num_pts_x_and_y - 1);
x = -5 : spacing : 5.01;
y = -5 : spacing : 5.01;
[xx,yy] = meshgrid(x,y);

% Sum up gaussian-shaped bumps to create desired surface.
% Pool located at (2,2).
z = -exp(-(((2-xx).^2)/4 + (2-yy).^2));

% Pit located at (3.5, -3.5).
z = z - exp(-((3.5-xx).^2 + (-3.5-yy).^2)/(1.5^2));

% Mesa located at (1,-1).
z = z + exp(-((1-xx).^2 + ((-1-yy).^2)/4));

% Drop off = edge of circle centered at (-3.5,13).
z = z - (sqrt((-3.5-xx).^2 + (13-yy).^2) < 9.5);

% Sand dune located at (-2.5,1.5)
z = z + exp(-(((-2.5-xx).^2)/4 + ((1.5-yy).^2)/4));

% Set spire locations (x,y) and width(squared).
spire_params_matrix ...
  = [-2, -1, 1/8; -2, -4, 1/8; -4, -3, 1/8; -1, -2, 1/8];
for spire_params = spire_params_matrix'
  z = z + exp(-(((spire_params(1)-xx).^2)/spire_params(3) ...
                + ((spire_params(2)-yy).^2)/spire_params(3)));
end
return
