Skip to content

Include a realistic treatment of solar light deflection for nearby objects - #10666

Merged
StuartLittlefair merged 5 commits into
astropy:masterfrom
StuartLittlefair:replace_ldsun
Oct 13, 2020
Merged

StuartLittlefair merged 5 commits into
astropy:masterfrom
StuartLittlefair:replace_ldsun

Conversation

@StuartLittlefair

Copy link
Copy Markdown
Contributor

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.atciqz and astropy.coordinates.builtin_frames.utils.aticq there 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.

from time import perf_counter
import numpy as np
from astropy.time import Time
import astropy.units as u
from astropy.coordinates import EarthLocation, CIRS, ICRS, get_body

times = Time('2020-08-21T00:00') + np.linspace(-3*u.hour, 3*u.hour, 400)
location = EarthLocation.of_site('lapalma')
t0 = perf_counter()
moon = get_body('moon', times, location=location)
moon_icrs = moon.transform_to(ICRS)
moon_cirs = moon.transform_to(CIRS(obstime=moon.obstime))
print(f'Transforms took {perf_counter() - t0:.2f} s')

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:

master_comparison
current_comparison

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

@StuartLittlefair
StuartLittlefair requested a review from mhvk August 18, 2020 13:34
@StuartLittlefair StuartLittlefair added this to the v4.0.2 milestone Aug 18, 2020
@pllim

pllim commented Aug 18, 2020 •

Copy link
Copy Markdown
Member

I don't understand the RTD warnings... Perhaps not related? Yeah, not related.

WARNING: the readthedocs_ext.readthedocs extension does not declare if it is safe for parallel reading, assuming it isn't - please ask the extension author to check and make it explicit
WARNING: doing serial read

@mhvk mhvk left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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...

Comment thread astropy/coordinates/builtin_frames/icrs_cirs_transforms.py
Comment thread astropy/coordinates/builtin_frames/icrs_cirs_transforms.py Outdated
Comment thread astropy/coordinates/builtin_frames/utils.py Outdated
Comment thread astropy/coordinates/builtin_frames/utils.py Outdated
Comment thread astropy/coordinates/builtin_frames/utils.py
@StuartLittlefair

Copy link
Copy Markdown
Contributor Author

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!

Jottings-103

@mhvk mhvk left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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!

Comment thread astropy/coordinates/builtin_frames/utils.py Outdated
# 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)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

@mhvk

mhvk commented Aug 19, 2020

Copy link
Copy Markdown
Contributor

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

@mkbrewer mkbrewer Aug 19, 2020 •

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For other than the time delay, that is.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 at t
  • e: vector from Sun to observer at the time of closest approach of the light ray to the Sun
  • q: vector from Sun to source at the time when the source emitted ray that reaches observer at t

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?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

@mkbrewer

Copy link
Copy Markdown
Contributor

I'll take a look at the data that you posted. That was from the the astropy_topo_test.py script that I sent you rather than the Skyfield script that you posted along with it, correct? While the Moon was an excellent target for testing topocentric positions, it's not very good for testing the issue at hand. I'll try to set up some more appropriate tests. I'm thinking the Venus transit and the conjunction of Jupiter around the 27th of December 2019. Clearly, your schedule on this doesn't jibe with mine, but I'll try to get something together soon.

@StuartLittlefair

Copy link
Copy Markdown
Contributor Author

I'll take a look at the data that you posted. That was from the the astropy_topo_test.py script that I sent you rather than the Skyfield script that you posted along with it, correct? While the Moon was an excellent target for testing topocentric positions, it's not very good for testing the issue at hand. I'll try to set up some more appropriate tests. I'm thinking the Venus transit and the conjunction of Jupiter around the 27th of December 2019. Clearly, your schedule on this doesn't jibe with mine, but I'll try to get something together soon.

Thank you - I really appreciate the guidance. The data I posted was from the astropy_topo_test.py script. I have implemented a 2012 Venus transit test in this PR, but if you have a precise and well known time and separation for the closest approach of Venus to the Sun's centre from a reliable source I could improve that test - the best I found was 554.4" from https://eclipse.gsfc.nasa.gov/OH/transit12.html

On schedule - even disregarding the UK/US timeshift I'm about to take 2 weeks holiday so please don't rush.

@mkbrewer

mkbrewer commented Aug 20, 2020 •

Copy link
Copy Markdown
Contributor

