%% Calculate Spherical Harmonics
% Spherical harmonics arise in the solution to Laplace's equation and are used to represent functions defined on the surface of a sphere. Use legendre to compute and visualize the spherical harmonic for Y_3^2.
% The equation for spherical harmonics includes a term for the Legendre function, as well as a complex exponential:
% Y_l^m(\theta,\phi) = \sqrt{\frac{\left(2l+1\right) \left(l-m\right)!}{4 \pi \left(l+m\right)!}} P_l^m(\cos \theta) e^{i m \phi}, \qquad -l\le m \le l.
%
% First, create a grid of values to represent all combinations of 0\le \theta \le \pi \; (colatitude angle) and 0\le \phi \le 2\pi (azimuthal angle). Here, the colatitude \theta ranges from 0 at the North Pole, to \pi/2 at the Equator, and to \pi at the South Pole.

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

%
% Calculate P_l^m \left(\cos \;\theta \;\right) on the grid for l=3.

l = 3;
Plm = legendre(l,cos(theta));

%
% Since legendre computes the answer for all values of m, Plm contains some extra function values. Extract the values for m=2 and discard the rest. Use the reshape function to orient the results as a matrix with the same size as phi and theta.

m = 2;
if l ~= 0
    Plm = reshape(Plm(m+1,:,:),size(phi));
end

%
% Calculate the spherical harmonic values for Y_3^2.

a = (2*l+1)*factorial(l-m);
b = 4*pi*factorial(l+m);
C = sqrt(a/b);
Ylm = C .*Plm .*exp(1i*m*phi);

%
% Convert the spherical coordinates to Cartesian coordinates. Here, \pi/2 - \theta becomes the latitude angle that ranges from \pi/2 at the North Pole, to 0 at the Equator, and to -\pi/2 at the South Pole. Plot the spherical harmonic for Y_3^2 using both the positive and negative real values.

[Xm,Ym,Zm] = sph2cart(phi, pi/2-theta, abs(real(Ylm)));
surf(Xm,Ym,Zm)
title('$Y_3^2$ spherical harmonic','interpreter','latex')

%
% Copyright 2018 The MathWorks, Inc.