code: def subpixel(superpix, nside_superpix, nside_subpix):
      
      if nside_superpix==nside_subpix: return superpix
      vec = hp.pix2vec(nside_superpix, superpix)
      radius = np.degrees(2. * hp.max_pixrad(nside_superpix))
      subpix = query_disc(nside_subpix, vec, radius)
      pix_for_subpix = superpixel(subpix,nside_subpix,nside_superpix)
      # Might be able to speed up array indexing...
      return subpix[pix_for_subpix == superpix]
docs: Return the indices of sub-pixels (resolution nside_subpix) within
      the super-pixel with (resolution nside_superpix).
      
      ADW: It would be better to convert to next and do this explicitly