Fix ecliptic transformations - #6508

Hi there @Juanlu001 👋 - thanks for the pull request! I'm just a friendly 🤖 that checks for issues related to the changelog and making sure that this pull request is milestoned and labeled correctly. This is mainly intended for the maintainers, so if you are not a maintainer you can ignore this, and a maintainer will let you know if any action is required on your part 😃.

Everything looks good from my point of view! 👍

If there are any issues with this message, please report them here.

Thanks! If this change is correct (I'll leave that decision up to @eteq and others), this needs a comment when calling the ERFA function to explain the choice of ERFA function (so people don't change it back), and a changelog entry.

Some of the failures are due to timeouts or jobs not starting, but some others are doctest precision differences with get_body_barycentric('moon', t, ephemeris='de432s'). I am not sure how those are related to this PR though.

Talking about the choice of ERFA function, in the absence of nutation I am not sure if we can keep calling these systems TrueEcliptic. Perhaps MeanEcliptic would be more appropriate @eteq?

I wonder if this means we should actually add frames called MeanEcliptic and leave the current ones as-is, and then use the mean ones for the coordinates benchmark?

I added yet another correction that I have been working on for several weeks now and that was causing us great pain in poliastro.

On the one hand, the vector arithmetic for the computation of the position of the Sun was flawed. I attach a sketch to support the corresponding change.

On the other hand, the equinox attribute of HeliocentricTrueEcliptic was being used both to compute the obliquity and the position of the Sun. I left the former, but for the latter I used a new attribute, obstime, consistent with what it's done with GCRS.

Thanks to these changes, we managed to read data from the JPL Small-Body Database Browser, set it to the corrected HeliocentricTrueEcliptic, convert it to ICRS using Astropy machinery and retrieve the same results that HORIZONS gives. You can imagine how proud we are! You can see our work in progress results here:

http://nbviewer.jupyter.org/github/Juanlu001/poliastro/blob/florence-example/docs/source/examples/Catch%20that%20asteroid!.ipynb

Surprisingly (at least on first thought), the coordinates-benchmarks results did not change. I am not sure about the reasons for this, but perhaps it is only checking the coordinates of very distant objects and traslations of ~800 000 km are not significant (whereas small changes in rotations have dramatic impacts). I would like to ask for guidance regarding this specific issue.

If the devs give their OK from an algorithmic point of view, I will continue adding tests, documentation and a changelog entry.

@Juanlu001 - this looks very promising indeed! I have no time to look at this myself (and am also not a real expert anyway), but also @StuartLittlefair, who wrote quite a bit of these parts of the code.

@Juanlu001 these fixes look great, and look right to me.

My only questions are about frame naming, and the separate equinox and obstime parameters.

With frame naming, ecliptic frames normally use either a mean equinox of a standard date (i.e J2000), the mean equinox of a specific date (i.e what you have here) or the true equinox of a specific date (precession + nutation). We definitely shouldn't call this TrueEcliptic any more! My preference is to add your barycentre fix to the TrueEcliptic frame and add a new MeanEcliptic frame which doesn't have nutation in it.

With respect to using separate equinox and obstime parameters I'm not sure I agree with what you've done here. The equinox of the frame should specify the position and obliquity of the Sun IMO. Given that, I'm not sure of the utility of a seperate obstime parameter. Perhaps @eteq can clarify, since he understands the distinction better than I.

My preference is to add your barycentre fix to the TrueEcliptic frame and add a new MeanEcliptic frame which doesn't have nutation in it.

Cool, will do that!

With respect to using separate equinox and obstime parameters I'm not sure I agree with what you've done here. The equinox of the frame should specify the position and obliquity of the Sun IMO.

I don't have an opinion at all on how these frames should be, so I tried to replicate what NASA and JPL do in their public APIs to suit my needs. In my view (after much struggling), Heliocentric* frames should be equivalent to Geocentric frames in this respect - having an osbtime attribute to tell the position of the central body, be it the Earth or the Sun. On the other hand, I cannot remove the equinox since the value for the obliquity is fixed in J2000 and not in the true date as I copied here:

