docs: def match(lon1, lat1, lon2, lat2, tol=None, nnearest=1):
      """
      Adapted from Eric Tollerud.
      Finds matches in one catalog to another.
   
      Parameters
      lon1 : array-like
          Longitude of the first catalog (degrees)
      lat1 : array-like
          Latitude of the first catalog (shape of array must match `lon1`)
      lon2 : array-like
          Longitude of the second catalog
      lat2 : array-like
          Latitude of the second catalog (shape of array must match `lon2`)
      tol : float or None, optional
          Proximity (degrees) of a match to count as a match.  If None,
          all nearest neighbors for the first catalog will be returned.
      nnearest : int, optional
          The nth neighbor to find.  E.g., 1 for the nearest nearby, 2 for the
          second nearest neighbor, etc.  Particularly useful if you want to get
          the nearest *non-self* neighbor of a catalog.  To do this, use:
          ``spherematch(lon, lat, lon, lat, nnearest=2)``
   
      Returns
      -------
      idx1 : int array
          Indices into the first catalog of the matches. Will never be
          larger than `lon1`/`lat1`.
      idx2 : int array
          Indices into the second catalog of the matches. Will never be
          larger than `lon2`/`lat2`.
      ds : float array
          Distance (in degrees) between the matches
      """
      from scipy.spatial import cKDTree
   
      lon1 = np.asarray(lon1)
      lat1 = np.asarray(lat1)
      lon2 = np.asarray(lon2)
      lat2 = np.asarray(lat2)
   
      if lon1.shape != lat1.shape:
          raise ValueError('lon1 and lat1 do not match!')
      if lon2.shape != lat2.shape:
          raise ValueError('lon2 and lat2 do not match!')
  
      rotator = SphericalRotator(0,0)
  
   
      # This is equivalent, but faster than just doing np.array([x1, y1, z1]).T
      x1, y1, z1 = rotator.cartesian(lon1.ravel(),lat1.ravel())
      coords1 = np.empty((x1.size, 3))
      coords1[:, 0] = x1
      coords1[:, 1] = y1
      coords1[:, 2] = z1
   
      x2, y2, z2 = rotator.cartesian(lon2.ravel(),lat2.ravel())
      coords2 = np.empty((x2.size, 3))
      coords2[:, 0] = x2
      coords2[:, 1] = y2
      coords2[:, 2] = z2
   
      tree = cKDTree(coords2)
      if nnearest == 1:
          idxs2 = tree.query(coords1)[1]
      elif nnearest > 1:
          idxs2 = tree.query(coords1, nnearest)[1][:, -1]
      else:
          raise ValueError('invalid nnearest ' + str(nnearest))
   
      ds = angsep(lon1, lat1, lon2[idxs2], lat2[idxs2])
   
      idxs1 = np.arange(lon1.size)
   
      if tol is not None:
          msk = ds < tol
          idxs1 = idxs1[msk]
          idxs2 = idxs2[msk]
          ds = ds[msk]
   
      return idxs1, idxs2, ds