module harmonic; % Solid & spherical harmonics.
% Author: Matthew Rebbeck, ZIB.
% Advisor: Rene' Grognard.
% Date: March 1994
% Version 0 (experimental)
% Solid Harmonics of order n (Laplace polynomials)
% are homogeneous polynomials of degree n in x,y,z
% which are solutions of Laplace equation:-
% df(P,x,2) + df(P,y,2) + df(P,z,2) = 0.
% There are 2*n+1 independent such polynomials for any given n >=0
% and with:-
% w!0 = z, w!+ = i*(x-i*y)/2, w!- = i*(x+i*y)/2,
% they are given by the Fourier integral:-
% S(n,m,w!-,w!0,w!+) =
% (1/(2*pi)) *
% for u:=-pi:pi integrate (w!0 + w!+ * exp(i*u) + w!- *
% exp(-i*u))^n * exp(i*m*u) * du;
% which is obviously zero if |m| > n since then all terms in the
% expanded integrand contain the factor exp(i*k*u) with k neq 0,
% S(n,m,x,y,z) is proportional to
% r^n * Legendre(n,m,cos theta) * exp(i*phi)
% Let r2 = x^2 + y^2 + z^2 and r = sqrt(r2).
% The spherical harmonics are simply the restriction of the solid
% harmonics to the surface of the unit sphere and the set of all
% spherical harmonics {n >=0; -n <= m =< n} form a complete orthogonal
% basis on it, i.e. <n,m|n',m'> = Kronecker_delta(n,n') *
% Kronecker_delta(m,m') using <...|...> to designate the scalar product
% of functions over the spherical surface.
% The coefficients of the solid harmonics are normalised in what
% follows to yield an ortho-normal system of spherical harmonics.
% Given their polynomial nature, there are many recursions formulae
% for the solid harmonics and any recursion valid for Legendre functions
% can be 'translated' into solid harmonics. However the direct proof is
% usually far simpler using Laplace's definition.
% It is also clear that all differentiations of solid harmonics are
% trivial, qua polynomials.
% Some substantial reduction in the symbolic form would occur if one
% maintained throughout the recursions the symbol r2 (r cannot occur
% as it is not rational in x,y,z). Formally the solid harmonics appear
% in this guise as more compact polynomials in (x,y,z,r2).
% Only two recursions are needed:-
% (i) along the diagonal (n,n);
% (ii) along a line of constant n: (m,m),(m+1,m),...,(n,m).
% Numerically these recursions are stable.
% For m < 0 one has:-
% S(n,m,x,y,z) = (-1)^m * S(n,-m,x,-y,z);
algebraic procedure SolidHarmonicY(n,m,x,y,z,r2);
begin scalar mp, v, Y0, Y1, Y2;
if not (fixp(n) and fixp(m)) then return
rederr " SolidHarmonicY : n and m must be integers";
if (n < 0) then return 0;
mp := abs(m);
if (n < mp ) then return 0;
Y0 := 1/sqrt(4*Pi);
if (n = 0) then return Y0;
if (mp > 0) then
<< if m > 0 then v:=x+i*y else v:=x-i*y;
for k:=1:mp do Y0 := - sqrt((2*k+1)/(2*k))*v*Y0;
if (n > mp) then <<
k := mp + 1;
Y1 := Y0;
Y0 := z*sqrt(2*k+1)*Y1;
if (n > mp + 1) then for k:=mp+2:n do <<
Y2 := Y1;
Y1 := Y0;
Y0 := z*sqrt((4*k*k-1)/(k*k-mp*mp))*Y1
-r2*sqrt(((2*k+1)*(k-mp-1)*(k+mp-1))/
((2*k-3)*(k*k-mp*mp)))*Y2 >>;
>>;
>> else << Y1 := Y0;
Y0 := z*sqrt(3)*Y1;
if n > 1 then for k:=2:n do <<
Y2 := Y1;
Y1 := Y0;
Y0 := ( z*sqrt(4*k*k-1)*Y1 - r2*(k-1)*
sqrt((2*k+1)/(2*k-3))*Y2)/k >>;
>>;
if m < 0 and not evenp mp then Y0 := - Y0;
return Y0
end;
algebraic procedure SphericalHarmonicY(n,m,theta,phi);
SolidHarmonicY(n,m,sin(theta)*cos(phi),
sin(theta)*sin(phi),cos(theta),1)$
endmodule;
end;