OK. If you have the time, it would be helpful for me if you could post the output of astropy_topo_test.py for the time of the Venus transit Jupiter conjunction. I could then compare your result to mine. Regarding this particular PR, after giving it more thought, I see two inaccuracies one inaccuracy. In the case of icrs_to_gcrs() most relevant to solar system bodies, you are deriving q by adding the heliocentric position of the Earth geocenter observer without time delay to the topocentric position of the target body. This contains the non-time delayed position of the Sun wrt to the solar system barycenter and also offsets the target body position by the position of the observer wrt the Earth geocenter. Instead of these rather complicated changes to atciqz(), here's what I would suggest:

newcart = icrs_coo.cartesian - astrom_eb

light_travel_time = newcart.norm() / speed_of_light
emitted_time = gcrs_frame.obstime - light_travel_time
sun = get_body_barycentric('sun', emitted_time)

qcart = icrs_coo.cartesian - sun

srepr = qcart.represent_as(SphericalRepresentation)
q_ra = srepr.lon.to_value(u.radian)
q_dec = srepr.lat.to_value(u.radian)

srepr = newcart.represent_as(SphericalRepresentation)
i_ra = srepr.lon.to_value(u.radian)
i_dec = srepr.lat.to_value(u.radian)

gcrs_ra, gcrs_dec = atciqz(i_ra, i_dec, q_ra, q_dec, astrom)

@mkbrewer

Copy link
Copy Markdown
Contributor

Of course, in the case of a star, the call to atciqz() would simply be:

gcrs_ra, gcrs_dec = atciqz(i_ra, i_dec, i_ra, i_dec, astrom)

from astropy.time import Time
from astropy.utils import iers
from astropy.utils.exceptions import AstropyWarning
from ..representation import CartesianRepresentation

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not needed anymore, I think.

@mkbrewer

Copy link
Copy Markdown
Contributor

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.
AstropyNewPR_MKB_v5
Here is the new data file.
moon_topo_4_6_2017_mkb_ap_v5_new_pr.txt

@mkbrewer

Copy link
Copy Markdown
Contributor

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.

@mkbrewer

Copy link
Copy Markdown
Contributor

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:
AstropyMaster_MKB_v5_SOFA_ld
Since this shows the difference between Astropy master and my program both using qco = pco, that means that the 0.3 mas offset in RA has nothing to do with light deflection. I have no idea what is causing this. The Moon moves at about 0.5 mas / ms mostly in RA, so this would represent about a 600 us offset in time. I am carrying your DUT1 to 1 us and your light time delay algorithm is set to converge to 10 ns, so those can't be it. I don't see this offset when I compare my program's output to that of Skyfield, so it looks like that's the best that we can do there. Since the only remaining issue is that of the time delay for the Sun and the Moon is obviously the wrong thing to look at there, I looked at the Jupiter conjunction on December 27th 2019. Here is a plot of the difference between using the light time from Jupiter to the observer to offset the Sun and no time delay at all as you are doing.
JupiterNoTimeDelay
Clearly, this is negligible, so I'm willing to concede on that point now also. I was, however still curious about this issue of using the light time from Jupiter to the observer to offset the Sun for q and no time delay for e versus using the time of closest approach of the light ray to the Sun for both q and e as seemingly suggested in the documentation for ld() (Note 3). The ES uses the first approach and they make no mention of any error whatsoever in using the light time delay from the target to the observer to offset the Sun for q. They only mention a slight error in using no time delay for e. I found this very curious as the time delay from Jupiter during this conjunction was 51 minutes as opposed to only 8.2 minutes for the time delay from the Sun. I went ahead and tested this by using the time delay from the Sun for both q and e and to my surprise found that there was absolutely no difference between the two approaches. I thought long and hard about this and finally realized that it's just geometry. To illustrate this, I made a drawing:

@mkbrewer

mkbrewer commented Aug 21, 2020 •

Copy link
Copy Markdown
Contributor

