∑MathAdvanced

Spherical Harmonics & SO(3) Representations

Key Points

•
Spherical harmonics are smooth wave patterns on the sphere that form an orthonormal basis, much like sine and cosine form a basis on the circle.
•
The rotation group SO(3) acts on spherical harmonics by mixing only the 2l+1 functions at the same degree l, which is the essence of irreducible representations.
•
Each band l is a (2l+1)-dimensional irreducible representation of SO(3) governed by Wigner D-matrices.
•
You can expand any square-integrable function on the sphere as a linear combination of spherical harmonics and rotate it by blockwise multiplying with $D^{l}$ .
•
In C++, you can evaluate spherical harmonics using associated Legendre polynomials and numerically stable normalization via gamma functions.
•
Rotations of spherical harmonic coefficients use Euler angles and Wigner D-matrices, which factor into complex exponentials and a real d-matrix.
•
For practical tasks like graphics lighting or geophysics, low-band expansions (e.g., $l \leq 3$ to 8) are often enough and computationally efficient.
•
Be consistent with conventions: physics (Condon–Shortley phase), angle order (θ is polar, φ is azimuth), and real vs complex bases.

Prerequisites

→Linear Algebra (vectors, matrices, eigenvalues) — Representations act on vector spaces, and rotations/orthogonality require understanding inner products and matrix multiplication.
→Calculus on the Sphere — Integration with the sin θ Jacobian and understanding eigenfunctions of differential operators appear in definitions and proofs.
→Complex Numbers and Euler's Formula — Spherical harmonics use e^{i m φ} and complex conjugation identities in projections and rotations.
→Special Functions (Legendre Polynomials) — Associated Legendre polynomials are the core building blocks of spherical harmonics.
→3D Rotations and Euler Angles — Wigner D-matrices are parameterized by Euler angles; understanding rotation parametrizations avoids convention errors.
→Numerical Methods (stability, recurrence relations) — Stable computation requires lgamma, recurrences, and careful summation to prevent overflow and loss of precision.
→Group Theory Basics — Irreducible representations and group actions on function spaces motivate blockwise rotation and the structure of D^l.

Detailed Explanation

Tap terms for definitions

01Overview

Spherical harmonics are special functions defined on the surface of a sphere that act like the spherical analog of sines and cosines. They provide a complete, orthonormal basis for representing square-integrable functions on the sphere, meaning any reasonable function on a sphere can be written as a weighted sum of these basis functions. The rotation group SO(3) describes all possible 3D rotations and has a deep relationship with spherical harmonics: each fixed degree l forms a (2l+1)-dimensional irreducible representation under rotations. This structure allows us to rotate spherical data (like environment maps, gravitational fields, or quantum wavefunctions) by applying small, well-structured linear operations within each degree l block rather than manipulating the function pointwise. Mathematically, spherical harmonics are eigenfunctions of the Laplace–Beltrami operator on the sphere with eigenvalues -l(l+1). In applications, we truncate the infinite expansion at some maximum degree L to obtain a band-limited approximation. Computation involves associated Legendre polynomials, careful normalization, and, for rotations, Wigner D-matrices parameterized by Euler angles. The combination of completeness, orthogonality, and clean rotation behavior makes spherical harmonics a powerful tool across physics, computer graphics, and signal processing on the sphere.

02Intuition & Analogies

Imagine painting patterns on a basketball. Some patterns have no wiggles (a flat color), some have a single north–south gradient, and others have more and more stripes that wrap around and between the poles. These patterns are the spherical equivalents of the familiar sine waves on a string. Spherical harmonics are precisely those building-block patterns: smooth ripples that fit perfectly on a sphere without edges or seams. Each pattern is labeled by two integers l and m. The number l counts how many bands of variation you see from pole to pole, while m counts how many waves wrap around the equator. As l increases, the patterns gain finer detail, just as higher frequencies do on a guitar string. Now consider rotating the basketball. If you rotate the sphere, you don’t create new kinds of ripples—you only rotate the existing pattern. In fact, all patterns with the same l behave like a tightly knit family: under rotation, they only mix among themselves. That family is an irreducible representation of the rotation group SO(3). The Wigner D-matrix is like a recipe that tells you how to remix the family members (the different m values) when you rotate the ball by certain angles. This viewpoint is powerful: instead of redrawing the painted pattern after a rotation, you just update the coefficients in the recipe by multiplying with a small matrix for each l. The result is faster, more stable, and perfectly faithful to rotation.

