Function solve_poisson_dct 

pub(crate) fn solve_poisson_dct(
    f: &[f64],
    nx: usize,
    ny: usize,
    nz: usize,
    vsx: f64,
    vsy: f64,
    vsz: f64,
) -> Vec<f64>

Solve ∇²u = f under a Neumann boundary condition on the array, in place, via DCT-II.

The DCT-II of a length-N signal is the FFT of its even extension restricted to the original samples, so this is the padded periodic solve with the extension never materialised: working memory is the volume itself plus one axis-length buffer, against 8x the volume in complex doubles for the explicit extension.

Eigenvalues of the second difference under this basis are 2(cos(πk/N) − 1)/h². The DC mode is set to zero, as in the periodic solve — the result is defined up to a constant.