Function surface_normal_kernel 

pub fn surface_normal_kernel(
    positions: &[[f64; 3]],
    densities: &[f64],
    masses: &[f64],
    h: f64,
) -> Vec<[f64; 3]>

Compute approximate surface normals using the color function gradient method.

n_i = sum_j (m_j / rho_j) * nabla W(r_ij, h)

The magnitude of n_i is proportional to the interface curvature.