# Associated legendre polynomials

**URL:** <https://discourse.julialang.org/t/associated-legendre-polynomials/132892>\
**Category:** General Usage\
**Tags:** question\
**Created:** [October 5, 2025, 5:38pm UTC](https://discourse.julialang.org/t/associated-legendre-polynomials/132892 "2025-10-05T17:38:37Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [October 5, 2025, 5:38pm UTC](https://discourse.julialang.org/t/associated-legendre-polynomials/132892/1 "2025-10-05T17:38:37Z")

</div>

Hi,

Is there a way to get the analytical expression of associated legendre polynomials in term of cos(theta) (see [here](https://en.wikipedia.org/wiki/Associated_Legendre_polynomials#Reparameterization_in_terms_of_angles)) ?

Thank you for your help!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 5, 2025, 11:14pm UTC](https://discourse.julialang.org/t/associated-legendre-polynomials/132892/2 "2025-10-05T23:14:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [October 5, 2025, 11:41pm UTC](https://discourse.julialang.org/t/associated-legendre-polynomials/132892/3 "2025-10-05T23:41:49Z")

</div>

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

```julia-auto
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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [October 5, 2025, 11:46pm UTC](https://discourse.julialang.org/t/associated-legendre-polynomials/132892/4 "2025-10-05T23:46:51Z")

</div>

This is what I tried otherwise

```julia-auto
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

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 5, 2025, 11:51pm UTC](https://discourse.julialang.org/t/associated-legendre-polynomials/132892/5 "2025-10-05T23:51:05Z")

</div>

> [@rveltz](#):
>
> I want to evaluate the spherical harmonics to compute a [_Yₗₘ series_]

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](https://github.com/lab-cosmo/sphericart), which includes a Julia implementation SpheriCart.jl, described in [Bigi et al (2023)](https://arxiv.org/abs/2302.08381), uses a few hard-coded spherical harmonics of low orders combined with a modified recurrence formula for higher orders.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/c/cc175981df0991e17756d91a27235290bc7c6572.png)

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [October 6, 2025, 6:12am UTC](https://discourse.julialang.org/t/associated-legendre-polynomials/132892/6 "2025-10-06T06:12:05Z")

</div>

Great link!, thank you