Observer.radec_of() incorrect?

[Christoph Deil]

Why doesn't separation2 in the following script have a value of 0.1?
Is this a bug or am I misunderstanding what obs.radec_of() is doing?

import ephem

obs = ephem.Observer()
#obs.lon, obs.lat, obs.elevation = 0, 0, 0
#obs.epoch = '2000'
#obs.date = '2001/01/01'

altitude = 0.1
az1, alt1 = 0, 0
az2, alt2 = 0, altitude

ra1, dec1 = obs.radec_of(az1, alt1)
ra2, dec2 = obs.radec_of(az2, alt2)
print 'ra1, dec1 = ', ra1, dec1
print 'ra2, dec2 = ', ra2, dec2

# I would expect separation == elevation, but this is not the case
print 'altitude = ', altitude
print 'separation1 = ', ephem.separation((az1, alt1), (az2, alt2)).real
print 'separation2 = ', ephem.separation((ra1, dec1), (ra2, dec2)).real

"""
Output:
ra1, dec1 = 0:13:46.67 89:22:17.4
ra2, dec2 = 12:13:42.27 84:29:00.3
altitude = 0.1
separation1 = 0.1
separation2 = 0.107251950977"""

[Brandon Rhodes (brandon-rhodes)]

Your code is asking "if one object appears at the horizon and another at 0.1 radians above it (about 5.729578°) then what is the REAL angle between them on a star chart?" And since the atmosphere artificially curves the image of objects upwards near the horizon, two objects must in reality be farther apart than 0.1 radians to APPEAR a mere 0.1 radians apart in altitude in the sky.

To turn off the effects of the atmosphere, set "obs.pressure = 0" and the separation2 should be very close to 0.1 — I am not finding that it is EXACTLY equal, but that might be because of the trig functions not wanting to work exactly backwards — or is there another slight effect besides atmospheric distortion that I am forgetting about? :)

[Christoph Deil (deil-christoph-s)]

Thanks for the explanation!

After setting pressure to zero I am finding a difference (separation1 - separation2) of 23 arcsec for the one example I gave above.
Is the accuracy of the trigonometry calculations really this bad?
Or is there another real effect?

I am trying to use pyephem to test results from another coordinate package, for that I have to be sure the pyephem results are correct approximately at the arcsec level.

Thank you for your help!

[Brandon Rhodes (brandon-rhodes)]

I will look at the code again and see if I can get exactly 0.1 to return from that function. :) In case I want to do my own comparisons: what is the other coordinate package that you are comparing against?

[Christoph Deil (deil-christoph-s)]

