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...
In the case of modeling I can tell you exactly--  For whatever reason `np.broadcast` only works on *existing* Numpy arrays.  I needed a utility to calculate broadcast shapes without necessarily having created all the arrays involved in the operation yet.  In other words, I needed a utility that could compute the shape of the resulting `np.broadcast` just given a list of array shapes, not the arrays themselves.

I think maybe the same thing was needed here.