Whoops. I accidentally clicked comment instead of attach file. Anyway, here's the drawing:
LightDeflectionGeometry
In this drawing, P is the target body, E is the Earth and D is the deflecting body moving along a path given by the dashed line. At the time that the light ray leaves P, D is at D', so q' is the vector from P to D' using the light time delay from the target. At the time that the light ray passes the deflecting body, the deflecting body is at D so q and e are the vectors at the the time of closest approach. Then when the light ray reaches Earth, the deflecting body is as D_0, so e_0 represents the non-time delayed vector to Earth. Assuming that D is moving at a constant velocity, it turns out that the ratio of q'to e_0 is exactly the same as that of q to e. Therefore the angle between q' and q is the same as that between e and e_0. Thus q' dot e_0 is equal to q dot e. This is clearly not the case for no time delay at all as q_0 dot e_0 is an entirely different angle.

@mkbrewer

Copy link
Copy Markdown
Contributor

btw, should you decide to do this rigorously as suggested in ld(), I figured out a way to do that without an additional ephemeris lookup. All you would need to do is modify prepare_earth_position_vel() to call get_body_barycentric_posvel('sun', time) and return pv for the Sun. Then you could use the velocity of the Sun to backtrack the Sun's position using the algorithm in ldn().

@mkbrewer

mkbrewer commented Aug 21, 2020 •

Copy link
Copy Markdown
Contributor

Well, not exactly, since the algorithm in ldn() won't work for target bodies that are closer to Earth than the Sun, but this does:

    // SOFA wants the light time delay from the Sun to be when the light ray makes its closest approach to the Sun
    // Topocentric direction and distance (unit vector and km) for the source.
    iauPn(source_pos, &source_dist, source_pos);
    // Light time to Source (days)
    tlight = source_dist / (C_LIGHT * DAYSEC);
    // Dot product of Sun's position wrt observer and source direction from observer.
    tlight_sun = iauPdp(sun_pos, source_pos) / (C_LIGHT * DAYSEC);
    // Need to be careful here as we have two times and want to choose the lesser one.
    tlight_sun = dmin(tlight, tlight_sun);
    // If the source elongation from the Sun is greater than 90 degrees, arrival at Earth is the closest approach.
    tlight_sun = dmax(0.0, tlight_sun); 

@mkbrewer

Copy link
Copy Markdown
Contributor

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

@mkbrewer

Copy link
Copy Markdown
Contributor

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 _apparent_position_in_true_coordinates(). Looking at the code for that, I realized that you are currently using an older version of the BPN matrix provided by pn00a() to do the transform to true equator and equinox of date rather than pnm06a() that we both use normally. I figured that might be the problem, so I modified my program to use pn00a() also. Here is the result:
AstropyMaster_MKB_v5_SOFA_ld_pn00a
There is still a small offset with a non-diurnal slope that I still don't understand, but about 90% of the apparent offset was due to the difference in BPNs. You really should update _apparent_position_in_true_coordinates() to use pnm06a() in order to be compatible with the rest of your code. Do I need to open an issue on this? Or can you just take care of this? It'll only take a few minutes to fix.

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.

@mhvk

mhvk commented Aug 24, 2020

Copy link
Copy Markdown
Contributor

@mkbrewer - thanks for the detailed sleuthing! I did a quick check on replacing pn00a with pn06a and there are no test failures, so at least we didn't bake in incorrect answers. Of course, if we make the correction, we should also include some comparisons with known-to-be-correct values (at a level where we would currently fail). Would you have a suggestion for a good test case?

@mkbrewer

