# In the matlab page of Associated Legendre Polynomials, is the example for the spherical harmonics correct?

9 views (last 30 days)
chaitanya acharya on 9 Dec 2019
Commented: chaitanya acharya on 17 Dec 2019
The '+ve' nodes will always be adjacent to the '-ve' modes in any vibration problem. In the example shown, upper 4 nodes are in one phase while the lower ones are in the opposite phase.

David Goodmanson on 9 Dec 2019
Edited: David Goodmanson on 10 Dec 2019
Hi chaitanya,
Good catch. Really good catch.
In addition to the signs of the lobes, I believe that Matlab's Y32 is functionally incorrect.
On Matlab's 'doc legendre' page there are a couple of issues. For spherical coordinates with Ylm's, most of the world uses
theta = elevation angle heading down from the z axis, 0 =< theta =< pi
Matlab is using
theta = elevation angle from the x-y plane, -pi/2 <= theta <= pi/2
Compared to the usual way of doing things, their choice of elevation angle has the effect of interchanging sin(theta) and cos(theta). Since they plug their cos(theta) into the associated legendre function, it appears that their result is not functionally correct.
They use
[Xm,Ym,Zm] = sph2cart(phi, theta, real(Y32))
The sp2cart function uses elevation from the x-y plane, which is appropriate for their theta. [In the sph2cart function itself the definition of phi and theta is swapped, which doesn't help]. But since real(Y32) can be negative, half the time that puts the point [x,y,z] on the opposite side of the origin from what you might think. That will have an effect on the signs of the lobes.
The code below doesn't use sph2cart and uses abs(real(Y32)) for the radius. The lobes do end up having alternating signs as a function of azimuthal angle phi.
delta = pi/60;
alt = 0:delta:pi;
az = 0:delta:2*pi;
[theta,phi] = meshgrid(alt,az);
l = 3;
Plm = legendre(l,cos(theta));
m = 2;
P32 = reshape(Plm(m+1,:,:), size(theta));
a = (2*l+1)*factorial(l-m);
b = 4*pi*factorial(l+m);
C = sqrt(a/b);
Y32 = C .* P32 .* exp(1i*m*phi);
R = abs(real(Y32));
Xm = R.*sin(theta).*cos(phi);
Ym = R.*sin(theta).*sin(phi);
Zm = R.*cos(theta);
surf(Xm,Ym,Zm,(real(Y32)))
title('\$Y_3^2\$ spherical harmonic','interpreter','latex')
xlabel('x')
ylabel('y')
absYsq = Y32.*conj(Y32);
normint = trapz(trapz(sin(theta).*absYsq))*delta^2 % should be 1
normint =
0.999999998921089
As a check, the normalization integral should equal 1, which it does.
If you ignore the signs of the lobes, do some appropriate rotations with the mouse and compare the shapes of Matlab's lobes vs. those in the code above, they're different. Also the normalization integral for Matlab, which includes a factor of cos(theta) instead of sin(theta), comes out wrong.
Plm = legendre(l,sin(theta));
they would have gotten the correct result for Ylm, although their method of plotting might still cause some problems.

chaitanya acharya on 13 Dec 2019
Hi David,
Sorry. I made a mistake. Here is the code. It gives funny results when m is odd.
%% visualization of spherical harmonics
% phi
clear;
clc;
dx = pi/60;
alt = 0:dx:2*pi;
az = 0:dx:pi;
[phi,theta] = meshgrid(az,alt);
n=1;
m=1; % m>0 && m<=n
leg_nm1=legendre(n,cos(theta));
leg_nm(:,:)=leg_nm1(m+1,:,:);
Ynm=sqrt(((2*n+1)*factorial(n-m))/(4*pi*factorial(n+m))).*leg_nm.*exp(1i*m*phi);
%% cartesian to spherical co-ordinates
%[Xm,Ym,Zm] = sph2cart(phi,theta,abs(Ynm));
R =abs(real(Ynm));
Xm = R.*sin(theta).*cos(phi);
Ym = R.*sin(theta).*sin(phi);
Zm = R.*cos(theta);
surf(Xm,Ym,Zm,real(Ynm))
David Goodmanson on 14 Dec 2019
Hi chaitanya,
Well, yes, since you are intentionally changing the anglular domain from
0<=theta<=pi
0<=phi<=2pi
to
0<=theta<=2pi
0<=phi<=pi 
there are going to be problems that really have nothing to do with how to calculate the Ylms. Both domains cover all possible points on a sphere, but the second one is not a continous map like the first one is. Consider a path going around the earth at constant latitude, say 45 degrees. Looking down at the path from above the north pole, you see a circle projected onto the xy plane. Let's say that the +x axis is to the right and the +y axis is up. As phi goes from 0 to pi, you go aound the top of the circle from the +x to the -x axis. Around the top, y has only positive values.
To go around the bottom half you need negative y values, but  above restricts sin(phi) to positive values. Since we have
Xm = R.*sin(theta).*cos(phi);
Ym = R.*sin(theta).*sin(phi);
Zm = R.*cos(theta);
the only way to make it happen is
sin(theta) picks up a minus sign
cos(theta) remains the same (still in the northern hemisphere)
This means that as you cross the -x axis to go around the bottom half of the circle, theta abruptly jumps from 45 degrees to 315 degrees. Theta remains there as phi jumps from pi down to to 0 and then goes from 0 to pi to trace out the bottom half of the circle.
The discontinouties cause problems, as you found. If you add
xlabel('x')
ylabel('y')
to the end of the code, you can see the discontinuities occur as you go from positive to negative y. With the original anglular domain there are no such problems.
chaitanya acharya on 17 Dec 2019
Hi David,
That was really helpful. I now understand the dicontinuties involved here much better.
Thank you.