function Ylm = sh(l, m, style);
%  function Ylm = sh(l, m);
%  Draw the equipotential surface of a multipole potential based on
%  Re Ylm, a spherical harmonic, degree l, order m
%  Also return numerical values evenly spaced in angle of Re Ylm.

%  if l < 0 draw instead the surface with r = Ylm(theta,phi)
%  style values
%  l = 0 - Contours
%  l = 1 - zonal, sectional, tesseral harmonics
%  l = 2 - 3D representation


nw =  6;  %  Number of samples per half wavelength
N = max([30, l*nw]);
theta = [0 : N]' *pi/N;
phi   =[-N : N]  *pi/N;

if (nargin < 3) style=2; end

%  Create the surface of Re Ylm(theta,phi)/r^(l+1) = 1;
%  Solves for r on a theta, phi grid
all = legendre(abs(l), cos(theta));
if (l == 0) all=all'; end         % Compensate for error in legendre 
Ylm = all(m+1, :)' * cos(m*phi);

%
if (style == 2)
  gamma = 1;
  if (l > 0) gamma = 1/(l+1); end
  r = abs(Ylm).^2; %.^ gamma;
else
  r=1;
end

%  Convert to Cartesian coordinates
X =  r.* (sin(theta)*cos(phi)) ;
Y =  r.* (sin(theta)*sin(phi));
Z =  r.* (cos(theta)*ones(size(phi)));

%  Color according to size and sign of Ylm
if (style == 1) 
    C = sign(Ylm);
    colormap copper 
    brighten(-0.5)
    material shiny 
else
    C = abs(Ylm);
    colormap jet
end

surf(X, Y, Z, C)
axis equal
if (style == 2) 
  ti=['Surface  Y_l^m/r^{l+1}=1  with l=' int2str(l) ', m=' int2str(m)];
else
  ti=['Contours of Re Y_l^m with l=' int2str(l) ', m=' int2str(m)];
  shading interp;
end
light
lighting phong
camlight('headlight')
axis tight equal off
%view(0,20)
view(0,10)
%camzoom(1.5)
%axis square
title(['l=' int2str(l) ', m=' int2str(m)]);