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

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.

Answer 1

[David Goodmanson]

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.

Had they used

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]

Hi David,

Thank you for clarifying my doubt. That was really helpful.

[chaitanya acharya]

The code below helps in visualization of the negative 'm' values as well.

%% visualization of spherical harmonics

% phi

clear;
clc;

dx = pi/60;
alt = 0:dx:pi;
az = 0:dx:2*pi;
[phi,theta] = meshgrid(az,alt);

n=1;
m=-1;
leg_nm=legendreP2(n,m,cos(theta));

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))

[chaitanya acharya]

If we vary alt from 0 to 2*pi and az from 0 to pi, the results come out wierd when m is odd. Why is such a weird behaviour observed?

[David Goodmanson]

Hi chaitanya, what is legendreP2 ?

[chaitanya acharya]

Hi David,

Its a code to choose the corresponding row after using legendre function in matlab.

You can replace it by this.

leg_nm1=legendre(n,cos(theta));
leg_nm=leg_nm1(m+1,:);

[David Goodmanson]

Hi chaitanya,

Using the two lines above, the code doesn't work. Could you supply code that works and gives a funny result?

[chaitanya acharya]

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]

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                            [1]

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 [1] 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]

Hi David,

That was really helpful. I now understand the dicontinuties involved here much better.

Thank you.