Observer.radec_of() applies corrections inconsistently

[keskitalo]

An observer object can translate local azimuth and altitude coordinates into astrometric RA/Dec by calling Observer.radec_of()

If the observer epoch is set to a special EOD (= epoch of the date) value, then radec_of() will skip trying to precess the RA/Dec to a desired epoch:

if (o->now.n_epoch != EOD) ap_as(&o->now, o->now.n_epoch, &ra, &dec);

The trouble is that ap_as() also attempts to apply stellar aberration and solar deflection corrections. Applying these corrections should not depend on the epoch. Additionally, ap_as() corrects for the nutation between present and requested epochs. The corrections are combined in obj_fixed(). This issue is made less significant by the fact pyephem does not appear to make EOD available directly, although using the numeric value from astro.h still works.

Furthermore, in this answer it is suggested that radec_of() should not apply the aberration correction, seemingly contradicting the behavior now in place. Currently, aberration, deflection and nutation are applied whenever the observer epoch is not EOD.

It was suggested in the above answer that Observer might support different flavors of radec_of(), some applying the aberration correction and some not. This way users would have control over which coordinate distortions are corrected for. An obvious specialization would be with and without the "fixed" source corrections (aberration and deflection). Both specializations should include nutation.

[drbitboy]

I thought stellar aberration depends on the solar system barycenter-relative velocity of the observer; why would that not depend on the epoch?

[keskitalo]

Apologies, I was trying to be concise and managed to confuse the subject. I meant that one should not apply the stellar aberration conditionally but rather every time, even if the epoch is set to EOD. I have revised the wording in the issue.

[drbitboy]

I'm the one that needs to apologize because I misread what you wrote, and in retrospect I see it is quite clear; apparently you write English better than I read it!

[brandon-rhodes]

@keskitalo — It has been a few years since I've looked at that code, so my understanding might be rusty, but: the behavior you are reading in the code is the same as that of the underlying "libastro" library's routine radec2ha() which you can find in libastro-3.7.7/misc.c. It reflects a standard practice in astronomy libraries of supporting two different coordinate systems that both use RA and Dec as their coordinates: "astrometric" RA and Dec, which by definition ignore aberration and deflection, and "apparent" RA and Dec, which include it. When the user asks "libastro" for an RA and Dec, its convention — which PyEphem follows — is to assume astrometric RA and Dec when a normal epoch is specified, but to switch to apparent mode if the special non-epoch EOD is provided.

For a modern description of the difference between astrometric and apparent coordinates, try reading the guide to the Naval Observatory's NOVAS library, which includes exactly the same logic as PyEphem in this case: one code path that includes aberration and deflection, and one which does not:

http://aa.usno.navy.mil/software/novas/novas_c/NOVAS_C3.1_Guide.pdf

Page C-18 is where you'll find the definition that "an astrometric place calculation can be used for some differential measurements; it is the same as virtual place except that light bending and aberration (and refraction) are not computed". There's also an extended discussion in PyEphem's manual, here:

http://rhodesmill.org/pyephem/radec

Let me know if these references help.

[keskitalo]

Thank you for the thorough explanation, @brandon-rhodes. The references were indeed helpful.

Applying the logic you outline there seems to be a difference between NOVAS and PyEphem. NOVAS defines the difference of apparent and astrometric systems as refraction, aberration and other effects (C16-C17). Observer.radec_of() will correct for refraction even if you ask for the apparent position by defining epoch=EOD. Of course one can still explicitly disable the refraction by setting the pressure to zero but should that not be automatic when epoch=EOD? Or is refraction considered special?

[brandon-rhodes]

Good question — that's indeed a difference!

[EDIT: as pointed out below, my next paragraph made no sense, because I said "refraction" but meant "deflection"]

The reason is that you can only include ~refraction~ deflection if you know the distance to the object, because ~refraction~ deflection only occurs as light passes a massive body. An object sitting right next to Saturn in the sky won't have its light bent by Saturn if it's on this side of Saturn from us — only if it's way out beyond Saturn so that its light comes past Saturn. And since radec_of() isn't given a distance, it can't safely assume it's being asked about an object far enough away to have experienced ~refraction~ deflection from major Solar System bodies.

[keskitalo]

I believe I have managed to introduce such a flurry of different effects that we are getting confused here. By refraction I mean the way our atmosphere bends light down, towards the horizon. It also seems that the relativistic deflection is only handled for the Sun which makes since. The effect is a thousand times weaker for even Jupiter. NOVAS lumps the effect with aberration and other effects but PyEphem singles it out.

I fully agree that without additional information, PyEphem cannot know if it should apply the aberration and deflection corrections. Having an extra keyword (e.g. "fixed" for objects beyond the solar system) argument could be useful. With this in mind the atmospheric refraction is indeed special because it should apply to all cases equally.

[brandon-rhodes]

Drat — I did indeed read your word “refraction” and started talking about deflection instead!

So, yes, refraction is indeed special because, as you say, it should be applied to all cases equally. Whether you're asking for the astrometric flavor of RA and dec that would match a sky chart, or the apparent flavor that accounts for the aberration, precession, and nutation induced by the real Earth's motion, you still typically want to step "outside the atmosphere" by reversing the effect of refraction, and if you for some reason don't want to, you can turn it off by setting the pressure to zero.

I note that NOVAS doesn't even seem to include a reverse version of its refract() — and also, to keep things straight, never applies refraction to RA and dec coordinates, but only ever to azimuth coordinates.

Neither library seems to include a routine that could reverse the effects of deflection. (I think I'm now using that word correctly? For light bending as is passes a mass?)

So — given all of the above, is there a particular combination of corrections you need, that you can't achieve in the library currently?

[keskitalo]

I agree with everything you said.

I ended up starting this thread when I was looking for a way to convert Az/Alt into astrometric RA/Dec without aberration and deflection corrections. This would be required if the object of concern was an Earth satellite.

Presently I can disable aberration and deflection with epoch=EOD but then I also loose nutation and precession into a specific coordinate system. One way to handle this would be to tell Obs.radec_of() if the coordinates refer to a fixed (distant) object or not. If this was a keyword argument, it could default to "fixed" to be backwards compatible.

If you think adding the keyword makes sense, I can submit a pull request for a review. I think it would require modest modifications to the libastro facilities. Would that be acceptable?

Additionally, the constant EOD could probably be made available through PyEphem so a user could access it without reading the numeric value from astro.h?

[brandon-rhodes]

Returning to this more than a year later, now that I'm catching back up: I would accept a pull request, yes, but in general find PyEphem to be enough of a jumble that I'd be even more excited about how we could accomplish the conversion over in Skyfield, which does a bit of a better job of exposing the different routines to Python.

I'm going to close this for now, but if you do come up with a pull request adding a keyword argument that makes things backwards-compatible, I'll be happy to look into merging it!