The surface normal methods smoothing length and kernel should be able to be set separately

Papers usually recommend using a higher smoothing_length or smoother kernel with surface normal methods.
Currently this can be only done using this quite inefficient approach:

    if smoothing_length != system.smoothing_length ||
       smoothing_kernel !== system.smoothing_kernel
        # TODO: this is really slow but there is no way to easily implement multiple search radia
        search_radius = compact_support(smoothing_kernel, smoothing_length)
        nhs = PointNeighbors.copy_neighborhood_search(nhs, search_radius,
                                                      nparticles(system))
        PointNeighbors.initialize!(nhs, system_coords, neighbor_system_coords)
    end

This was removed from #539 #584 and #666

Start by locating the surface normal methods and reviewing the removed approaches referenced in #539, #584, and #666. Determine how neighbor search radii could support independent smoothing_length and smoothing_kernel settings without the inefficient neighborhood copy; done means those settings work separately with appropriate surface normals.