Spherical Harmonics

For a full description of the code, please see:

Associated Legendre Polynomials and Spherical Harmonics Computation for Chemistry Applications (2014). Taweetham Limpanuparb and Josh Milthorpe. arXiv: 1410.1748 [physics.chem-ph]

Quick start

The normalized associated Legendre polynomials for an angle θ for all l in 0 <= l <= lmax and all m in -l <= m <= l may be generated using the signature computePlm(θ; lmax), eg.

julia> P = computePlmcostheta(pi/2, lmax = 1)
3-element SHArray(::Vector{Float64}, (ML(0:1, 0:1),)):
  0.3989422804014327
  4.231083042742082e-17
 -0.4886025119029199

The polynomials are ordered with m increasing faster than l, and the returned array may be indexed using (l,m) Tuples as

julia> P[(0,0)]
0.3989422804014327

julia> P[(1,1)] == P[3]
true

Spherical harmonics for a colatitude θ and azimuth ϕ may be generated using the signature computeYlm(θ, ϕ; lmax), eg.

julia> Y = computeYlm(pi/3, 0, lmax = 1)
4-element SHArray(::Vector{Complex{Float64}}, (ML(0:1, -1:1),)):
  0.2820947917738782 + 0.0im
  0.2992067103010745 - 0.0im
 0.24430125595146002 + 0.0im
 -0.2992067103010745 - 0.0im

The returned array may be indexed using (l,m) Tuples as well, as

julia> Y[(1,-1)]
0.2992067103010745 - 0.0im

julia> Y[(1,-1)] == Y[2]
true

Special angles SphericalHarmonics.NorthPole() and SphericalHarmonics.SouthPole() may be passed as θ to use efficient algorithms.

Increased precision

Arguments of BigInt and BigFloat types would increase the precision of the result.

julia> Y = computeYlm(big(pi)/2, big(0), lmax = 1) # Arbitrary precision
4-element SHArray(::Vector{Complex{BigFloat}}, (ML(0:1, -1:1),)):
    0.2820947917738781434740397257803862929220253146644994284220428608553212342207478 + 0.0im
    0.3454941494713354792652446460318896831393773703262433134867073548945156550201567 - 0.0im
 2.679783085063171668225419916118067917387251852939708540164955895366691604430101e-78 + 0.0im
   -0.3454941494713354792652446460318896831393773703262433134867073548945156550201567 - 0.0im

Semi-positive harmonics

For real functions it might be sufficient to compute only the functions for m >= 0. These may be computed by passing the flag m_range = SphericalHarmonics.ZeroTo.

julia> computeYlm(pi/3, 0, lmax = 1, m_range = SphericalHarmonics.ZeroTo)
3-element SHArray(::Vector{Complex{Float64}}, (ML(0:1, 0:1),)):
  0.2820947917738782 + 0.0im
 0.24430125595146002 + 0.0im
 -0.2992067103010745 - 0.0im

Real harmonics

It's also possible to compute real spherical harmonics by passing the flag SHType = SphericalHarmonics.RealHarmonics(), eg.

julia> Y = computeYlm(pi/3, pi/3, lmax = 1, SHType = SphericalHarmonics.RealHarmonics())
4-element SHArray(::Vector{Float64}, (ML(0:1, -1:1),)):
  0.2820947917738782
 -0.3664518839271899
  0.24430125595146002
 -0.21157109383040865

These are faster to evaluate and require less memory to store.

See also

FastTransforms.jl: The function FastTransforms.sphevaluate is faster at evaluating real spherical harmonics for a single mode.

julia> @btime FastTransforms.sphevaluate($(big(pi)/3), $(big(pi)/3), 100, 100)
  153.142 μs (1336 allocations: 72.64 KiB)
-3.801739606943941485088961175328961189010502022112528054751517912248264631529766e-07

julia> @btime SphericalHarmonics.sphericalharmonic($(big(pi)/3), $(big(pi)/3), 100, 100, SphericalHarmonics.RealHarmonics())
  165.932 μs (1439 allocations: 78.01 KiB)
-3.801739606943941485088961175328961189010502022112528054751517912248264631529107e-07

This difference might reduce in the future.