Repository navigation
Include a realistic treatment of solar light deflection for nearby objects - #10666
Conversation
|
I don't understand the RTD warnings... Perhaps not related? Yeah, not related. |
mhvk
left a comment
There was a problem hiding this comment.
A few initial comments - I very much like the cleanup in the actual transformations, but have some worries about the new implementation in utils (but have not quite wrapped my head around it). Sending it now in case it helps...
|
Thanks for the review @mhvk - it made me realise I was using the wrong type of representation to make this work more easily! It took a while to get my head around this - maybe in reviewing you'll find the attached diagram helpful? It will certainly give you an insight into how bad my spatial visualisation is! |
mhvk
left a comment
There was a problem hiding this comment.
This indeed looks a lot better! I do wonder whether the utilities now should not just return SphericalRepresentation, but think that might also get done as part of #10647, so as you see fit.
One in-line comment about perhaps not using representation/quantity (can it really be me!?).
On the more prosaic side: this does need a changelog entry!
| # Find BCRS direction of Sun to object. | ||
| # astrom['eh'] and astrom['em'] contain Sun to observer unit vector, | ||
| # and distance, respectively. | ||
| eh = astrom['em'] * CartesianRepresentation(astrom['eh'], unit=u.au, xyz_axis=-1, copy=False) |
There was a problem hiding this comment.
I love representations and quantities, but I think here I would not use them, since everything is basically in the right format already. So, just
q = (astrom['em'][..., np.newaxis] * astrom['eh']
+ srepr.distance[..., np.newaxis].to_value(u.au) * before)
sundist = np.linalg.norm(q, axis=-1, keepdims=True)
q /= sundist
q[sundist < 1.e-10] = before
(and, as noted above, maybe get the distance out on top, perhaps even in AU?)
There was a problem hiding this comment.
Thanks for the suggestion - and the solution! It really does take me three goes at least to get shapes right 🥇
I stuck with np.where for the last line, since the direct mask struggled with trailing dimensions.
|
p.s. Bug fix for 4.0.2 does seem appropriate, |
| # and distance, respectively. | ||
| eh = astrom['em'] * CartesianRepresentation(astrom['eh'], unit=u.au, xyz_axis=-1, copy=False) | ||
| # apply parallax to find unit vector from Sun to object | ||
| q = eh + astrometric |
There was a problem hiding this comment.
Sorry to be so late here. I just noticed that you had already made an attempt at this. The problem here is that eh is the non-light delayed position of the Earth wrt the Sun. The Earth moves fairly quickly in its orbit, so it isn't valid to combine the the light delayed position of the target body wrt to the observer and the non-light time delayed position of the Earth wrt to the Sun. For the proper way to do this please see my (late) comments at #10291.
There was a problem hiding this comment.
To do this properly from a topocentric perspective, you'd need the position of the observer wrt to the Sun time delayed by the time it takes the light to get from the target body to the observer. This can be large for a target body such as Jupiter near conjunction with the Sun, so it's much better to do this from a barycentric perspective. That takes the Earth's orbital velocity out of the equation.
There was a problem hiding this comment.
For other than the time delay, that is.
There was a problem hiding this comment.
Thanks for the comments @mkbrewer - I must admit I'm getting quite confused by all the different possible time delays. I'm also confused by what you mean when you say a "topocentric perspective". All the vectors in erfa.ld are barycentric - it's just a question of at what time these positions correspond to.
My understanding is below. Assuming the light ray from the body reaches the observer (Earth) at time t:
p: vector from observer to location of source when it emitted ray that reaches observer atte: vector from Sun to observer at the time of closest approach of the light ray to the Sunq: vector from Sun to source at the time when the source emitted ray that reaches observer att
In my code I'm ignoring the motion of the sun around the Barycentre and using the position at time t, but in the terms I use above, could you point out where I'm going wrong?
There was a problem hiding this comment.
Sorry, my comment concerning the Earth's orbital velocity isn't relevant. Both eh and eb are calculated without any time delay. I hadn't thought about this since May and was not fully back up to speed when I said that. Regarding e, the time of closest approach is mainly for gravitational deflection by planets. Jupiter, for instance, can move 12 arcseconds wrt the solar system barycenter in the time that it takes for light to travel from it to Earth. The Sun moves much more slowly. In fact, it is stated in the Explanatory Supplement that just using eh without any time delay can only cause an error in the deflection of up to 0.1 mas which is less than the 0.5 mas accuracy of the formula itself. My main concern here is the lack of time delay in the position of the Sun used in q. Just scaling up that 0.1 mas from 8.3 minutes to the 4 hours that it takes light to travel from Neptune yields 2.9 mas. What I meant by topocentric perspective is that in computing q you are using two vectors that involve the position of the Earth rather than the solar system barycenter. One of these, eh is geocentric, but the other astrometric (now srepr) can be topocentric which can displace the target body by up to 6378 km. For the Moon which has an orbital velocity of 1 km / s about the Earth, that translates into 106 minutes of time or about 53 arcminutes in position relative to the Earth. I don't know by how much these approximations affect the gravitational deflection, but you can find out yourself by implementing the algorithm that I suggested below. It is the algorithm that both the ES and I use, so comparing the results of that algorithm to the results from this PR should give you a good idea of the error.
There was a problem hiding this comment.
I must apologize once more for my error. I forgot that eh is converted to topocentric in apcs(), so my objection to combining that with eb is withdrawn. The only remaining issue is the proper light time delay. For my thoughts on that, see below.
|
I'll take a look at the data that you posted. That was from the the |
Thank you - I really appreciate the guidance. The data I posted was from the On schedule - even disregarding the UK/US timeshift I'm about to take 2 weeks holiday so please don't rush. |
|
OK. If you have the time, it would be helpful for me if you could post the output of |
|
Of course, in the case of a star, the call to |
| from astropy.time import Time | ||
| from astropy.utils import iers | ||
| from astropy.utils.exceptions import AstropyWarning | ||
| from ..representation import CartesianRepresentation |
|
I'm not sure what file of mine that you compared your output to. I was unable to reproduce your plot. I could come close by using the v4 version of my program though. The main difference between v4 and v5 is that v5 uses TDB whereas v4 uses TT. Using TT was an issue for Brandon Rhodes so I updated to TDB. Anyway here is what I get using v5 after making a slight adjustment to my EOP as your DUT1 seems to have changed. I also took out the polar motion adjustment to LST that I use from USNO Circular 179 as you don't use that and it makes for a slight difference in our LSTs which are not used in any of the calculations anyway. This new PR is a major improvement as there is only a small residual in RA left that may be due to your approximations. The residual is actually within the stated 0.5 mas accuracy of the algorithm. I still think that what you are doing is over complicated and that you should simplify it by using my suggested algorithm which would also be much more clear to the user. |
|
Of course, the real tests would be the Venus transit and the Jupiter conjunction on December 27th 2019. If you can provide me with data for those two, I can get a much better idea of the error introduced by your approximations. |
|
The plot that I posted above looked familiar. In fact it is the same as the one I originally posted at issue #10292. Here it is again with the y axis reversed and the y axis label corrected: |
|
btw, should you decide to do this rigorously as suggested in |
|
Well, not exactly, since the algorithm in |
|
With the errors in GCRS pretty much wrapped up, it's time to move on to CIRS. For that I'll post under issue #10356 |
|
When I set up a version of my program to duplicate the behavior of Astropy in transforming to Alt/Az for issue #10356, I was able to get agreement to 20 uas. This made me think some more about the apparent 0.3 mas offset in topocentric apparent RA. In order to get that, the script that I gave you calls btw, those jumps every 2.5 hours are something on your side. I have no idea what, but nothing in my code has any periodicity like that. |
|
@mkbrewer - thanks for the detailed sleuthing! I did a quick check on replacing |
|
Sorry, for some reason, I didn't receive notification of your comment. Unfortunately, I don't have a suggestion for a third party test of this. JPL Horizons isn't a good test as their apparent right ascension is off by about 53 mas and Skyfield doesn't currently handle polar motion correctly (I raised an issue there concerning this, but it has yet to be addressed). However, you can test for internal consistency by extracting the equation of the origins ( |
|
Note that testing this is purely a formality. If you use |
|
Very nice to see such good agreement. But I must admit that I've been focussing too much on a couple of papers for the last week that I at least have lost track of whether this comparison is against the current PR, or against one that also makes changes like you recommended in #10666 (comment). Or was that recommendation just a suggested simplification? |
|
The comparison is against the data that @StuartLittlefair posted along with this PR. However, while the Moon is an excellent source for checking most things, it's rather useless in the context of this PR. The light travel time from the Moon is so small that it makes no difference whatsoever. A much better test for this PR is the Jupiter conjunction on December 27th 2019. During that conjunction, Jupiter passed within 0.1 degrees of the center of the Sun. Here is my best projection of the error in the position of Jupiter using this PR and the outdated BPN that you currently use to calculate apparent position. |
|
Given the magnitude of the effect, I think it is perfectly reasonable to make that approximation. If you do, you should probably note that in the documentation in order to avoid future confusion, however. |
|
Yes, let's get it in 4.1 for sure. On 4.0, I think mis-interpreting a distance is not really that different from a constant being wrong, leap seconds being misapplied, or some calculation only going wrong in the scalar case. So, still seems a bug to me, and worth back-porting. |
|
I agree with @mhvk I'm afraid |
|
Alright, as an immediate item I'll backport this to 4.1 and make a quick PR noting the change in the 4.1 what's new. |
|
(I'm also going to update the milestone while I'm doing that, but we can change it back to 4.0.x if that's the consensus!) |
Include a realistic treatment of solar light deflection for nearby objects
|
I am fairly neutral on backporting to 4.0.x, but if I was forced to pick, I would agree with @mhvk and @StuartLittlefair! |
|
I'm definitievly outvoted, so I can accept these arguments! That said, after backporting to 4.1.x I'm now also concerned there might be some serious backport issues going even further back to 4.0.x. But I'm happy saying our plan is to do 4.0.x and only back out of that if the backport is too complex to be practical. |
|
If it does appear to be too complicated, do ping - I'm happy to re-apply the commits (the one rebase I did was quite unpleasant, but not difficult to keep the implementations apart). |
|
Thanks, @mhvk - actually I'm now seeing a bunch of failures on v4.1.x that I think may be due to #10647 conflicts with this PR? (see e.g. https://travis-ci.org/github/astropy/astropy/jobs/736086646) - the commit of relevance in that branch is 30d8abb . If you see where the problem lies, a fix would be wonderful! (you can either PR straight into the v4.1.x branch or point me to a corrected backport in your own fork). I was just about to dive in but I'll give a bit of time for you to take a look then. |
|
(note this fixing this breakage and backporting #10881 is now the sole blocker for v4.1... At least until someone finds a new bug I need to backport 😉) |
|
OK, I'll have a go... |
|
@eteq - I have to remove the failed merge in 4.1.x; is it better to push that up and then make a new PR against it, or just to replace it with a single commit that does all the backporting and push directly? |
…e_ldsun" This reverts commit 30d8abb.
|
#10885 was great, thanks @mhvk! That said, it brought my original concern about this in an LTS into sharper focus - https://github.com/astropy/astropy/pull/10885/files is basically the diff of what would go into 4.0.x, and some of the subtle changes (particularly I'm thinking of https://github.com/astropy/astropy/pull/10885/files#diff-cf7d1f9b0f2d0f47c0cfec77e1cffb8f417bbac1d46146c4e0dd78913077a0f1R226) that alter behavior in unexpected ways that might make downstream users' lives painful - e.g., I know I've worked around quirks of some of this code in the past, and I wouldn't want to find my workaround suddenly breaking in an LTS... And it might apply even for people who aren't anywhere close to the sun since So to me that's new information I didn't have any more, and renews my concern. But @mhvk, @StuartLittlefair, and @adrn, if you still think I'm over-worrying feel free to just 👎 my comment here and we can stick with the backport-to-4.0.x plan since #10885 should make that easy for the v4.0.x releaser. |
|
I'm still 👍 on backporting... (and think that just using the second commit from #10666 is easier, as it removes one layer of git cherry-pick just doing something arbitrary...) |
|
I am still relatively ambivalent. This is a minor bug fix that will only affect a few calculations and only at a level that hardly any users care about, so I’m not desperate to backport. In principle I think we should but I guess the issue is what if there is some edge case we’ve missed where this code does something unexpected? |
Astropy 4.1 causes some Solar System objects to shift by 5 mas in the body unit tests. Since the tolerance is 1 mas, this causes the tests to fail. This is most likely due to astropy/astropy#10666, which improves the treatment of solar light deflection for nearby objects. Take the new positions as gospel, since the old behaviour is considered to be a bug. Now that Astropy 4.1 is out we can also remove a workaround for 4.1rc1 (fixed in astropy/astropy#10658).









Description
This pull request addresses an issue where the Solar light deflection was always being applied assuming the source was an infinite distance away, so that the observer -> source vector equalled the Sun -> source vector. In the absolute worst case, this causes errors of up to 3 arc seconds. See #10291 for full details.
Implementation was more awkward than I expected, since we need special case treatments for objects which appear to actually be the Sun itself, and objects at the origin of their frames. I've also special cased the example where there is no distance, to reduce performance impact.
Despite the added complexity in
astropy.coordinates.builtin_frames.utils.atciqzandastropy.coordinates.builtin_frames.utils.aticqthere is no major performance impact. Obtaining the GCRS position of the Moon at 400 times and transforming to ICRS and thence to CIRS took no detectable extra time compared to astropy master - as determined using the test code below.Despite these changes, I have not been able to reproduce without any discrepancies the Moon positions supplied by @mkbrewer in #10291. The discrepancies are however reduced:
I've attached the output of @mkbrewer s test script and the script itself:
moon_topo_astropy_thispr.txt
skyfield_topo_test.py.txt
If @mkbrewer has any input as to the remaining discrepancies I'd be extremely grateful to receive it!
Fixes #10291