Difficulties with RA/Dec and Alt/Azi conversions with pyEphem

Question

I'm trying to go from alt/azi to RA/Dec for a point on the sky at a fixed location, trying out pyEphem. I've tried a couple of different ways, and I get sort of the right answer, within a couple of degrees, but I'm expecting better, and I can't figure out where the problems lie.

I've been using Canopus as a test case (I'm not after stars specifically, so I can't use the in-built catalogue). So in my case, I know that at

stn = ephem.Observer()
# yalgoo station, wa
stn.long = '116.6806'
stn.lat = '-28.3403'
stn.elevation = 328.0
stn.pressure = 0 # no refraction correction.
stn.epoch = ephem.J2000
stn.date = '2014/12/15 14:32:09' #UTC

Stellarium, checked with other web sites tell me Canopus should be at azi, alt '138:53:5.1', '56:09:52.6' or in equatorial RA 6h 23m 57.09s/ Dec. -52deg 41' 44.6" but trying:

cano = ephem.FixedBody()
cano._ra = '6:23:57.1'
cano._dec = '-52:41:44.0'
cano._epoch = ephem.J2000
cano.compute( stn)
print( cano.az, cano.alt)
>>>(53:22:44.0, 142:08:03.0)

about 3 degrees out. I've also tried the reverse,

ra, dec = stn.radec_of('138:53:5.1', '56:09:52.6')
>>>(6:13:18.52, -49:46:09.5)

where I'm expecting 6:23 not 6:13. Turning on refraction correction makes a small difference, but not enough, and I've always understood aberration and nutation were much smaller effects than this offset as well? As a follow up, I've tried manual calculations, based on 'Practical Astronomy with your calculator'; so for dec:

LAT = math.radians(-28.340335)
LON = math.radians(116.680621667)
ALT = math.radians(56.16461)
AZ = math.radians(138.88475)

sinDEC = (math.sin( LAT)*math.sin( ALT) 
        + math.cos( LAT)*math.cos( ALT)*math.cos( AZ) )
DEC = math.asin( sinDEC)
DEC_deg = math.degrees(DEC)
print( 'dec = ', DEC_deg )
>>>('dec = ', -49.776032754148986)

again, quite different from '56:09:52.6', but reasonably close to pyEphem - so now I'm thoroughly confused! So now I'm suspecting the problem is my understanding, rather than pyEphem - could someone enlighten me about the correct way to go do RADEC/ALTAZI conversions, and why things are not lining up?!

Answer 1

First some notes

1. Atmospheric scattering and relative speed between observer and object

   have the maximal error (near horizon) up to 0.6 degree which is nowhere near your error.

2. how can altitude be over 90 degrees?

   you got swapped data for azimut and altitude

3. I put your observer data into mine program and result was similar to yours

   but I visually search for that star instead of putting the coordinates. Result was also about 3-4 degrees off in RA axis

   RA=6.4h Dec=-52.6deg
   azi=142.4deg alt=53.9deg
   

   mine engine is in C++, using Kepler's equation

Now what can be wrong:

1. mine stellar catalog can be wrongly converted

   rotated wrongly with some margin but strongly doubt that it is 3 degrees. Also perspective transforms can add some error while rendering at 750AU distance from observer. I never tested for Southern sky (not visible from mine place).

2. we are using different Earth reference frame then the data you comparing to

   I found out that some sites like NASA Horizon use different reference frame which does not correspond with mine observations. Look here

   • calculate the time when the sun is X degrees below/above the Horizon

   at the start of the answer is link to 2 sites with different reference frames when you compare the result they are off. The second link is corresponding with mine observations the rest is dealing (included source code) with Kepler's equation based Solar system simulation. The other sublinks are also worth looking into.

3. I could have a bug in mine simulation/data

   I have referenced data to this engine which could partially hide any computation errors from mine observer position so handle all above text with taken that it mind.

4. you could use wrong time/Julian date to stellar time conversions

   if your time is off then the angles will not match...

How to resolve this?

pick up your Telescope, set up equatoreal coordinate system/mount to it and measure Ra/Dec Azi/Alt for known (distant) object in reality and compare with computed positions. Only this way you can decide which value is good or wrong (for reference frame you are using). Do this on star not planet !!! Do this on high altitude angles not near Horizon !!!

How to transform between azimutal and equatoreal coordinates

I compute transform matrix Eath representing earth's coordinate system (upper right) in heliocentric coordinate system as global coordinate system (left) then I compute another matrix NEH representing observer on Earth's surface (North,East,High/Altitude ... lower right).

After this it is just a matter of matrix and vector multiplications and conversion between Cartesian and spherical coordinate systems look here:

• Representing Points on a Circular Radar Math approach

for more insight to azimutal coordinates. if you use just that simple equation like in your example then you do not account for many things... The Earth position is computed by Kepler's equation, rotation is given by daily rotation, nutation and precession included.

I use 64 bit floating point values which can create round errors but not that high ...

I use geometric North Pole as observer reference (this could add some serious error near poles).

The biggest thing that can affect this is the speed of light but that account for near earth 'moving' objects like planets not stars (except Sun) because their computed position is visible after some time ... For example Sun-Earth distance is about 8 light minutes so we see the Sun where it was 8 minutes ago. If the effemerides data is geometrical only (not account for this) then this can lead to high errors if not computed properly.

Newer effemerides models use gravity integration instead of Kepler so their data must be geometrical and the final output is then corrected by the time shift ...