% random spherical harmonic, hacked from Denise L. Chen 9-1-93.
% Barnett 2/25/10. Fixed normalization to prevent overflow, 9/3/13

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

degree = 40;       % overall degree
a = 0*p;       % accumulation array
Y = legendre(degree,cos(theta),'norm'); % Calculate the harmonic, normalized
for order = 0:degree
  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 = 1; %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
end
  
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 three ways

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

figure; surf(x,y,z,'cdata', real(a)); colormap(jet(256));
light, lighting gouraud, shading interp, axis tight equal off; view(170,30)
camzoom(1.5);

r = 1+0.01*real(a);
figure; surf(x.*r,y.*r,z.*r, 1*(real(a)>0)); colormap(jet(256));
light, lighting gouraud, shading interp, axis tight equal off; view(170,30)
camzoom(1.5);
