fixed bugs in get_sun - #4952

Very puzzling about the tests passing by hand vs on travis. Might there be some astrometric state that gets set up in running previous tests that actually adds a wrong path? cc @eteq

Forget the tests passing by hand; this was confusion on my part about which version of astropy I had installed. You should be able to see an updated version of the comment above but if not - the test failures are caused by the GCRS-ICRS conversion behaving oddly for coordinates very near the Sun.

Perhaps this is due to asymmetric handling of light deflection?

On this PR in a larger scope: Since the previous routine was clearly wrong, and since answers will therefore change anyway, how about just making get_sun an alias for get_body('sun', ...)? That way at least one possible source of inconsistency is removed. (and the big scary note in solarsystem.rst can be removed...).

Set milestone to 1.2 at first, but then realised this affects 1.0 too. Obviously, for 1.0.10 we cannot use get_body... I do still think, that we should use get_body for 1.2, which may mean two separate PRs... Maybe this one as is first, then another to that lets get_sun use get_body and is for 1.2 only?

@eteq: what do you think?

How about making this PR use the same code as the default implementation in your PR?

This PR could then be back ported as a bug fix.

In your PR you could then make get_sun into an alias for the routines in solar_system.py. The bug in older versions would get fixed but the results wouldn't suddenly change across versions.

Did I get that right?

Yes, makes sense.

@StuartLittlefair - the light deflection should operate exactly the same for the different routes, as it's the same underlying ERFA function. But that doesn't mean it actually is... Or maybe some part of the transform graph is now sending it through some frame that results in it being inverted (which is not the same)? Will try to investigate that ASAP.

@mhvk - I very much disagree that get_sun should be replaced by get_body. Can discuss more about that later.

Wait, I'm confused about the "your implementation". Is that a reference to something in #4920 ?

Looking into it, i think we are just pushing the limits of the atciqz and aticq algorithms in erfa.

Round trip calculations for directions pointing straight at the sun don't work using these algorithms, but atic13 and atci13 do. These routines are 10 times slower though! Maybe I can come up with a quick test to see if the accurate routines need to be used.

@eteq no, that's not it. Using the more accurate routines doesn't seem to help much.

@eteq -> "your implementation" indeed refers to my erfa implementation in #4920, which does the aberration indirectly by taking into account that the light-travel travel time; with this, the position of the sun is still calculated using epv00.

OK - I don't know why I keep going round in circles on this but I think using atic13 etc fixes this problem. Is there any chance that one of you could test this? If so, I can implement a fix tomorrow where we use the accurate transforms close to the Sun

OK. After fairly extensive testing I now understand that this is caused by the light deflection routine in erfa ldsun. For directions very near to the sun this seems to be numerically unstable which means that the algorithm used to make atciqz and aticq reversible does not succeed.

The following snippet shows this:

In [1]: from astropy.coordinates.builtin_frames.utils import get_jd12

In [2]: from astropy.time import Time

In [3]: from astropy import _erfa as erfa

In [4]: import numpy as np

In [5]: pv = np.array([[0,0,0], [0,0,0]])

In [6]: t = Time("2014-06-25T00:00")

In [7]: jd1, jd2 = get_jd12(t, 'tdb')

In [8]: astrom = erfa.apcs13(jd1, jd2, pv)

In [9]: icrs_ra, icrs_dec = 1.6313043965326768, 0.40839069815051204

In [10]: gcrs_ra, gcrs_dec = erfa.atciqz(icrs_ra, icrs_dec, astrom)

In [11]: icrs_ra_return, icrs_dec_return = erfa.aticq(gcrs_ra, gcrs_dec, astrom) 
In [12]: coo = coord.SkyCoord(icrs_ra, icrs_dec, unit='rad')

In [13]: coo_return = coord.SkyCoord(icrs_ra_return, icrs_dec_return, unit='rad')

In [14]: coo.separation(coo_return).to(u.arcsec)
Out[14]: <Angle 62.246915327364164 arcsec>

There is a mechanism in erfa to avoid this issue. ldsun is a wrapper for ld, which takes a "deflection limiter" parameter. This is set in ldsun so that for directions less than 9 arcsecs from the centre of the sun, the deflection is artificially reduced and drops to zero. If we loosen this limit so that the deflection starts to drop for separations smaller than 90 arcseconds, we ensure ICRS<->GCRS round-trips accurate to around 0.1 arcseconds.

In the latest commits I have relaxed this parameter in erfa. I wonder whether this counts as a "bug" that should be reported upstream to SOFA? I imagine it is both machine and implementation specific so I am thinking perhaps not.

Note that I did not re-write get_sun to use exactly the same code as the default in #4920; this would have meant re-writing a version of _get_apparent_body_position which felt inelegant. The code here and the approximate ephemeris in #4920 give the same answer to much better than one arc second anyway.

