Add search_around_* functions/methods to match coordinates by radius

This PR adds two more "matching" functions and the related `SkyCoord` convinience methods: `search_around_sky` and `search_around_3d`.  These are different from the `match_coordinates_*` functions in that the latter always find the nearest match of one coordinate in another, while these instead find *all* matches within a certain distance.

Note that in #1846 @taldcroft suggested using `Table.group_by` as a convinient way to deal with the fact that there is no longer a one-to-one correspondence between the input coordinates and the output arrays (as in `match_coordinates_*`).  I tried this, but the problem I ran into was that the 2d and 3d distance outputs need to have units, and to be similar to `match_coordinates_*` they should behave like `Quantity`s.  If I use `Table`, it seems to be impossible to do things like:
```
tab = coo1.search_around_sky(coo2, somedistance)
tab['distance'] > 1*u.kpc
```
because table columns don't yet act like `Quantity`s in comparisons.  It's basically a one-line change to switch the output to a `Table`, though, so I'd say merge this now, and then we can add that in later if #1839 makes the above code work (which hopefully it will!).

If merged, this closes #2876 and also closes #1846.

cc @astrofrog @taldcroft @EmilRex
Are the pairs returned in any specific order?