Use new ERFA ecliptic transformation machinery?

[eteq]

As liberfa/erfa#34 points out there's a new SOFA release, which means we'll soon have a new ERFA release to include in astropy.

This new release includes a new pair of functions to go from equatorial to ecliptic and back again. Probably we will want to use this as a new implementation for the ecliptic transformations? An open question is whether we want to just use these in place of the current machinery or if there's some subtle difference in implementation. Will have to see how it goes once a new liberfa is out and packaged up.

If it turns out this gives substantially different answers we may want to consider trying to get this in for 1.2. Won't have time to do that by Friday, though... @astrofrog, do you agree this is worth exploring post feature-freeze? (given that it may actually be fixing a "bug" in some sense of the word in the ecliptic transformation)

cc @mhvk

[astrofrog]

This is a similar discussion to the one in #3344 about Galactic coordinates - whether the ERFA implementation should replace or be an additional one compared to what we already have.

I thought that our current Ecliptic frame is still experimental (we have a note about the accuracy not being tested) so I think we can replace the implementation anytime if it gives better results.

Exploring this post-feature freeze sounds good to me.

[astrojuanlu]

I am willing to do this, and check whether the precision improves (ref #6461). Any advice on how to proceed? Something like coordinates-benchmark_in_the_loop?

[astrojuanlu]

As of today, Astropy uses the precession-nutation matrix (including frame bias) of IAU 2006 (erfa.pnm06a) and the mean obliquity of IAU 2006 (erfa.obl06):

https://github.com/astropy/astropy/blob/v2.0.1/astropy/coordinates/builtin_frames/ecliptic_transforms.py#L24-L28

Whereas SOFA/ERFA uses the precession-bias matrix of IAU 2006 (erfa.pmat06) (obliquity is the same):

https://github.com/liberfa/erfa/blob/v1.4.0/src/ecm06.c#L69-L73

(I'm ignoring the long-term model)

Perhaps this is the source of the error?

[astrofrog]

@Juanlu001 - it is possible to run the coordinates-benchmark locally, so if you make local changes to astropy you could re-run the benchmarks and see if the precision improves. If you have any issues running the coordinates-benchmark, just ping me.

[astrojuanlu]

Going back to this question, which I think is somewhat related to #4268 as well.

As a general course of action, should we

1. use ERFA as much as possible in frame transformations,
2. avoid ERFA where possible in frame transformations, or
3. study case by case?

I guess there are reasons for both 1 and 2, and 3 is in a way equivalent of not making any decision.

[mhvk]

My sense is to use erfa wherever possible. It is meant to be (from) "Standards of Astronomy" after all! And anything else we do should at least match erfa - if not, it is a bug either in our code or in theirs, and should be resolved.

[bsipocz]

And I think we can more prominently make it clear that we use SOFA/ERFA, as it regularly comes up (last time at scipy, but I can't recall any more who asked it, they wasn't an astronomer) whether we compare our solutions to theirs and people are surprised that under the hood we actually use it.

[astrojuanlu]

By the way, for what it's worth, to close this issue all it's needed is to replace

astropy/astropy/coordinates/builtin_frames/ecliptic_transforms.py

Lines 34 to 36 in a68bee5

	rbp = erfa.pmat06(jd1, jd2)
	obl = erfa.obl06(jd1, jd2)*u.radian
	return matrix_product(rotation_matrix(obl, 'x'), rbp)

with

    rmat = erfa.ecm06(jd1, jd2)
    return rmat

The other family of ecliptic frames *TrueEcliptic uses a different nutation model and there's no such function in ERFA.

[homeboy445]

I would like to tackle this task if it's still active?

[astrojuanlu]

The task is still active, and my last comment describes what should be done.