03Formal Definition

For integers l

\geq

0 and m

\in

[-l, l], the (complex) spherical harmonics

Y_{l}^{m}

(

θ

ϕ

) are defined by

Y_{l}^{m}

(

θ

ϕ

) =

N_{l m}

P_{l}^{m}

(

cos

θ

)

e^{im ϕ}

, where

θ

\in

[0,

π

] is the polar angle from the positive z-axis and

ϕ

\in

[0, 2

π

) is the azimuth. Here

P_{l}^{m}

denotes the associated Legendre function with the Condon–Shortley phase, and

N_{l m}

\frac{2 l + 1}{4 π} \frac{( l - m )!}{( l + m )!}

is the normalization constant ensuring orthonormality. The set {

Y_{l}^{m}

} over all l and m forms a complete orthonormal basis for

L^{2}

(

S^{2}

), so any f

\in

L^{2}

(

S^{2}

) can be expanded as f(

θ

ϕ

) =

\sum_{l = 0}^{\infty}

\sum_{m = - l}^{l}

a_{l m}

Y_{l}^{m}

(

θ

ϕ

), with coefficients

a_{l m}

\int_{S^{2}}

θ

ϕ

)

\overline{Y_{l}^{m} (θ, ϕ)}

\; d

Ω

. The irreducible unitary representations of SO(3) are indexed by l and act on the (2l+1)-dimensional space spanned by {

Y_{l}^{m}

}_{m=-l}^{l}. Under a rotation R with Euler angles (

α

β

γ

), the basis transforms as

Y_{l}^{m}

(

R^{- 1}

\hat{r}

) =

