Units issue with FK4 transforms with and without e terms - #4293

Huh, I thought that we had this right since the coordinates benchmarks were giving good results. This reminds me, we should run those benchmarks again! (and consider automating the process)

I'm surprised - we have accuracy tests against AST here with a tolerance of 10µas - but if the units really are so wrong above, I'm surprised we agree this well! It would be great if someone could dig into this further.

@astrofrog - note that there are two errors, one in that r is not properly normalized, and one that k has degrees rather than radians as a unit. The former will only affect coordinates with an actual distance -- if most tests are with coordinates on the unit sphere, then those will be fine. The unit of k seems more surprising...

@eteq - while you're on coordinates again, could you have a look at this issue? I'm fairly convinced there is a bug -- in particular in the units of the constant of aberration -- but as @astrofrog notes, it is odd that it would pass all the benchmarks and comparisons (unless of course others are making the same mistake...)

Everything @mhvk has said seems right to me... But I confess I don't really understand the E-terms very well - I think @astrofrog or maybe @adrn actually wrote it? It was definitely prior to Quantity, so that's why the units aren't used there. (I know "git blame" says it was me, but that's from various layers of re-organization...)

And re: the degree/radians question: is it possible that it's just that the terms are small enough that even magnified by a factor of 360/2pi, they're close enough?

Looking over old issues/PRs: I checked the suggestion of @eteq that perhaps the changes introduced by treating k as radians while it was meant to be in degrees were too small to be noticed. Unfortunately, this is not true: test_fk4_no_e_fk4.py and test_galactic_fk4.py fail if I convert k to radians (independent of whether I normalise r properly; curiously, test_fk4_no_e_fk5.py works without problem). Could it be that AST does it wrong as well?

See https://github.com/astropy/astropy/compare/master...mhvk:abberation-e-terms?expand=1 for my trials

@astrofrog, do you have any insight on this from your recent asv/coordinate-benchmark fiddlings?

@eteq - not yet, I haven't had much of a chance to dig into this yet.

@eteq, @astrofrog - looking with fresh eyes after a long time at this FK4 e-terms issue, I see that previous results were mostly correct because two errors compensated each other: while the unit of the constant of aberration was inconsistently taken to be in radians, the eccentricity of the earth's orbit was inconsistently taken to be in degrees. Since only the product of the two enters, the final answer was correct. Anyway, I now made the simplest possible correction (still without explicitly using units), and also ensured that the normalisation of distances was done correctly. Please have a look.

p.s. Since the normalisation was a bug that could have affected coordinates with non-unit distances, I'll add a note to the changelog.

@StuartLittlefair - just realised that since you fixed some aberration issues earlier, you might also be a logical one to have a look at this PR.

@mhvk this looks right to me.

As a double check I tested the implementation here against the transformations on pp170 of Explanatory Supplement to the Astronomical Almanac (Seidelmann, 2005) which is an independent test of sorts since his formulae work directly on the ra and dec properties whilst the astropy routine transforms the CartesianRepresentation . The two methods agree so I think this is 👍

I include my test below as a separate comment in case you want to add as a regression test...

Here's my test:

from astropy import coordinates as coord, units as u
from astropy.tests.helper import assert_quantity_allclose as assert_allclose

fk4 = coord.FK4(10.00039633*u.deg, 40.99492542*u.deg) # M31
equinox = fk4.equinox

Dc = -0.065838*u.arcsec
Dd = 0.35299*u.arcsec
# obliquity of Earth's orbit, B1950, radians
o = 0.40920613517217536

fk4noe_dec = fk4.dec - (Dd*np.cos(fk4.ra) - Dc*np.sin(fk4.ra))*np.sin(fk4.dec) -Dc*np.tan(o)*np.cos(fk4.dec)
fk4noe_ra = fk4.ra - (Dc*np.cos(fk4.ra) + Dd*np.sin(fk4.ra))/(np.cos(fk4.dec)) 

fk4noe = fk4.transform_to(coord.FK4NoETerms)
assert_allclose(fk4noe.ra, fk4noe_ra, rtol=5e-7)
assert_allclose(fk4noe.dec, fk4noe_dec)

We should make sure we update the coordinates-benchmark after this.

@StuartLittlefair - thanks! I added the test. @eteq, @astrofrog, or @adrn - could you please check that the normalisation in the distances is now done correctly? The original equation was clearly wrong, but I'm not sure that my "correction" is OK, since it effectively assumes large distance.

I rebased to get rid of the unrelated failures from a new scipy version. I also looked up the explanatory supplement to convince myself that any conversion from FK4 involving e-terms is assumed to be on the unit sphere, which means that the current code should be correct. @eteq, @astrofrog, or @adrn - would you please do a final check?

Could someone change the milestone? 1.0.10 is already released. (Or maybe the friendly astropybot does check it, too?)

(There are also two other unclosed issues in the milestone, ideally those should be moved on as well.)

@bsipocz - thanks for noticing. It is actually a bit odd that 1.0.10 got released without this being noted! Anyway, I moved the entry in the changelog and rebased for good measure.

@eteq, @astrofrog, or @adrn - I think this really should be OK, but could one of you do a final check?

Also: did you mean to update the astropy_helpers commit?

@mhvk Code looks fine to me -- a few minor notes above. Thanks!

@adrn - thanks for the review, and for noticing the "abritrarily set" precision of the test. This gives a bit more confidence in the routine itself too.

Thanks all for the reviews! Will now merge.