9 views (last 30 days)

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.

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.

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 [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.

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.