#6461 (comment)

I don't have a problem in removing obstime and using equinox for finding the position of the Sun, but that will render these systems unusable for me and I will have to redefine them.

@Juanlu001

Sorry for the slow conversation; things are extremely hectic for me at the moment. If I understand correctly then, you wish for the '''MeanEcliptic''' frame to use the same epoch for the obliquity as Horizons, but a user-specified epoch for the central body's position? However, there's no such requirement for the '''TrueEcliptic''' frame?

Would the frames be useable if they had one parameter (e.g '''obstime'''), and the TrueEcliptic frame used this for both the obliquity and the position of the Sun, but the MeanEcliptic frame used a reference epoch (not given as a parameter, but documented) for the position and/or obliquity? If not, perhaps if you can show a use case that would fail it would help me understand why both '''equinox''' and '''obstime''' are needed.

Sorry for the slow conversation; things are extremely hectic for me at the moment.

No need to apologize :)

If I understand correctly then, you wish for the '''MeanEcliptic''' frame to use the same epoch for the obliquity as Horizons, but a user-specified epoch for the central body's position?

Correct. In fact, I don't even need a variable epoch for the obliquity - they just use the value of J2000 always. What I do need though is a way to specify the position of the Sun.

However, there's no such requirement for the '''TrueEcliptic''' frame?

I have no way to validate this. My entire journey has been to replicate what JPL Horizons gives (you can do the experiments by yourself), and the web interface only has one ecliptic frame to select:

Ecliptic and Mean Equinox of Reference Epoch

Reference epoch: J2000.0
XY-plane: plane of the Earth's orbit at the reference epoch
          Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76)
X-axis  : out along ascending node of instantaneous plane of the Earth's
          orbit and the Earth's mean equator at the reference epoch
Z-axis  : perpendicular to the xy-plane in the directional (+ or -) sense
          of Earth's north pole at the reference epoch.

Probably because it's the only modern ecliptic frame defined in SPICE:

https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html#Appendix.%20%60%60Built%20in%27%27%20Inertial%20Reference%20Frames

There were already concerns about validation of ecliptic frames (see original issue #6461 and discussion in #3749). This pull request can be considered another data point.

Would the frames be useable if they had one parameter (e.g '''obstime'''), and the TrueEcliptic frame used this for both the obliquity and the position of the Sun,

I won't use the TrueEcliptic so I'm not very concerned about this :)

but the MeanEcliptic frame used a reference epoch (not given as a parameter, but documented) for the position and/or obliquity?

Yes, definitely. If we have a reference epoch (say, J2000) for the obliquity and an obstime parameter for the position of the Sun, this would work perfectly for me. If you agree I can implement exactly this, and leave TrueEcliptic with just one parameter.

Ok, excellent. I would keep the single parameter's name the same as before to minimise API changes. So I would implement two frames, and have '''equinox''' set the position and obliquity in '''TrueEcliptic''' but only the position in '''MeanEcliptic'''

The name MeanEcliptic seems the right one, but I am less sure about using equinox to then define the position of the Sun, especially as in HCRS the position of the Sun is determined from obstime (as noted by @Juanlu001 above).

But perhaps I should first ensure I really understand what these coordinates are:

• For BarycentricMeanEcliptic, only the equinox matters, as this gives the obliquity and thus the orientation of the coordinates relative to ICRS, defining the coordinate system fully, as there are offsets involved.
• For Helio/GeocentricMeanEcliptic, the obliquity defines the plane, but the offset within that plane depends still on the time via the positions of the Sun and Earth.
• For BarycentricTrueEcliptic, the obliquity now depends on time, so this has variable orientation (but no offset).
• For Helio/GeocentricTrueEcliptic, both orientation and offset now depend on time.

With this list, it seems fairly obvious that equinox has to be the time when the obliquity is taken at, as it defines orientation. However, offsets in other frames like GCRS and HCRS are inferred from obstime, so to me that seems most logical here too, i.e., I would thus advocate the *MeanEcliptic frames to have both equinox and obstime (i.e., they become what @Juanlu001 implemented here).

For the *TrueEcliptic frames, it would seem that one has to insist that obstime and equinox are in fact identical (which I think would need to be implemented).

Separately, looking at the links @Juanlu001 posted, the splice reference frames include:

      17  ECLIPJ2000  Ecliptic coordinates based upon the
                      J2000 frame.
 
                      The value for the obliquity of the
                      ecliptic at J2000 is taken from page 114
                      of [7] equation 3.222-1. This agrees with the
                      expression given in [5].
 
      18  ECLIPB1950  Ecliptic coordinates based upon the B1950
                      frame.
 
                      The value for the obliquity of the ecliptic at
                      B1950 is taken from page 171 of [7].

It seems the frame actually wanted here is the equivalent of ECLIPJ2000, i.e., it has a fixed equinox. I'm not quite sure how ECLIPB1950 would fit in here, as it changes both obliquite and reference frame for that obliquity.

Let me also cc @scottransom, @paulray, and @luojing1211 since they may have some insight on ecliptic reference frames via pulsar timing.

@mhvk, you might also want input from @demorest and David Nice (who may not be on github)

Just to clarify the terminology @mhvk: acoording to the appendix of the SOFA Tools for Earth Attitude cookbook, the terms "true" and "mean" refer to the inclusion of the nutation terms in the rotation matrix. As @eteq explained in his original pull request (see #3749 (comment)), he ended up adding the "true" versions and not the "mean" versions.

Fast forward until two weeks ago, I realized that neither HeliocentricTrueEcliptic nor BarycentricTrueEcliptic were suiting my needs. In the beginning I thought it was a problem of accounting for the nutation (see #4873 (comment)), so I fixed it. The effects of this were:

• The coordinates-benchmark numbers improved by several orders of magnitude. This means that the other tools we are comparing to are not accounting for the nutation either. Notice that, without looking at the code, it's not obvious which precession model they are using (if any) since, at J2000, pre and post IAU 2000 conventions agree (by design).
• These frames can no longer be called *TrueEcliptic since, as per the standard naming, they don't account for nutation anymore.
• My problems were still there. When computing the position of the Florence asteroid and comparing it to the Horizons system, I was having errors "of the same order of magnitude than the norm of the position vector of the Sun with respect to the Solar System Barycenter" (see Coordinates from NEOS are in some ecliptic reference frame poliastro/poliastro#222 (comment)) and the nutation change was of little effect.

So I went on investigating, and found that the position of the Sun was the main driver of the errors - in Astropy it was being computed using the equinox parameter (default to J2000), but apparently in NASA/SPICE/Horizons it was computed at the observation time, independently from the value of the obliquity. I fixed the computation of the position of the Sun, and the effects of this were:

• Absolutely no changes in coordinates-benchmark - not a single decimal digit. I still wonder why, but I guess it has to do with the fact that origin offsets do not affect much angular observations.
• The result of fetch Florence position from NASA in "their" ecliptic frame + assign it to our NewHeliocentricMeanEcliptic + convert to Astropy ICRS was exactly equal to fetch Florence position from NASA in "their" ICRS, so I called it a day.

After all this journey, apart from being precise with the True/Mean naming, I don't have a strong opinion on whether equinox and/or obstime should be used to set the position of the Sun (center of the Heliocentric frames) and the plane of the ecliptic (be it mean of date, mean of epoch, true of date or true of epoch). I feel that using equinox for everything might be a bit confusing, but of course it's a minor issue and I would not complain about that.

Of course, by fixing my needs I don't want to screw up the rest of the folks that use this, so I'm happy they are all brought to the conversation :)

@Juanlu001 - Thanks for the summary! Reading back the comments in #3749, it seems to me that we indeed want the Mean variants and that those should have both equinox and obstime for Geo/Helio. For the True variants, I think it simply means equinox is always equal to obstime.

But really @eteq should pipe in as well..

Separately, shouldn't Heliocentric*Ecliptic connect to HCRS, just like Geocentric*Ecliptic connects to GCRS? (and then any future B1950 ecliptic coordinate system would connect to FK4...).

Separately, shouldn't Heliocentric*Ecliptic connect to HCRS, just like Geocentric*Ecliptic connects to GCRS?

+1!

@mhvk @Juanlu001 I think I'm clearer on the issue here now. I think it's fair to say that the only difference between TrueEcliptic and MeanEcliptic frames is whether one includes nutation or not.

On reflection, I can see no problem with including equinox and obstime parameters for mean frames, but perhaps they need very carefully documenting of what each parameter does, since these terms are different to the "epoch" and "date" terminology usually adopted for ecliptic frames?

It is probably also worth documenting somewhere which equinox/epoch to use to get the same results as Horizons....

I think it's fair to say that the only difference between TrueEcliptic and MeanEcliptic frames is whether one includes nutation or not.

Correct!

On reflection, I can see no problem with including equinox and obstime parameters for mean frames, but perhaps they need very carefully documenting of what each parameter does, since these terms are different to the "epoch" and "date" terminology usually adopted for ecliptic frames?

In fact, "of epoch" equals "of date when date is J2000", at least in the books and papers I have read (perhaps others think otherwise).

Anyway, I agree with documenting this very carefully, it's a confusing topic and some pointers to "ground truth" references (like the SOFA cookbook, also USNO Circular 179) are totally worth it.

It is probably also worth documenting somewhere which equinox/epoch to use to get the same results as Horizons....

With these changes, we can note in the docstring of HeliocentricMeanEcliptic that it is the closest one to what Horizons uses. However, for utter precision we should also remove the precession, since they use a pure rotation around the X axis. Using the SpiceyPy library to directly read the NASA SPICE kernels and compute the rotation matrix:

In [14]: tdb = '2017-09-01 12:05:50 TDB'  # Any date will yield the same rotation matrix!

In [15]: et = spice.str2et(tdb)

In [16]: et
Out[16]: 557539550.0

In [17]: rot = spice.pxform('ECLIPJ2000', 'J2000', et)  # After loading the kernels, etc

In [18]: rot
Out[18]: 
array([[ 1.        ,  0.        ,  0.        ],
       [ 0.        ,  0.91748206, -0.39777716],
       [ 0.        ,  0.39777716,  0.91748206]])

In [19]: np.arccos(rot[1, 1]) * u.rad
Out[19]: <Quantity 0.40909280422232897 rad>

In [20]: _.to(u.arcsec)
Out[20]: <Quantity 84381.448 arcsec>

Which is the value that _erfa.obl80 gives:

In [14]: equinox = Time("J2000", scale="tdb")

In [15]: jd1, jd2 = get_jd12(equinox, 'tdb')

In [16]: obl80 = _erfa.obl80(jd1, jd2) * u.radian

In [17]: obl80.to(u.arcsec)
Out[17]: <Quantity 84381.448 arcsec>

In [18]: obl = _erfa.obl06(jd1, jd2) * u.radian

In [19]: obl.to(u.arcsec)
Out[19]: <Quantity 84381.406 arcsec>

I understand however that this might be too application-specific. If we do not include it in Astropy I will include it in poliastro anyway :)

@Juanlu001 - in principle, it seems to me that if it is used in the SPICE kernels, astropy should support it, or at least we should document how to get it. Just to be sure, in the frame defined for SPICE, which is J2000, both precession and nutation should be irrelevant, correct? Hence, the only term left is the "frame bias" - so is it that term that is not applied in ECLIPJ2000? It might still be good to define corresponding ecliptic frames, though I'm a bit unsure about the names. Maybe *centricJ2000Ecliptic? This then would have no equinox, but still have an obstime.

Just as another reference point, looking at PINT's pulsar_ecliptic.py, I see that their (barycentric) ecliptic frame also implements a transformation to ICRS that worries about obliquity only, i.e., that ignores the frame bias. (Not sure if that is in fact intended; @paulray, @demorest?)

@mhvk -- I haven't had time to really read/understand this discussion fully yet. As far as PINT's implementation, the only intent was to be compatible with how TEMPO/TEMPO2 have traditionally defined ecliptic coordinates. As you say, this is a simple rotation using obliquity only; there are several choices for what obliquity value is used, there may be some useful historical info in $TEMPO/ephem/ecliptic.dat.

Just to be sure, in the frame defined for SPICE, which is J2000, both precession and nutation should be irrelevant, correct?

To be clear: for the transformation between J2000 and ECLIPJ2000 either way, as defined by SPICE, both precession and nutation are irrelevant.

Hence, the only term left is the "frame bias" - so is it that term that is not applied in ECLIPJ2000?

I think that correction was introduced in the new IAU resolutions, and from what I understand they are using the old models.

It might still be good to define corresponding ecliptic frames, though I'm a bit unsure about the names. Maybe *centricJ2000Ecliptic? This then would have no equinox, but still have an obstime.

My proposal:

• BarycentricTrueEcliptic & HeliocentricTrueEcliptic with IAU 2006 precession, nutation and obliquity of equinox (default to J2000 TDB) + obstime parameter for the position of the Sun (Heliocentric only). These frames should not be used in coordinates-benchmark, and it should be noted that there are no sources for validation¹.
• BarycentricMeanEcliptic & HeliocentricMeanEcliptic with IAU 2006 precession and obliquity of equinox (default to J2000 TDB) + obstime parameter for the position of the Sun (Heliocentric only). This is what this PR introduces. These frames should be used in coordinates-benchmark².
• EclipticJ2000 with IAU 1976 precession³ and IAU 1980 obliquity of J2000 (no equinox parameter), centered in the Sun with obstime parameter, explicitly for Horizons compatibility. This frame can be validated copying values from SpiceyPy and hardcoding them in the tests.

¹In principle I would be against including frames that cannot be externally validated, but I understand that maybe others want them there anyway.
²It is possible that the other libraries in coordinates-benchmark use the IAU 1976 precession (see ³) and IAU 1980 obliquity, in which case EclipticJ2000 would perform better. It should be straightforward to check that.
³The IAU 1976 precession matrix at J2000 is the identity, which confirms that SPICE uses the pre IAU 2000 models:

In [13]: equinox = Time("J2000", scale="tdb")

In [14]: jd1, jd2 = get_jd12(equinox, 'tdb')

In [15]: _erfa.pmat76(jd1, jd2)
Out[15]: 
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

Does this make sense?

(For those following this discussion by email -- I made some tweaks to my last comment, please open GH issue)

@Juanlu001 - on the Astropy 3.0 changes break out the various frames, I think you should go full steam ahead. I think I also tilt towards obliquity as a frame attribute if we can make it work (you might find you need to add a custom frame attribute class if it's supposed to accept either a string or a degree value, though...) - it's not the best because I think different obliquity models conceptually map better onto different frames, but this is probably one of those "practicality beats purity" cases, given just what you said about how confusing this all is already.

But I still think having one actual True (i.e. obstime=equinox) equinox is necessary, because it's the one that maps best conceptually to what less-informed people think an equatorial system is. It may be more a nod to pedagogy than to precision, but that's still useful.

Thank you very much @eteq and everybody else that engaged in the discussion! I will be glad to start working on this, hopefully before the end of the year.

A year has passed but I did not forget about this :) I hope I can put something together before Astropy 3.2 feature freeze.

Thank you, @Juanlu001 ! Looking forward to it. Happy New Year!

For the interested readers, here's my followup #8394