Function angular 

fn angular<const N: usize>(
    l: i32,
    m: i32,
    theta: Jet<N>,
    phi: Jet<N>,
    tangential: bool,
) -> [Jet<N>; 3]

The normalized scalar spherical harmonic Y_lm(theta, phi) and, if tangential and l > 0, the angular functions pi = m Y_lm / sin(theta) and tau = dY_lm / dtheta, both divided by sqrt(l (l + 1)): [Y, pi, tau], or [Y, 0, 0] otherwise.

The recurrence runs over normalized associated Legendre functions: it starts from the normalized diagonal term of degree |m| and steps the degree up to l. Every term stays of order one, which avoids the factorial overflow of the unnormalized functions at high orders. For m = 0 the formula for tau used at other orders would divide by sin(theta), so a second recurrence gives the derivative with respect to cos(theta) instead.