%
%  w14spher2.m -- (djm:  06 apr 2004)
%    - spherical harmonics

ll = 3;

%  colormap
Nc = 20;
s = transpose([0:1/20:1 ones(1,20)]);  fs = flipud(s);
red2blu = [s min(fs,s) fs].^2;

%  spherical coordinates
Ns = 64;
[X,Y,Z]   = sphere(Ns);
[Ph,Th,R] = cart2sph(X,Y,Z);

%  quadrature tests
trap = ones(size(X));
trap(1,:) = 1/2;  trap(end,:) = 1/2;  trap(:,end) = 0;
trap = trap .* sin(Th+(pi/2)) * (2*pi/Ns) * (pi/Ns);
normcheck = [];

for mm = 0:ll
    figure(mm+1);  clf
    Plms = legendre(ll,cos(Th+(pi/2)),'norm');
    Ylm  = squeeze(Plms(abs(mm)+1,:,:)).*exp(i*mm*Ph)/sqrt(2*pi);

    surf(X,Y,Z,real(Ylm))
    axis equal
    caxis([-1 1]*max(abs(Ylm(:))));
    shading interp;  colormap(red2blu);  colorbar
    light
    
    title(['\bf m=' num2str(mm) ' mode'])
    
    normcheck = [normcheck sum(sum( (abs(Ylm).^2) .* trap ))];
end

normcheck
ortho = 2*ones(ll+1,ll+1);

for mm = 0:ll
    Plms = legendre(ll,cos(Th+(pi/2)),'norm');
    Ylm  = squeeze(Plms(abs(mm)+1,:,:)).*exp(i*mm*Ph)/sqrt(2*pi);

    for jj = 0:mm
        Ylm0  = squeeze(Plms(abs(jj)+1,:,:)).*exp(i*jj*Ph)/sqrt(2*pi);
        
        ortho(mm+1,jj+1) = sum(sum( (Ylm .* conj(Ylm0) .* trap )));
        ortho(jj+1,mm+1) = ortho(mm+1,jj+1);
    end
end

ortho