def METHOD_NAME(par, create_data3D):
    """WavefieldDecomposition and PressureToVelocity reconstruction of 3d data"""
    p3d, vz3d, p3d_minus, p3d_plus = create_data3D
    # decomposition
    p3d_minus_est, p3d_plus_est = WavefieldDecomposition(
        p3d,
        vz3d,
        par["nt"],
        (par["ny"], par["nx"]),
        par["dt"],
        (par["dy"], par["dx"]),
        rho_sep,
        vel_sep,
        nffts=(nfftk, nfftk, nfftf),
        kind=par["kind"],
        critical=critical * 100,
        ntaper=ntaper,
        dottest=True,
        dtype="complex128",
        **dict(damp=1e-10, iter_lim=10, atol=1e-8, btol=1e-8, show=2)
    )
    assert np.linalg.norm(p3d_minus_est - p3d_minus) / np.linalg.norm(p3d_minus) < 3e-1
    assert np.linalg.norm(p3d_plus_est - p3d_plus) / np.linalg.norm(p3d_plus) < 3e-1
    # reconstruction
    PtoVop = PressureToVelocity(
        par["nt"],
        (par["ny"], par["nx"]),
        par["dt"],
        (par["dy"], par["dx"]),
        rho_sep,
        vel_sep,
        nffts=(nfftk, nfftk, nfftf),
        critical=critical * 100.0,
        ntaper=ntaper,
        topressure=False,
    )
    vz3d_plus_est = (PtoVop * p3d_plus_est.ravel()).reshape(
        par["ny"], par["nx"], par["nt"]
    )
    vz3d_minus_est = (PtoVop * p3d_minus_est.ravel()).reshape(
        par["ny"], par["nx"], par["nt"]
    )
    vz3d_est = np.real(vz3d_plus_est - vz3d_minus_est)
    assert np.linalg.norm(vz3d_est - vz3d) / np.linalg.norm(vz3d) < 3e-1