Function sh_d_impl 

pub fn sh_d_impl<P, E, V, const L: usize, const N: usize, const CS: bool>(
    x: V,
    y: V,
    z: V,
    out: &mut [V; N],
    ddx: &mut [V; N],
    ddy: &mut [V; N],
    ddz: &mut [V; N],
)
where
    P: Policy,
    E: FloatElement + ShConsts<L, N, CS>,
    V: FloatVector<Element = E> + CoreMath,

sh_impl plus the ambient Cartesian gradient of every harmonic.

out receives the values exactly as sh_impl produces them, and ddx/ddy/ddz receive $\partial Y_{\ell m}/\partial\{x,y,z\}$ of the polynomial form at the given (unit) input. See the module docs for what that means off the sphere and how to project to the tangential gradient.

Two passes over an internal q scratch: the pure z-recurrence first, then one combining sweep that emits values and all three derivatives from tabulated ratios, with no recurrences beyond those sh_impl already runs.