code: def interpolate_bilinear_lonlat(lon, lat, values, order='ring'):
      
      nside = npix_to_nside(values.shape[0])
      indices, weights = bilinear_interpolation_weights(lon, lat, nside, order=order)
      values = values[indices]
      # At this point values has shape (N, M) where both N and M might be several
      # dimensions, and weights has shape (N,), so we need to transpose in order
      # to benefit from broadcasting, then transpose back so that the dimension
      # with length 4 is at the start again, ready for summing.
      result = (values.T * weights.T).T
      return result.sum(axis=0)
docs: Interpolate values at specific longitudes/latitudes using bilinear interpolation
  
      Parameters
      ----------
      lon, lat : :class:`~astropy.units.Quantity`
          The longitude and latitude values as :class:`~astropy.units.Quantity` instances
          with angle units.
      values : `~numpy.ndarray`
          Array with the values in each HEALPix pixel. The first dimension should
          have length 12 * nside ** 2 (and nside is determined automatically from
          this).
      order : { 'nested' | 'ring' }
          Order of HEALPix pixels
  
      Returns
      -------
      result : float `~numpy.ndarray`
          The interpolated values