Function zernike_basis_d_impl 

pub fn zernike_basis_d_impl<P, E, V, const L: usize, const NORM: u8, const N: usize>(
    x: V,
    y: V,
    out: &mut [V; N],
    ddx: &mut [V; N],
    ddy: &mut [V; N],
)
where
    P: Policy,
    E: FloatElement,
    V: FloatVector<Element = E> + CoreMath,

zernike_basis_impl plus the Cartesian gradient of every mode.

out receives the values exactly as zernike_basis_impl produces them, and ddx/ddy receive $\partial Z_n^m/\partial\{x,y\}$ at the same point.

This is what a Shack-Hartmann reconstruction integrates against: the sensor measures wavefront slopes, so the fit matrix is built from the gradient basis rather than the value basis. Prefer it over seeding a Dual<V, 2> and calling the value form, which carries two derivative components through every operation of the whole ladder, where this shares the Q recurrence between the value and both gradients and differentiates only the two factors that actually depend on the point.

The gradient is finite everywhere including the pupil centre, which is the practical payoff of the Cartesian formulation: the polar $\partial_\theta Z/\rho$ is singular there.