def bilinear_interpolation_weights(lon, lat, nside, order='ring'):
    """
    Get the four neighbours for each (lon, lat) position and the weight
    associated with each one for bilinear interpolation.

    Parameters
    ----------
    lon, lat : :class:`~astropy.units.Quantity`
        The longitude and latitude values as :class:`~astropy.units.Quantity`
        instances with angle units.
    nside : int
        Number of pixels along the side of each of the 12 top-level HEALPix tiles
    order : { 'nested' | 'ring' }
        Order of HEALPix pixels

    Returns
    -------
    indices : `~numpy.ndarray`
        2-D array with shape (4, N) giving the four indices to use for the
        interpolation
    weights : `~numpy.ndarray`
        2-D array with shape (4, N) giving the four weights to use for the
        interpolation
    """

    lon = lon.to_value(u.rad)
    lat = lat.to_value(u.rad)

    _validate_nside(nside)

    nside = np.asarray(nside, dtype=np.intc)

    result = _core.bilinear_interpolation_weights(lon, lat, nside)
    indices = np.stack(result[:4])
    weights = np.stack(result[4:])

    if _validate_order(order) == 'nested':
        indices = ring_to_nested(indices, nside)

    return indices, weights