\sum_{m^{'} = - l}^{l}

D_{m m^{'}}^{l}

(

α

β

γ

)

Y_{l}^{m^{'}}

(

\hat{r}

), where

D^{l}

is the Wigner D-matrix with entries

D_{m m'}^{l} = e^{-i

α

}

d_{m m^{'}}^{l}

(

β

)

e^{- i m^{'} γ}

04When to Use

Use spherical harmonics when your data naturally lives on the sphere or when rotation invariance/equivariance matters. In computer graphics, they compress and rotate low-frequency environment lighting for real-time shading; a band limit of l \le 3 often suffices for diffuse effects. In geophysics and planetary science, gravitational and magnetic fields around planets are expanded in spherical harmonics to analyze global variations and to filter noise. In quantum mechanics, they are the angular part of the hydrogen atom’s wavefunctions and represent angular momentum eigenstates; rotation properties via Wigner D-matrices are central to coupling and selection rules. In 3D signal processing and robotics, spherical harmonics enable rotation-equivariant neural operators or descriptors for point clouds. Use the representation theory viewpoint (SO(3) irreps) when you need to rotate coefficient vectors without resampling the function. This is valuable for fast operations: rotating a band-limited field reduces to independent small matrix multiplications for each l, which is more stable and efficient than rotating samples on a mesh. If your function is smooth and dominated by low frequencies, a low-band SH expansion gives an excellent approximation with far fewer numbers than a dense sampling.

⚠️Common Mistakes

• Mixing conventions: Some sources define Y_{l}^{m} without the Condon–Shortley phase or swap \theta and \phi. Decide on one convention (commonly physics: \theta is polar, \phi is azimuth, with the phase) and stick to it across code and formulas. • Degrees vs radians: Trigonometric functions in C++ take radians; passing degrees silently corrupts results. • Real vs complex SH: Graphics often use real SH, while physics uses complex SH. Rotations are simplest with complex SH via Wigner D; for real SH, you need real-valued rotation blocks obtained by basis change. • Normalization errors: Forgetting the \sqrt{\frac{2l+1}{4\pi} \frac{(l-m)!}{(l+m)!}} factor or double-counting the (−1)^{m} phase leads to incorrect magnitudes and broken orthonormality. • Numerical overflow: Factorials grow fast; use log-gamma (lgamma) or recurrence relations to keep values in range. • Sampling weights: When projecting from samples, include the Jacobian sin \theta and use adequate quadrature in \theta; uniform sampling in \theta without weights biases polar regions. • Rotation ordering: Euler angle conventions (ZYZ vs ZXZ, active vs passive rotations) change D^{l}; match your math to your 3D engine’s convention. • Expecting rotations to mix different l: Rotations never mix different degrees; if your code does, the basis or indexing is inconsistent. • Aliasing from too low band limit L: Sharp features require higher L; truncation smears or rings (Gibbs-like artifacts). • Indexing m from 0..2l incorrectly: Remember m runs from −l to l; a common storage is m_to_index = m + l.

Key Formulas

Spherical Harmonic Definition

Y_{l}^{m} (θ, ϕ) = N_{l m} P_{l}^{m} (cos θ) e^{im ϕ}, N_{l m} = \frac{2 l + 1}{4 π} \frac{( l - m )!}{( l + m )!}

Explanation: This defines complex spherical harmonics in terms of associated Legendre functions and a normalization factor. It ensures orthonormality over the unit sphere using the physics convention.

Orthonormality

\int_{S^{2}} Y_{l}^{m} (θ, ϕ) \overline{Y_{l^{'}}^{m^{'}} (θ, ϕ)} d Ω = δ_{l l^{'}} δ_{m m^{'}}

Explanation: Different spherical harmonics are orthogonal, and each has unit energy over the sphere. This property allows clean coefficient extraction via inner products.

SH Expansion

f (θ, ϕ) = l = 0 \sum \infty m = - l \sum l a_{l m} Y_{l}^{m} (θ, ϕ), a_{l m} = \int_{S^{2}} f (θ, ϕ) \overline{Y_{l}^{m} (θ, ϕ)} d Ω

Explanation: Any square-integrable function on the sphere can be written as a sum of spherical harmonics. Coefficients are computed by projecting the function onto each basis function.

Eigenvalue Equation

\nabla_{S^{2}}^{2} Y_{l}^{m} (θ, ϕ) = - l (l + 1) Y_{l}^{m} (θ, ϕ)

Explanation: Spherical harmonics are eigenfunctions of the spherical Laplacian with eigenvalues -l(l+1). This characterizes their 'frequency' on the sphere.

Conjugation Symmetry

Y_{l}^{- m} = (- 1)^{m} \overline{Y_{l}^{m}}

Explanation: Negative-order harmonics are proportional to the complex conjugates of positive-order ones. This identity is useful to reduce computation and storage.

Wigner D Factorization

D_{m m^{'}}^{l} (α, β, γ) = e^{- im α} d_{m m^{'}}^{l} (β) e^{- i m^{'} γ}

Explanation: The full rotation matrix at degree l factors into simple complex exponentials and the real Wigner small-d matrix. This makes implementation and analysis easier.

Wigner small-d (closed form)

d_{m m^{'}}^{l} (β) = k = k_{m i n} \sum k_{m a x} (- 1)^{k - m + m^{'}} \frac{( l + m )! ( l - m )! ( l + m ^{'} )! ( l - m ^{'} )!}{( l + m - k )! k ! ( m ^{'} - m + k )! ( l - m ^{'} - k )!} (cos \frac{β}{2})^{2 l + m - m^{'} - 2 k} (sin \frac{β}{2})^{m^{'} - m + 2 k}

Explanation: This series computes the small-d elements with finite summation limits $k_{min}$ = $max$ (0,m-m') and $k_{ma x}$ = $min$ (l+m, l-m'). It is practical for moderate l with log-gamma for stability.

Associated Legendre Recurrences

P_{m}^{m} (x) = (- 1)^{m} (2 m - 1)!! (1 - x^{2})^{m /2}, P_{m + 1}^{m} (x) = x (2 m + 1) P_{m}^{m} (x), P_{l}^{m} (x) = \frac{( 2 l - 1 ) x P _{l - 1}^{m} ( x ) - ( l + m - 1 ) P _{l - 2}^{m} ( x )}{l - m}

Explanation: These recurrences compute $P_{l}^{m}$ (x) efficiently without large factorials. They are numerically robust for typical degrees used in applications.

Addition Theorem

m = - l \sum l Y_{l}^{m} (Ω) \overline{Y_{l}^{m} (Ω^{'})} = \frac{2 l + 1}{4 π} P_{l} (cos γ)

Explanation: The sum over orders at fixed l depends only on the angle $γ$ between $Ω$ and $Ω$ '. This underlies rotationally invariant kernels and fast algorithms.

Parseval/Plancherel on S^2

l = 0 \sum \infty m = - l \sum l ∣ a_{l m} ∣^{2} = \int_{S^{2}} ∣ f (Ω) ∣^{2} d Ω

Explanation: Energy is preserved between spatial domain and SH coefficients. It helps verify implementations and choose truncation levels.

Complexity Analysis

Let L be the band-limit (maximum degree). The number of coefficients scales as

\sum_{l = 0}^{L}

(2l+1) = (L+1)^{2}, so both storage and bandwidth for SH coefficients are O(

L^{2}

). Evaluating a single

Y_{l}^{m}

(

θ

ϕ

) using associated Legendre recurrences costs O(l), dominated by building

P_{l}^{m}

; computing all harmonics up to L at one direction costs O(

L^{2}

) if you reuse intermediate

P_{l}^{m}

. Projecting a sampled function onto SH via numerical quadrature with N samples is O(N

L^{2}

) time and O(

L^{2}

) space, since for each sample you update all coefficients. More specialized quadrature (e.g., Gauss–Legendre in

θ

with FFT in

ϕ

) can reduce practical constants and enable O(

L^{2}

lo g

L) transforms, but the straightforward approach is O(N

L^{2}

). Rotations act block-diagonally in l. Building a full Wigner

D^{l}

by the closed-form small-d series for all m,m' takes O(

l^{3}

) time because there are O(

l^{2}

) entries each with an O(l) summation. Applying

D^{l}

to a vector of size (2l+1) takes O(

l^{2}

). If you must rotate all bands up to L and you rebuild

D^{l}

each time, the worst-case total cost is

\sum_{l = 0}^{L}

l^{3}

) = O(

L^{4}

). In many applications, L is small (e.g., L

\leq

8), making this tractable. Precomputing and caching

D^{l}

for common angles avoids recomputation; for varying angles, use stable recurrences or factorization to reduce overhead. Numerically, using log-gamma to form factorial-based coefficients keeps values in range and maintains accuracy for moderate L (typically up to 50–100 in double precision).

Code Examples

Evaluate complex spherical harmonics Y_l^m(θ, φ) with stable normalization

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 // Utility: constants
5 static const double PI = 3.141592653589793238462643383279502884;
6 
7 // Compute associated Legendre P_l^m(x) for l >= 0, |m| <= l, with Condon–Shortley phase
8 // Uses stable recurrences. Returns P_l^m(x) for m >= 0. For m < 0, use relation outside.
9 static double associated_legendre(int l, int m, double x) {
10     m = abs(m);
11     if (l < m) return 0.0;
12     // P_m^m(x)
13     double pmm = 1.0;
14     if (m > 0) {
15         double somx2 = sqrt(max(0.0, 1.0 - x*x));
16         double fact = 1.0;
17         for (int i = 1; i <= m; ++i) {
18             pmm *= -fact * somx2; // includes Condon–Shortley (-1)^m and (2m-1)!! via recurrence
19             fact += 2.0;
20         }
21     }
22     if (l == m) return pmm;
23     // P_{m+1}^m(x) = x (2m+1) P_m^m(x)
24     double pmmp1 = x * (2*m + 1) * pmm;
25     if (l == m + 1) return pmmp1;
26     // Upward recurrence for l
27     double pll = 0.0;
28     for (int ll = m + 2; ll <= l; ++ll) {
29         pll = ((2*ll - 1) * x * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
30         pmm = pmmp1;
31         pmmp1 = pll;
32     }
33     return pll;
34 }
35 
36 // Normalization N_{lm} = sqrt((2l+1)/(4pi) * (l-m)!/(l+m)!) using log-gamma to avoid overflow
37 static double sh_normalization(int l, int m) {
38     int am = abs(m);
39     double lognum = log(2.0 * l + 1.0) - log(4.0 * PI) + lgamma(l - am + 1.0);
40     double logden = lgamma(l + am + 1.0);
41     return exp(0.5 * (lognum - logden));
42 }
43 
44 // Complex spherical harmonic Y_l^m(theta, phi) using physics convention (theta: polar, phi: azimuth)
45 static complex<double> spherical_harmonic(int l, int m, double theta, double phi) {
46     double x = cos(theta);
47     if (m >= 0) {
48         double Plm = associated_legendre(l, m, x);
49         double Nlm = sh_normalization(l, m);
50         complex<double> eimphi = polar(1.0, m * phi); // e^{i m phi}
51         return Nlm * Plm * eimphi;
52     } else {
53         int mp = -m;
54         // Use Y_l^{-m} = (-1)^m * conj(Y_l^{m})
55         complex<double> Ypos = spherical_harmonic(l, mp, theta, phi);
56         double phase = (mp % 2 == 0) ? 1.0 : -1.0; // (-1)^m with m=mp
57         return phase * conj(Ypos);
58     }
59 }
60 
61 int main() {
62     // Example: evaluate Y_3^m at theta=1.0 rad, phi=2.0 rad for m=-3..3
63     double theta = 1.0, phi = 2.0;
64     int l = 3;
65     cout.setf(std::ios::fixed); cout << setprecision(8);
66     for (int m = -l; m <= l; ++m) {
67         complex<double> Y = spherical_harmonic(l, m, theta, phi);
68         cout << "Y_" << l << "^" << m << "(theta,phi) = "
69              << Y.real() << " + i" << Y.imag() << "\n";
70     }
71     // Quick sanity: Y_l^{-m} ≈ (-1)^m * conj(Y_l^{m})
72     complex<double> Ypos = spherical_harmonic(l, 2, theta, phi);
73     complex<double> Yneg = spherical_harmonic(l, -2, theta, phi);
74     complex<double> rhs = ((2 % 2 == 0) ? 1.0 : -1.0) * conj(Ypos);
75     cout << "Check conjugation symmetry (m=2): |Y_-2 - (-1)^2 conj(Y_2)| = " << abs(Yneg - rhs) << "\n";
76     return 0;
77 }
78

This program evaluates complex spherical harmonics using the physics convention with the Condon–Shortley phase. It computes associated Legendre polynomials via stable recurrences, uses log-gamma for normalization to avoid overflow, and enforces the symmetry for negative m. A simple check verifies Y_l^{-m} = (-1)^m conj(Y_l^{m}).

Time: O(l) for a single Y_l^m; O(L^2) if you compute all m for l up to L at one directionSpace: O(1) additional space

Rotate SH coefficients within a fixed band l using Wigner D-matrix

1 #include <bits/stdc++.h>
2 using namespace std;
3 static const double PI = 3.141592653589793238462643383279502884;
4 
5 // Log-factorial via lgamma
6 static inline double logfact(int n) { return lgamma(n + 1.0); }
7 
8 // Compute Wigner small-d element d^l_{m m'}(beta) using the closed-form sum with log-gamma stabilization
9 static double wigner_small_d(int l, int m, int mp, double beta) {
10     // Summation limits
11     int kmin = max(0, m - mp);
12     int kmax = min(l + m, l - mp);
13     double cb2 = cos(0.5 * beta);
14     double sb2 = sin(0.5 * beta);
15 
16     // Precompute log of the big sqrt factorial term
17     double logsqrt = 0.5 * (logfact(l + m) + logfact(l - m) + logfact(l + mp) + logfact(l - mp));
18 
19     // Accumulate in double using exp of logs; Kahan summation for better accuracy
20     double sum = 0.0, c = 0.0; // compensation
21     for (int k = kmin; k <= kmax; ++k) {
22         int a = l + m - k;
23         int b = k;
24         int c1 = mp - m + k;
25         int d = l - mp - k;
26         if (a < 0 || b < 0 || c1 < 0 || d < 0) continue; // safety
27         double sign = ((k - m + mp) % 2 == 0) ? 1.0 : -1.0;
28         double logden = logfact(a) + logfact(b) + logfact(c1) + logfact(d);
29         int p = 2*l + m - mp - 2*k; // power of cos
30         int q = mp - m + 2*k;       // power of sin
31         double term = sign * exp(logsqrt - logden) * pow(cb2, p) * pow(sb2, q);
32         // Kahan summation
33         double y = term - c;
34         double t = sum + y;
35         c = (t - sum) - y;
36         sum = t;
37     }
38     return sum;
39 }
40 
41 // Build D^l matrix for given Euler angles (alpha, beta, gamma) in ZYZ convention
42 static vector<vector<complex<double>>> wigner_D_matrix(int l, double alpha, double beta, double gamma) {
43     int n = 2*l + 1;
44     vector<vector<complex<double>>> D(n, vector<complex<double>>(n));
45     for (int m = -l; m <= l; ++m) {
46         for (int mp = -l; mp <= l; ++mp) {
47             double dl = wigner_small_d(l, m, mp, beta);
48             complex<double> left = polar(1.0, -m * alpha);
49             complex<double> right = polar(1.0, -mp * gamma);
50             D[m + l][mp + l] = left * dl * right;
51         }
52     }
53     return D;
54 }
55 
56 // Multiply D^l by coefficient vector a_{l, m} (m=-l..l) to rotate coefficients
57 static vector<complex<double>> rotate_band_l(int l, const vector<complex<double>>& a, double alpha, double beta, double gamma) {
58     auto D = wigner_D_matrix(l, alpha, beta, gamma);
59     int n = 2*l + 1;
60     vector<complex<double>> out(n, {0.0, 0.0});
61     for (int i = 0; i < n; ++i) {
62         complex<double> acc = 0.0;
63         for (int j = 0; j < n; ++j) acc += D[i][j] * a[j];
64         out[i] = acc;
65     }
66     return out;
67 }
68 
69 int main() {
70     int l = 2; // example band
71     int n = 2*l + 1;
72     vector<complex<double>> a(n, {0.0, 0.0});
73     // Put unit energy in m=0 component: a_{2,0} = 1
74     a[l + 0] = 1.0;
75 
76     // Rotate by Euler angles (alpha, beta, gamma) in radians (ZYZ convention)
77     double alpha = 0.3, beta = 0.7, gamma = -0.2;
78     auto a_rot = rotate_band_l(l, a, alpha, beta, gamma);
79 
80     cout.setf(std::ios::fixed); cout << setprecision(6);
81     cout << "Rotated coefficients for l=2 (m=-2..2):\n";
82     for (int m = -l; m <= l; ++m) {
83         auto v = a_rot[m + l];
84         cout << "m=" << m << ": " << v.real() << " + i" << v.imag() << "\n";
85     }
86     return 0;
87 }
88

This code constructs the Wigner D^l matrix from Euler angles using the closed-form small-d series and applies it to rotate a vector of SH coefficients within degree l. It uses log-gamma to keep factorial-related quantities stable and the ZYZ Euler convention. Each l-block rotates independently; here we demonstrate l=2.

Time: O(l^3) to build D^l (O(l^2) entries each summing O(l)) and O(l^2) to apply it to a vectorSpace: O(l^2) to store D^l and O(l) for vectors

Project a function f(θ, φ) onto SH coefficients by numerical quadrature

1 #include <bits/stdc++.h>
2 using namespace std;
3 static const double PI = 3.141592653589793238462643383279502884;
4 
5 // From Example 1: associated Legendre and SH (copied for self-containment)
6 static double associated_legendre(int l, int m, double x) {
7     m = abs(m);
8     if (l < m) return 0.0;
9     double pmm = 1.0;
10     if (m > 0) {
11         double somx2 = sqrt(max(0.0, 1.0 - x*x));
12         double fact = 1.0;
13         for (int i = 1; i <= m; ++i) {
14             pmm *= -fact * somx2;
15             fact += 2.0;
16         }
17     }
18     if (l == m) return pmm;
19     double pmmp1 = x * (2*m + 1) * pmm;
20     if (l == m + 1) return pmmp1;
21     double pll = 0.0;
22     for (int ll = m + 2; ll <= l; ++ll) {
23         pll = ((2*ll - 1) * x * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
24         pmm = pmmp1;
25         pmmp1 = pll;
26     }
27     return pll;
28 }
29 
30 static double sh_normalization(int l, int m) {
31     int am = abs(m);
32     double lognum = log(2.0 * l + 1.0) - log(4.0 * PI) + lgamma(l - am + 1.0);
33     double logden = lgamma(l + am + 1.0);
34     return exp(0.5 * (lognum - logden));
35 }
36 
37 static complex<double> spherical_harmonic(int l, int m, double theta, double phi) {
38     double x = cos(theta);
39     if (m >= 0) {
40         double Plm = associated_legendre(l, m, x);
41         double Nlm = sh_normalization(l, m);
42         complex<double> eimphi = polar(1.0, m * phi);
43         return Nlm * Plm * eimphi;
44     } else {
45         int mp = -m;
46         complex<double> Ypos = spherical_harmonic(l, mp, theta, phi);
47         double phase = (mp % 2 == 0) ? 1.0 : -1.0;
48         return phase * conj(Ypos);
49     }
50 }
51 
52 // Example target function: f(θ, φ) = 1 + 0.5 * cos θ
53 static double f_example(double theta, double phi) {
54     (void)phi; // no dependence on phi
55     return 1.0 + 0.5 * cos(theta);
56 }
57 
58 int main() {
59     int L = 3;                 // band-limit
60     int nTheta = 100;          // polar samples
61     int nPhi = 200;            // azimuth samples
62 
63     // Coefficients a_{l m} for 0<=l<=L, -l<=m<=l stored in a flat vector
64     vector<complex<double>> coeffs; coeffs.resize((L+1)*(L+1)); // (L+1)^2 entries
65 
66     auto idx = [&](int l, int m){ return l*l + l + m; }; // maps (l,m) to [0,(L+1)^2)
67 
68     // Quadrature: uniform in phi, uniform in theta with sin(theta) weight (simple Riemann rule)
69     double dphi = 2.0 * PI / nPhi;
70     double dtheta = PI / nTheta;
71 
72     for (int it = 0; it < nTheta; ++it) {
73         // midpoint rule in theta for better accuracy
74         double theta = (it + 0.5) * dtheta;
75         double sinth = sin(theta);
76         for (int ip = 0; ip < nPhi; ++ip) {
77             double phi = (ip + 0.5) * dphi;
78             double val = f_example(theta, phi);
79             for (int l = 0; l <= L; ++l) {
80                 for (int m = -l; m <= l; ++m) {
81                     complex<double> Y = spherical_harmonic(l, m, theta, phi);
82                     coeffs[idx(l,m)] += val * conj(Y) * (sinth * dtheta * dphi);
83                 }
84             }
85         }
86     }
87 
88     cout.setf(std::ios::fixed); cout << setprecision(8);
89     // Print non-negligible coefficients
90     for (int l = 0; l <= L; ++l) {
91         for (int m = -l; m <= l; ++m) {
92             auto a = coeffs[idx(l,m)];
93             double mag = abs(a);
94             if (mag > 1e-4) {
95                 cout << "a_{" << l << "," << m << "} = " << a.real() << " + i" << a.imag() << "\n";
96             }
97         }
98     }
99 
100     // Compare with analytic values for f = 1 + 0.5 cos(theta):
101     // a_{0,0} = sqrt(4*pi); a_{1,0} = 0.5 * sqrt(4*pi/3); others zero.
102     cout << "Expected a_{0,0} ~ " << sqrt(4.0*PI) << "\n";
103     cout << "Expected a_{1,0} ~ " << 0.5 * sqrt(4.0*PI/3.0) << "\n";
104 
105     return 0;
106 }
107

This example projects a simple function f(θ, φ) = 1 + 0.5 cos θ onto SH up to L=3 using straightforward quadrature. It accounts for the sin θ Jacobian and uses midpoint sampling. The printed coefficients should show a_{0,0} and a_{1,0} close to their analytic values and others near zero, validating orthonormality and implementation.

Time: O(nTheta * nPhi * L^2) for naive quadrature and evaluationSpace: O(L^2) to store coefficients

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	// Utility: constants
5	static const double PI = 3.141592653589793238462643383279502884;
6
7	// Compute associated Legendre P_l^m(x) for l >= 0, \|m\| <= l, with Condon–Shortley phase
8	// Uses stable recurrences. Returns P_l^m(x) for m >= 0. For m < 0, use relation outside.
9	static double associated_legendre(int l, int m, double x) {
10	m = abs(m);
11	if (l < m) return 0.0;
12	// P_m^m(x)
13	double pmm = 1.0;
14	if (m > 0) {
15	double somx2 = sqrt(max(0.0, 1.0 - x*x));
16	double fact = 1.0;
17	for (int i = 1; i <= m; ++i) {
18	pmm = -fact somx2; // includes Condon–Shortley (-1)^m and (2m-1)!! via recurrence
19	fact += 2.0;
20	}
21	}
22	if (l == m) return pmm;
23	// P_{m+1}^m(x) = x (2m+1) P_m^m(x)
24	double pmmp1 = x * (2m + 1) pmm;
25	if (l == m + 1) return pmmp1;
26	// Upward recurrence for l
27	double pll = 0.0;
28	for (int ll = m + 2; ll <= l; ++ll) {
29	pll = ((2ll - 1) x * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
30	pmm = pmmp1;
31	pmmp1 = pll;
32	}
33	return pll;
34	}
35
36	// Normalization N_{lm} = sqrt((2l+1)/(4pi) * (l-m)!/(l+m)!) using log-gamma to avoid overflow
37	static double sh_normalization(int l, int m) {
38	int am = abs(m);
39	double lognum = log(2.0 * l + 1.0) - log(4.0 * PI) + lgamma(l - am + 1.0);
40	double logden = lgamma(l + am + 1.0);
41	return exp(0.5 * (lognum - logden));
42	}
43
44	// Complex spherical harmonic Y_l^m(theta, phi) using physics convention (theta: polar, phi: azimuth)
45	static complex<double> spherical_harmonic(int l, int m, double theta, double phi) {
46	double x = cos(theta);
47	if (m >= 0) {
48	double Plm = associated_legendre(l, m, x);
49	double Nlm = sh_normalization(l, m);
50	complex<double> eimphi = polar(1.0, m * phi); // e^{i m phi}
51	return Nlm * Plm * eimphi;
52	} else {
53	int mp = -m;
54	// Use Y_l^{-m} = (-1)^m * conj(Y_l^{m})
55	complex<double> Ypos = spherical_harmonic(l, mp, theta, phi);
56	double phase = (mp % 2 == 0) ? 1.0 : -1.0; // (-1)^m with m=mp
57	return phase * conj(Ypos);
58	}
59	}
60
61	int main() {
62	// Example: evaluate Y_3^m at theta=1.0 rad, phi=2.0 rad for m=-3..3
63	double theta = 1.0, phi = 2.0;
64	int l = 3;
65	cout.setf(std::ios::fixed); cout << setprecision(8);
66	for (int m = -l; m <= l; ++m) {
67	complex<double> Y = spherical_harmonic(l, m, theta, phi);
68	cout << "Y_" << l << "^" << m << "(theta,phi) = "
69	<< Y.real() << " + i" << Y.imag() << "\n";
70	}
71	// Quick sanity: Y_l^{-m} ≈ (-1)^m * conj(Y_l^{m})
72	complex<double> Ypos = spherical_harmonic(l, 2, theta, phi);
73	complex<double> Yneg = spherical_harmonic(l, -2, theta, phi);
74	complex<double> rhs = ((2 % 2 == 0) ? 1.0 : -1.0) * conj(Ypos);
75	cout << "Check conjugation symmetry (m=2): \|Y_-2 - (-1)^2 conj(Y_2)\| = " << abs(Yneg - rhs) << "\n";
76	return 0;
77	}
78