Function wrapped_laplacian_neumann 

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

Wrapped Laplacian under a Neumann boundary: at each array face the missing neighbour is the sample itself, so the wrapped difference across the face is zero. This is exactly the half-sample even extension the DCT-II assumes, so pairing it with solve_poisson_dct reproduces the even-extended periodic solve without building the 2x-per-axis extension.