Copy link
Copy Markdown
Contributor

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 (erfa.eors()). Adding this to the apparent right ascension returned from _apparent_position_in_true_coordinates() should agree exactly with your CIRS right ascension. Of course, if you use pn06a() rather than pnm06a() as I strongly suggest (see utils.get_cip()), there will be a slight difference due to the different nutation model that is used. One caveat is that you must test geocentric coordinates only as your topocentric CIRS coordinates are currently incorrect (see issue #10356).

@mkbrewer

mkbrewer commented Aug 31, 2020 •

Copy link
Copy Markdown
Contributor

Note that testing this is purely a formality. If you use pnm06a(), I guarantee that it will work.

@mkbrewer

mkbrewer commented Sep 3, 2020

Copy link
Copy Markdown
Contributor

I managed to get rid of the last tiny offset. This was on my side. My program is normally used to process data sampled at 200 Hz, so I only update the nutations once every 60 seconds. For these tests against astropy, that is the sampling interval, but the problem was that I use a 30 second look ahead so that the nutations are calculated at the midpoint of the interval. Eliminating the look ahead got rid of that offset.
AstropyNewPR_MKB_v6

@mhvk

mhvk commented Sep 5, 2020

Copy link
Copy Markdown
Contributor

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?

@mkbrewer

mkbrewer commented Sep 6, 2020 •

Copy link
Copy Markdown
Contributor

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.
JupiterNoTimeDelayWrongBPN
Here is a script that you can run using this PR to check my prediction: astropy_jupiter_topo_test.py.txt
If you run this script and post the output, I can generate a data file for you to use as a comparison. As you can see from the plot, most of the error is due to the incorrect BPN, but there is also a small signal that is due to the lack of time delay for the Sun. If you choose to correct the lack of time delay, then you can either use the algorithm that I posted in the comment that you quote (which I have edited to correct my error that I noted here: #10666 (comment)), or you can use the SOFA method as I described here: #10666 (comment). The back tracked position of the Sun (using tlight_sun as defined here: #10666 (comment), should be applied to the observer's heliocentric position (astrom['eh']) before calculating q as you do in this PR, as the light delay applies to both the source and the observer.

@StuartLittlefair

Copy link
Copy Markdown
Contributor Author

@mkbrewer - many thanks for your detailed investigations whilst I was away. Given that the stated accuracy of the algorithm used for light deflection is 0.5 mas, I'm tempted to not include the time delay calculations in this pull request. I'm open to being talked round though. Comments @mhvk ?

@mkbrewer

Copy link
Copy Markdown
Contributor

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.

@mhvk

mhvk commented Oct 13, 2020

Copy link
Copy Markdown
Contributor

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.

@StuartLittlefair

Copy link
Copy Markdown
Contributor Author

I agree with @mhvk I'm afraid

@eteq

eteq commented Oct 15, 2020

Copy link
Copy Markdown
Member

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.

@eteq

eteq commented Oct 15, 2020

Copy link
Copy Markdown
Member

(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!)

@eteq eteq modified the milestones: v4.0.4, v4.1 Oct 15, 2020
eteq pushed a commit that referenced this pull request Oct 15, 2020
Include a realistic treatment of solar light deflection for nearby objects
@adrn

adrn commented Oct 15, 2020

Copy link
Copy Markdown
Member

I am fairly neutral on backporting to 4.0.x, but if I was forced to pick, I would agree with @mhvk and @StuartLittlefair!

@eteq

eteq commented Oct 16, 2020 •

Copy link
Copy Markdown
Member

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.

@mhvk

mhvk commented Oct 16, 2020

Copy link
Copy Markdown
Contributor

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).

@eteq

eteq commented Oct 16, 2020 •

Copy link
Copy Markdown
Member

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.

@eteq

eteq commented Oct 16, 2020 •

Copy link
Copy Markdown
Member

(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 😉)

@mhvk

mhvk commented Oct 16, 2020

Copy link
Copy Markdown
Contributor

OK, I'll have a go...

@mhvk

mhvk commented Oct 16, 2020

Copy link
Copy Markdown
Contributor

@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?

mhvk added a commit to mhvk/astropy that referenced this pull request Oct 16, 2020
mhvk pushed a commit to mhvk/astropy that referenced this pull request Oct 16, 2020
@mhvk mhvk mentioned this pull request Oct 16, 2020
@mhvk

mhvk commented Oct 16, 2020

Copy link
Copy Markdown
Contributor

@eteq - see #10885 for the backport. I'm happy to do it more bluntly in v4.1.x itself (i.e., rebase and get rid of your earlier merge instead of reverting it), but we may as well check whether the tests pass.

@eteq

eteq commented Oct 19, 2020

Copy link
Copy Markdown
Member

#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 aticq/aticz are used in several places (I think?)

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.

@mhvk

mhvk commented Oct 19, 2020

Copy link
Copy Markdown
Contributor

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...)

@StuartLittlefair

Copy link
Copy Markdown
Contributor Author

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?

@mkbrewer

Copy link
Copy Markdown
Contributor

Again, I would remind you that this is not a small error. Here is the error for the transit of Mercury on November 11, 2019. The maximum error is 6", whereas the angular diameter of Mercury at the time is 10". I should think that people who study these transits would care about an error of this size.
MercuryTransit
Regarding unexpected behavior, it seems to me that the sooner this is distributed to the wider community, the sooner you will hear about something like that.

vla22 pushed a commit to ska-telescope/katpoint that referenced this pull request Nov 3, 2020
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).
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Account for distance properly in aticq and atciqz

8 participants