@StuartLittlefair - great that you found the origin of the problem. Just so I understand: given that we look exactly at the centre of the Sun, light deflection should be zero. Is the problem is that the code currently gets fed a coordinate which is not exactly zero, and therefore calculates a light deflection that is far too large?

Also, more for practical purposes: I don't really like the idea of changing erfa, but if this is just for 1.0.10, maybe it is OK. For 1.2, we can use get_body without changing erfa, correct?

The problem is slightly more subtle than that.

Basically, the routine ldsun calculates the deflection due to the sun alone. As well as the astrometry context it takes as an argument a vector defining the coordinate direction. When the coordinate direction is very close to the sun's centre the natural direction output by ldsun is extremely sensitive to the input vector.

The routine atciq simply applies ldsun to calculate the light deflection. However, aticqz is supposed to be its inverse; i.e it wants to start with the natural direction and undo the deflection to get the coordinate direction. To "undo" the light deflection it starts by using the natural direction as a guess for the coordinate direction and iterates a few times, updating its guess.

However, the output of ldsun varies a lot, even when the changes in the coordinate direction are comparable to the numerical accuracy of my mac. This means that the algorithm above does not converge, and the light deflection is different between atciq and aticqz.

To get round this issue, ldsun is set to "give up" when the coordinate direction comes too close to the sun, and start artificially reducing the light deflection. It shouldn't actually be non-zero unless you're looking at the Sun bang-on, but this a way to ensure numerical stability. In erfa as it stands this is set to 9 arcseconds but this doesn't seem to be enough.

In 1.2 we will still have this issue. Although the coordinates returned by get_body for the sun are accurate, this is not the problem (the coordinates return by get_sun are accurate in this PR too). The issue is that if you converted them to ICRS and back again you would not end up where you started.

Since I share your caution about changing ERFA I've checked that C code using the SOFA library gives the same bug (it does). I've also contacted SOFA with details of this PR to see if they believe it is a bug that needs fixing and if they approve of the change I've made.

@StuartLittlefair - Thanks for investigating this - the meat of this looks good to me.

I am a bit uneasy about changing ldsun for this because it's a substansive ERFA<->SOFA change. But I looked and see nowhere else that it's used. It will be useful to hear from the SOFA board just what they would want to do about it, but it's hard to know how long that will take... Ideally we would instead have a parameter that could control which version is used, but looking at how ld/ldsun/atciq interconnect I think that would be really awkward so this is still probably best.

But be sure to make a PR against https://github.com/liberfa/erfa that replicates this behavior. That's the "official" copy of ERFA, so if you don't include it there it will be removed from Astropy the next time there's an ERFA update to Astropy. And we definitely don't want a "special Astropy" version of ERFA.

Oh, but just to be clear: the fixes that do not include the ERFA modifications are the most important fixes to get in that area clearly bugs (rather than debatable "features", given that the light deflection thing is still right going in one direction). Am I understanding right that the erfa changes are only for that one test, @StuartLittlefair ? If so, it could be we should merge only the first part of this now (and relax the tolerance on the test_gcrs_altaz_bothroutes test) and wait until we hear from SOFA on the ldsun changes? (Or it comes time for a v1.2 release candidate... whichever comes first)

@eteq - that's right. The erfa changes are only necessary for that one test. We could change the object for the test_gcrs_altaz_bothroutes to be something other than the sun and keep the same tolerance.

Also, an alternative to changing erfa itself would be to create our own versions of atciq and aticqz. Since these are just convenience wrappers to call a few erfa subroutines this would not be hard. We could replace the call to ldsun in those routines with an ld call; this is all ldsun does behind the scenes anyway.

@StuartLittlefair - oh, good point re: making our own at... functions. That seems a bit safer as it doesn't require any touching of ERFA. The main downside is that if SOFA comes back with "oh, right, we totally agree with your change", then it was wasted effort to wrap them... But perhaps if they don't get back to us we should do that? (And in that case it still makes sense to split this PR into the erfa-changing and rest part)

Alright, created #4979 for the first part (which should be backported to 1.0.x), and I re-assigned this one to 1.2. Once #4979 is merged, you can rebase on master.

I just merged #4979, so fine to rebase on master now.

Alright, tests are passing (the coverage decrease is in another subpackage so not a real issue). Merging now - thanks @StuartLittlefair!

P.S. re: @mhvk's

OK, yes, should have noticed that. Just goes to show we should probably expose all of erfa!

My tentative plan is to move the python-wrapping stuff (probably as a separate pyerfa or similar) over to the erfa repository, along side the next ERFA release (probably post Astropy v1.2) - then the bundled machinery in astropy can all be tied to a specific ERFA release, too. That will be a good time to make sure everything in ERFA gets python-wrapped.

@eteq - does this need a changelog entry? Or is it affects-dev? (in the latter case it should be given that label).

Good catch, @astrofrog - it does indeed need a changelog entry. Added to master in a48ea34, and 1.2.x in 34eb2fa