Geodesic Acceleration 3D function

[cdelv]

Hi, I want to use geodesic acceleration for a 3D function. Actually, it's a special case of a 2D Gaussian:

function twoD_Gaussian(xyz, p)
    x = xyz[:, 1]
    y = xyz[:, 2]
    z = xyz[:, 3]
    amplitude, x0, y0, sigma2 = p
    g = @. amplitude * exp(-0.5*(((x - x0))^2/sigma2 + ((sqrt(y^2 + z^2) - y0))^2/sigma2 ))
    return g[:]
end

However, following the example of the README using Zygot I haven't been able to make it work. This is the best of my attempts

    function Avv!(dir_deriv,p,v)
        for i=1:length(xyz)
            X = xyz[i,:]
            dir_deriv[i] = transpose(v) * Zygote.hessian(z->twoD_Gaussian(X[:,:]',z),p) * v
        end
        return nothing
    end

Can you tell me how to do it for multidimensional functions?

Best regards.