## SphericalHarmonics.jl

Associated Legendre Polynomials and Spherical Harmonics in Julia
Author jishnub
Popularity
8 Stars
Updated Last
7 Months Ago
Started In
August 2019

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

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.