FYI: ERFA wrapper factorization

This is an attempt to factorize most of the code on the python side of the ERFA wrapper. This follows what @mdboom has done in #3164. It reduces the memory footprint, but may have a small performance penalty per ERFA call. I have no precise quantification of these two aspects.

PS: I will be out of touch for the next four weeks. Thought I would push this out, in case somebody wants to look into it. This is a very low priority in comparison to all the other ERFA related PRs...
I'm a little confused why `np.broadcast` is not used (both here and in `modeling.utils`). As that is written in `C`, it should nominally be faster, and that seems to be confirmed by a simple test::
```
In [1]: import numpy as np

In [2]: from astropy.modeling.utils import check_broadcast

In [3]: a = np.arange(125).reshape(5,5,5)

In [4]: b = np.arange(5).reshape(1,5,1)

In [5]: np.broadcast(a,b)
Out[5]: <numpy.broadcast at 0x3651dc0>

In [6]: %timeit np.broadcast(a, b).shape
1000000 loops, best of 3: 481 ns per loop

In [7]: %timeit check_broadcast(a.shape, b.shape)
100000 loops, best of 3: 4.14 µs per loop
```