def get_skydir_distance_mask(src_skydir, skydir, dist, min_dist=None,
                             square=False, coordsys='CEL'):
    """Retrieve sources within a certain angular distance of an
    (ra,dec) coordinate.  This function supports two types of
    geometric selections: circular (square=False) and square
    (square=True).  The circular selection finds all sources with a given
    angular distance of the target position.  The square selection
    finds sources within an ROI-like region of size R x R where R
    = 2 x dist.

    Parameters
    ----------

    src_skydir : `~astropy.coordinates.SkyCoord` 
      Array of sky directions.

    skydir : `~astropy.coordinates.SkyCoord` 
      Sky direction with respect to which the selection will be applied.

    dist : float
      Maximum distance in degrees from the sky coordinate.

    square : bool
      Choose whether to apply a circular or square selection.

    coordsys : str
      Coordinate system to use when applying a selection with square=True.

    """

    if dist is None:
        dist = 180.

    if not square:
        dtheta = src_skydir.separation(skydir).rad
    elif coordsys == 'CEL':
        dtheta = get_linear_dist(skydir,
                                 src_skydir.ra.rad,
                                 src_skydir.dec.rad,
                                 coordsys=coordsys)
    elif coordsys == 'GAL':
        dtheta = get_linear_dist(skydir,
                                 src_skydir.galactic.l.rad,
                                 src_skydir.galactic.b.rad,
                                 coordsys=coordsys)
    else:
        raise Exception('Unrecognized coordinate system: %s' % coordsys)

    msk = (dtheta < np.radians(dist))
    if min_dist is not None:
        msk &= (dtheta > np.radians(min_dist))
    return msk