Include polar angle motion correction in erfa LST calculation?

[bhazelton]

I am a little concerned that the erfa calculation is not quite right because it doesn't include correction for polar motion in the TIO locator (see https://github.com/astropy/astropy/blob/073d359ebd31c6a3cd10755a9717ac6243da7dd5/astropy/time/core.py#L2329). It seems to be a very small effect (~8.4 microarcseconds for the times I tested with) but I think we should probably include it.

This plot shows the difference between the astropy & erfa calculation without the TIO correction:

And this is with that correction included (basically copied it from astropy, we might want a fix that is more pure erfa). I added this block to the erfa LST calculation:

        # Correct for TIO
        from astropy.coordinates.builtin_frames.utils import get_polar_motion
        from astropy.coordinates.matrix_utilities import rotation_matrix
        sp = erfa.sp00(times.tt.jd1, times.tt.jd2)
        xp, yp = get_polar_motion(times)
        # Form the rotation matrix, CIRS to apparent [HA,Dec].
        r = (
            rotation_matrix(longitude, "z")
            @ rotation_matrix(-yp, "x", unit=units.radian)
            @ rotation_matrix(-xp, "y", unit=units.radian)
            @ rotation_matrix(gast_array + sp, "z", unit=units.radian)
        )
        # Solve for angle.
        angle = np.arctan2(r[..., 0, 1], r[..., 0, 0]) << units.radian

        lst_array = np.mod(angle.value, 2.0 * np.pi)[reverse_inds]

[bhazelton]

@kartographer Interested in your thoughts here.

[bhazelton]

I think we should also reconsider our defaults for the astrometry library.