% random spherical harmonic, hacked from Denise L. Chen 9-1-93.
% barnett 2/25/10; 3/13/12 sums over all sph harms up to given order.

delta = pi/1000;         % Create the grid
theta = 0 : delta : pi; % altitude
phi = 0 : 2*delta : 2*pi; % azimuth
[p,t] = meshgrid(phi,theta);

a = 0*p;       % accumulation array

tic;  % it's now slow to run.
maxdeg = 80;  % overall max degree
for degree=0:maxdeg  % each degree
Y = legendre(degree,cos(theta)); % Calculate the harmonic
for order = 0:degree
  fprintf('degree, order = (%d,%d)\n', degree, order)
  yc = repmat(Y(order+1,:)', [1 numel(phi)]) .* ...
        repmat(cos(order*phi), [size(Y,2), 1]);
  ys = repmat(Y(order+1,:)', [1 numel(phi)]) .* ...
        repmat(sin(order*phi), [size(Y,2), 1]);
  nrm = sqrt(delta^2*sum(yc(:).^2 .* sin(t(:))));  % find L^2 norm of cos term
  a = a + randn(1)/nrm  * yc + randn(1)/nrm * ys;  % sin term zero for order=0
%  a = a + 1/nrm  * yc + 1/nrm * ys;  % sin term zero for order=0
end
%delta^2*sum(a(:).^2 .* sin(t(:)))  % check norm^2 = 2p+1 when all coeffs=1
end
toc

rho = 1;
% Apply spherical coordinate equations
r = rho.*sin(t);
x = r.*cos(p);    % spherical coordinate equations
y = r.*sin(p);
z = rho.*cos(t);

% Plot the surface
figure; surf(x,y,z,'cdata', 1*(real(a)>0)); colormap(gray(2)); % nodal
%figure; surf(x,y,z,'cdata', real(a)); colormap(jet(256)); % value
light, lighting gouraud, shading interp, axis tight equal off; view(170,30)
camzoom(1.5);

%print -djpeg rpw_sphere_nodal_1000_deg40.jpg
%print -djpeg rpw_sphere_nodal_500_deg0-40.jpg
print -djpeg rpw_sphere_nodal_1000_deg0-80.jpg
