spherical harmonics

spherical_harmonics

function Ymn = spherical_harmonics(degree, order, theta, phi)

% Compute the spherical harmonics of given degree and order at theta,phi
%
% Parameters
% ----------
% degree: int, the degree of the underlying Legendre polynomial
% order: int, the order of the r
%
% theta: bounded between 0 and pi
% phi: bounded between 0 and 2pi
%
% Returns
% -------
% Ymn: the value of this spherical harmonic in the points given by theta,
%      phi
%
% Notes
% -----
% Based on: http://en.wikipedia.org/wiki/Spherical_harmonics
% With the Condon-Shortley phase normalization scheme (http://mathworld.wolfram.com/Condon-ShortleyPhase.html)

size_theta = size(theta);
size_phi = size(phi);

if size_theta ~= size_phi
    error('Inputs must have the same size');
end

% Linearize
theta = theta(:);
phi = phi(:);

Lmn = legendre_lm(degree,order,cos(theta));
% The
Ymn = Lmn' .* exp(1i.*order.*phi);

% The normalization factor:
N = sqrt(((2*degree+1) * factorial(degree-order)) / (factorial(degree+order) * (4 * pi)));

Ymn = Ymn * N;

% Return something that conforms to the shape of the original inputs:
Ymn = reshape(Ymn, size_theta);

end

% Helper functions:
function L = legendre_lm(l,m,x)
% Function to get the legendre polynomial of x with degree l and order m:

L = legendre(l, x); 

% By convention this is the symmetry across the 0 point:
if m<0
    m = -1 * m;
    % By convention:
    L = (-1) ^ m .* (factorial(l-m)/factorial(l+m)) .* L;
    
end

L = L(m+1,:);

end