Repository navigation
Account for distance properly in aticq and atciqz #10291
Description
Both of these routines include a call to erfa.ld to account for light deflection from the Sun. However, in both routines we ignore any distance information and therefore treat objects in front of the Sun as if they were at large distances.
These routines need updating to fix this problem, and as a result the GCRS-ICRS transforms may need re-working, since the first step in this chain is to drop distance information and treat all objects as if they are at large distances, handling parallax separately.
Activity
mhvk commentedon May 4, 2020
Thanks for raising this to keep track! We may have to start to try to keep
UnitSphericala bit more consistently in place, so that we can go for simpler transformations as needed. Though I guess that is really a separate issue - #10294.mkbrewer commentedon May 4, 2020
I'm not sure what you are getting at here. The only input to ld() that contains a distance is em which you have. p, q and e are all unit vectors. If this is about the difference between p and q, then you should be able to handle that. Taking the example of atciqz(), you are passing it the topocentric position (p), but you also have astrom_eb available to calculate the barycentric position (q). This could be done inside atciqz() or you could modify atciqz() to pass both positions to it. The only issues are that q is actually supposed to be the heliocentic position and the position of the Sun is supposed to be light time delayed to the time that the light ray from the source makes its closest approach. These are really very minor issues however. The Sun is quite close to the solar system barycenter and it doesn't move much in the time that the light ray passes the Sun and arrives at the position of the observer. If you want to be rigorous about this, you'd need to also get the barycentric position of the Sun in order to calculate the correct heliocentric position of the source, but that wouldn't be very difficult either.
mkbrewer commentedon May 4, 2020
Oh never mind. I do see your problem. You have dropped the distance before calling atciqz() so you would need to alter it to take both positions. If you are going to do that, you might as well calculate the proper heliocentric position for the source. It might make a difference for Mercury and Venus. It would most certainly make a difference for the Parker Solar Probe.
mhvk commentedon May 4, 2020
Seeing @mkbrewer's comments, it might make sense to try to be as rigorous as we can at first (something I feel we've been attempting generally in astropy) and see how much worse in performance that is than the quicker route (which perhaps we can special-case).
mkbrewer commentedon May 11, 2020
Just copying this result from issue #10292.
Difference in apparent position of the Moon between my program and astropy without correcting light deflection from the Sun:
Difference in apparent position of the Moon after correcting light deflection from the Sun in my program compared to astropy:
StuartLittlefair commentedon May 12, 2020
great detective work @mkbrewer - thanks a lot!
mkbrewer commentedon May 12, 2020
Just letting you know. In ld(), SOFA says:
I suspected that the Sun moves so slowly wrt the SSB that this would make no difference at all. I went ahead and tested it for October 27, 2017 when Jupiter was 1 degree from the Sun using the following code (pardon my C):
Just as I suspected,
there was absolutely no difference to 9 decimal placesfrom just settingtlight_sun = 0.0.Oh, I found a problem with my plot. There was actually a difference of up to 1 in the 9th decimal place (3.6 uas).
mkbrewer commentedon May 13, 2020
Lest you think this is a small effect, here is the error in the standard SOFA Ldsun() for the transit of Venus on June 6th, 2012.
StuartLittlefair commentedon Aug 14, 2020
Finding time to look at this again and want to check that I've made the right plan before I start hacking.
@mkbrewer would you mind commenting on the following?
First the problem:
aticqandatciqzboth callerfa.ldsun(which assumes a distant star) and not the more widely applicableerfa.ld.aticqandatciqzroutines making a fix difficult.My proposed solution:
atciqandatciqzthat callerfa.ldusing distances if present.Does this sound sensible?
@mhvk these routines could be C code and go as additions to liberfa, to be wrapped by pyerfa. Would that violate the "no divergence from SOFA" philosophy? If so, we can add them as Python code to astropy(we had similar routines in the past).
mhvk commentedon Aug 14, 2020
@StuartLittlefair - looking at
atciq, it just calls 6 other erfa routines in series. I think from a speed perspective, there may not be that much difference between doing that in C or in python (at least for large arrays). So, I would definitely start with that. If it does turn out to make a large difference, my tendency would be to just include a bit of C code in astropy and compile it there (and if that isn't easy to do with the current pyerfa, we'll make it easy!). I would like to avoid adding anything to erfa that is not in sofa.mkbrewer commentedon Aug 19, 2020
This is really quite simple to do. In
icrs_to_gcrs()you just need to add a call toThen just subtract this from the light time delayed position of the solar system body to yield the light time delayed position of the solar system body wrt the Sun. Once you've done this, simply add an additional argument toget_body()to get the light time delayed position of the Sun wrt to the observer.atciqz()to provide this position to use as the third argumentqofld()in place ofpthat is used byldsun().mkbrewer commentedon Aug 19, 2020
Oh, you actually just need to call
_get_apparent_body_position()rather thanget_body()for the light time delayed position of the Sun wrt the solar system barycenter.mkbrewer commentedon Aug 19, 2020
Sending it
obsgeolocfrom the solar system body, of course.mkbrewer commentedon Aug 19, 2020
There's a lot of folderol about the light time delay being the time of closest approach to the Sun. I can give you code to do that if you really want, but I've found that it doesn't make a lick of difference as the Sun's motion wrt to the solar system barycenter in the time it takes light to get from the Sun to Earth is entirely negligible.
mkbrewer commentedon Aug 19, 2020
I see that you are using a modified version of
atciqz()fromutilsthat already callsld(), so you are all set. Just add an additional argumentqcofor the position of the source wrt the Sun and use this as the third argument instead ofpco.mkbrewer commentedon Aug 19, 2020
Of course, you could also convert this position to
ra, decand then convert it back to a unit vector withinatciqz()to be compatible with the way you are doing things already. That way, for a star, you can just callatciqz()with bothra, decarguments set to that of the star.mkbrewer commentedon Aug 20, 2020
Oh darn. It's been so long since I looked at this that I have forgotten. The proper light time delay for the position of the Sun wrt the barycenter is that from the target body to the observer, so it's more complicated. The way that I do this in my code is to look up the position of the target body and the Sun simultaneously using the same time delay for both. This means that you'd need to modify
get_body()to return the light time or add this as an attribute to the target body position. Then you'd need to add an optional argument to_get_apparent_body_position()to provide the light time delay to use when looking up the position of the Sun.mkbrewer commentedon Aug 20, 2020
Another approach would to simply do as I do and look up the heliocentric time delayed position of the target body whenever
get_body()is called an then add that unit vector as an attribute to the returned target body position.mkbrewer commentedon Aug 20, 2020
Yet another approach would be to ignore the light delay and assume that the Sun doesn't move significantly wrt the solar system barycenter during the light travel time from the target body. This is, of course, a simplification, but it would make things much easier as you'd only need to call
get_body_barycentric()for the Sun. I haven't tested this, but I can.mkbrewer commentedon Aug 20, 2020
Oh never mind. I'm just being dumb here. I have a bad habit of posting train of thought messages that I should try to avoid. You have the light delayed distance from the source to the observer from which you can derive the light travel time. So all you need to do is call
get_body_barycentric()with the light time delay to get the proper position of the Sun wrt the solar system barycenter. Then subtracting this from the light delayed position of the target body wrt the solar system barycenter gives youq.