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]
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.
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
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
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.
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.