Hi,
Is there a way to get the analytical expression of associated legendre polynomials in term of cos(theta) (see here) ?
Thank you for your help!
What do you mean by “get the analytical expression”? What form of answer are you looking for, and for what purpose?
If you just want to write a function that evaluates P_\ell^m(x) for x = \cos \theta, I think the normal method would be to use a three-term recurrence (i.e. a Clenshaw algorithm). Of course, there is also the explicit series formula (which may be less numerically stable?). (If you want this in symbolic form, you could evaluate the recurrence or series symbolically using a CAS.)
If you want to use the P_\ell^{m}(\cos\theta) = (-1)^m (\sin \theta)^m\ \frac{d^m}{d(\cos\theta)^m}\left(P_\ell(\cos\theta)\right) formula to get an expression in terms of both sines and cosines, couldn’t you just plug it into P_\ell(x)=2^{\ell}\sum_{k=0}^{\ell} x^k\binom{\ell}{k}\binom{(\ell+k-1)/2}{\ell} and get
P_\ell^{m}(\cos\theta) = (-1)^m (\sin \theta)^m 2^{\ell}\sum_{k=m}^{\ell} \frac{k!}{(k-m)!}(\cos \theta)^{k-m}\binom{\ell}{k}\binom{(\ell+k-1)/2}{\ell}
(As usual, in a computer implementation of a series like this, you won’t want to explicitly call factorial functions. Instead, analytical cancel numerators and denominators of the coefficients, then write each term as a recurrence from the previous term involving a small number of multiplications/additions. Of course, if you are doing it in a CAS, with arbitrary-precision arithmetic so you aren’t worrying about overflow, and don’t care about efficiency, you can just evaluate the sum naively, with each term evaluated independently straight from the formula. Conversely, for a numerical implementation, I still suspect you’d be better off with a three-term recurrence.)
I want to evaluate the spherical harmonics to compute the following
\sum\limits_{l=1, |m|\leq l}^{lmax} a_{2l, m}Y_{2l,m}
and its derivative.
I tried several very nice packages (SatelliteToolboxLegendre, AssociatedLegendrePolynomials) which are quite fast but they are 3x slower than a boring but tedious
function foo(theta, phi, Y)
st, ct = sincos(theta)
# sin(phi), cos(phi)
sp, cp = sincos(phi)
s2p, c2p = sincos(2f0 * phi)
s3p, c3p = sincos(3f0 * phi)
s4p, c4p = sincos(4f0 * phi)
s5p, c5p = sincos(5f0 * phi)
s6p, c6p = sincos(6f0 * phi)
s7p, c7p = sincos(7f0 * phi)
s8p, c8p = sincos(8f0 * phi)
st2 = st * st
st3 = st2 * st
st4 = st3 * st
st5 = st4 * st
st6 = st5 * st
st7 = st6 * st
st8 = st7 * st
ct2 = ct * ct
ct3 = ct2 * ct
ct4 = ct3 * ct
ct5 = ct4 * ct
ct6 = ct5 * ct
ct7 = ct6 * ct
ct8 = ct7 * ct
stct = st * ct
ylm[1] = 0.28209479177387814f0
ylm[2] = 0.54627421529f0 * st2 * s2p
ylm[3] = -1.0925484305920792f0 * stct * sp
ylm[4] = 0.31539156525252005f0 * (3f0 * ct2 - 1f0)
ylm[5] = -1.0925484305920792f0 * stct * cp
ylm[6] = 0.54627421529f0 * st2 * c2p
....
end
I want to automate this.
This is what I tried otherwise
import SatelliteToolboxLegendre as STL
function _update_sph_cache!(cache::CacheSPH, θ::𝒯, ::Val{lmax}) where {𝒯, lmax}
STL.legendre!( Val(:full), cache.Λ, θ, lmax; ph_term = true)
STL.dlegendre!(Val(:full), cache.dΛ, θ, cache.Λ, lmax; ph_term = true)
end
function get_sph_acc(θ::𝒯, ϕ, V, cache::CacheSPH, _lmax::Val{lmax}) where {𝒯, lmax}
_update_sph_cache!(cache, θ, _lmax)
Λ = cache.Λ
dΛ = cache.dΛ
Yₗₘ, ∂θYₗₘ, ∂ϕYₗₘ = zero(𝒯), zero(𝒯), zero(𝒯)
# cosmϕ = @SVector [cos(m*ϕ) for m=0:8]
# sinmϕ = @SVector [sin(m*ϕ) for m=0:8]
cosmϕ = cache.cosm
sinmϕ = cache.sinm
cosmϕ[1] = one(𝒯)
sinmϕ[1] = zero(𝒯)
sϕ, cϕ = sincos(ϕ)
for m in 1:lmax
sinmϕ[m+1] = sinmϕ[m]*cϕ + cosmϕ[m]*sϕ
cosmϕ[m+1] = cosmϕ[m]*cϕ - sinmϕ[m]*sϕ
end
s2 = sqrt(𝒯(2))
s4 = 1/sqrt(𝒯(4pi))
s8 = 1/sqrt(𝒯(8pi))
n = 1
@inbounds for l = 0:2:lmax, m = -l:l
factor = (m == 0 ? s4 : s2*s8)
p1 = m < 0 ? sinmϕ[-m+1] : cosmϕ[m+1]
p2 = m < 0 ? cosmϕ[-m+1] : sinmϕ[m+1]
am = abs(m)
Vₙ = V[n]
d = dΛ[l+1, am+1]
L = Λ[l+1, am+1]
Yₗₘ += Vₙ * L * factor * p1
∂θYₗₘ += Vₙ * d * factor * p1
∂ϕYₗₘ -= Vₙ * L * factor * p2 * m
n += 1
end
return Yₗₘ, ∂θYₗₘ, ∂ϕYₗₘ
end
I think the typical thing would be to use a three-term recurrence, so that you get each term of the series from a few operations on the preceding ones, similar to how one evaluates Chebyshev polynomials.
For example, the sphericart software, which includes a Julia implementation SpheriCart.jl, described in Bigi et al (2023), uses a few hard-coded spherical harmonics of low orders combined with a modified recurrence formula for higher orders.
Great link!, thank you