The other coordinate package I am writing unit tests for is the HESS ( http://www.mpi-hd.mpg.de/hfm/HESS/ ) data analysis software package, which is not available outside the HESS collaboration.

I have also started to look at pytmp ( http://phn.github.com/pytpm/ ), which might be a good choice to compare pyephem results against.
TPM / PyTPM is very powerful, has unit tests and good documentation, but it is not as easy to use as xephem / pyephem.

Another test for pyephem I have stumbled upon is https://gist.github.com/1672906, where the pyephem result is compared to GPL Horizon.

If you are interested in adding unit tests to pyephem, I could help with that in the next weeks as I debug / add tests to the HESS software.

[Brandon Rhodes (brandon-rhodes)]

PyEphem already has a few tests written, if you will look inside of its "tests" directory — most of the larger tests compare its output against JPL predictions for the same objects.

[Brandon Rhodes (brandon-rhodes)]

By the way, as I work on this issue, I want to double-check our results. I am getting 4.7 arcseconds of difference here in the separation, not the 23 that you got. Could you share your code? I got "0:00:04.7" by adjusting the last few lines to this code:

s1 = ephem.separation((az1, alt1), (az2, alt2)).real
s2 = ephem.separation((ra1, dec1), (ra2, dec2)).real
print ephem.degrees(s2 - s1)

[Brandon Rhodes (brandon-rhodes)]

Hah! No wonder you got a different number of arcseconds: the number of arcseconds changes through the day, even as I sit here in the Common Grounds coffee shop in Bluffton, Ohio, tracking down this issue! So we would have had to run your test program at exactly equivalent times of day to get the same result.

The two effects that cause points 0.1 radian apart in the apparent sky to be a different distance apart in the "real" sky lying behind it are "nutation" and "aberration".

Aberration accounts for the fact that the Earth's velocity through space makes objects to either side of us look like they are "more in front" of us than they really are, just like driving through snow makes the flakes look like they are mostly in front of you. The amount by which aberration affects any given coordinate differs: a point directly ahead or directly behind retains its coordinate exactly; and any point directly abeam is maximally affected. Unless your two points 0.1 radian apart happened to be along a circle of points that all have the same aberration, you will find that the "real" points behind those apparent points of starlight are some slightly different distance apart than 0.1.

Nutation puzzles me more — I am not sure why it should matter. It involves the fact that the earth's axis wobbles slightly all the time, and so converting an alt-az to a "real star atlas" RA and dec need to edit out this wobble since star atlases pretend that our axis precesses in a uniform circle without wobbles. But why should this affect a pair of coordinates whose time (mjd) is exactly the same? I am not sure. At the moment the effect is huge — about 21.1 arcseconds as I type this. I will have to look into this later, but at least I can now name the two effects that you are seeing. If you comment out one or both of these two lines in "circum.c" "obj_fixed()" and recompile, you should be able to drop the difference between your separations nearly to zero:

 nut_eq(mjed, &ra, &dec);
 ab_eq(mjed, lsn, &ra, &dec);

[Christoph Deil (deil-christoph-s)]

Brandon, thanks for all your work identifying the issue!

So the fact that nutation makes the separation between two points in the sky different, depending on whether the separation is computed in (RA, Dec) or (Alt, Az) spherical coordinates is a bug?

[Christoph Deil (deil-christoph-s)]

Hi Brandon,

there is still the unresolved issue about nutation ("I will have to look into this later ..." in comment #7).
Should I file a new ticket on bitbucket or can you reset this ticket's status?

Thanks!
Christoph

[Christoph Deil (deil-christoph-s)]

One more usage question:
I have an Observer (location and time) set up and now want to compute the alt, az for a given GLON, GLAT coordinate.
Observer does have radec_of(alt, az), but not altaz_of(ra, dec).

I tried to create a FixedBody and then call body.compute(observer); print body.alt, body.az on it, but I didn't know how to set the RaDec or Galactic coordinates of the FixedBody.

Computing Alt / Az for a given Ra / Dec or GLON / GLAT seems a very common thing to do with pyephem, but I couldn't find an example in the docu.

[Brandon Rhodes (brandon-rhodes)]

Christoph, a user has just asked a question on Stack Overflow which also involves radec_of(), but the situation is much simpler: I am going to have to figure out why the south celestial pole is not coming out exactly at the south celestial pole! So that should give me a more constrained problem to solve, that then might fix radec_of in a way that will improve the answer you are getting here.

As to your usage question: is there any chance that you could ask on Stack Overflow? That gives me the chance to keep improving and editing my answer — as well as to have other people weigh in with answers of their own in case I get something wrong — and also will help other people to find the answer themselves in the future. Once we have figured out how to do your coordinate transforms, I can also add some instructions to this document, which I think should include information on taking an RA and dec and changing them into altitude and azimuth:

http://rhodesmill.org/pyephem/coordinates.html

Thanks!

[Brandon Rhodes (brandon-rhodes)]

Oh, and for reference, here is the Stack Overflow question that I just mentioned:

http://stackoverflow.com/questions/11161028/pyephem-angle-clarification/

[Christoph Deil (deil-christoph-s)]

I've made a new question here:

http://stackoverflow.com/questions/11169523/how-to-compute-alt-az-for-given-galactic-coordinate-glon-glat-with-pyephem

There are similar questions already for (ra, dec) <-> (alt, az) instead of (glon, glat) <-> (alt, az):

http://stackoverflow.com/questions/8962426/convert-topocentric-coordinates-azimuth-elevation-to-equatorial-coordinates
http://stackoverflow.com/questions/6870743/calculation-star-position-in-the-sky-pyephem

Apart from the fact that pyephem seems to give incorrect results at the moment I think the docs and / or API for coordinate transformations could be improved a bit.
I've made tickets on bitbucket to discuss this because I presume you want to turn off launchpad at some point (there's no link to it any more from the docs)

https://bitbucket.org/brandon/pyephem/issue/7/documentation-improvements-more-docstrings
https://bitbucket.org/brandon/pyephem/issue/8/add-convenience-methods-altaz_of-or

Thanks for looking at this issue again!

[Brandon Rhodes (brandon-rhodes)]

Thank you for this collection of links, and for the issues to help me keep up!