I cannot resolve why I am getting a different azimuth/altitude for a comet calculation between PyEphem and Skyfield, yet the right ascension and declination match.
Code example:
#!/usr/bin/env python3
# External verification:
# https://in-the-sky.org/data/object.php?id=0088P
# https://theskylive.com/88p-info
# https://heavens-above.com/comet.aspx?cid=88P
# https://www.calsky.com/cs.cgi/Comets/3?
import datetime, ephem, io, skyfield.api, skyfield.constants, skyfield.data.mpc
now = datetime.datetime.strptime( "2020-07-22", "%Y-%m-%d" )
latitude = -33
longitude = 151
cometName = "88P/Howell"
# https://minorplanetcenter.net/iau/Ephemerides/Comets/Soft03Cmt.txt
cometDataPyEphem = "88P/Howell,e,4.3838,56.6855,235.9159,3.105737,0.1800765,0.56433120,347.8225,07/21.0/2020,2000,g 11.0,6.0"
# https://minorplanetcenter.net/iau/Ephemerides/Comets/Soft00Cmt.txt
cometDataSkyfield = "0088P 2020 09 26.6241 1.353073 0.564331 235.9159 56.6855 4.3838 20200721 11.0 6.0 88P/Howell MPEC 2019-JE2"
print( "PyEphem:", ephem.__version__ )
print( "Skyfield:", skyfield.__version__ )
# PyEphem
observer = ephem.Observer()
observer.date = ephem.Date( now )
observer.lat = latitude
observer.lon = longitude
body = ephem.readdb( cometDataPyEphem )
body.compute( observer )
print( "PyEphem comet", cometName,
"\n\tLat:", observer.lat,
"\n\tLon:", observer.lon,
"\n\tAlt:", body.alt,
"\n\tRA:", body.ra,
"\n\tDec:", body.dec )
# Skyfield
timeScale = skyfield.api.load.timescale( builtin = True )
topos = skyfield.api.Topos( latitude_degrees = latitude, longitude_degrees = longitude )
ephemeris = skyfield.api.load( "de421.bsp" )
with io.BytesIO( cometDataSkyfield.encode() ) as f:
dataframe = skyfield.data.mpc.load_comets_dataframe( f ).set_index( "designation", drop = False )
sun = ephemeris[ "sun" ]
earth = ephemeris[ "earth" ]
body = sun + skyfield.data.mpc.comet_orbit( dataframe.loc[ cometName ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )
t = timeScale.utc( now.year, now.month, now.day, now.hour, now.minute, now.second )
alt, az, bodyDistanceToEarth = ( earth + topos ).at( t ).observe( body ).apparent().altaz()
ra, dec, bodyDistanceToEarth = ( earth + topos ).at( t ).observe( body ).radec()
print( "Skyfield comet", cometName,
"\n\tAz:", az,
"\n\tAlt:", alt,
"\n\tRA:", ra,
"\n\tDec:", dec )
The RA/Dec match between PyEphem and Skfield and also with several websites (see the comment section at the top). Whilst the Skyfield values for Az/Alt are somewhat close to that given by Heavens Above, I cannot determine why the discrepancy between PyEphem and Skyfield for Az/Alt.
EDIT: Amend the script to print out lat/lon and sample (correct) results:
PyEphem: 3.7.6.0
Skyfield: 1.24
PyEphem comet 88P/Howell
Lat: -33:00:00.0
Lon: 151:00:00.0
Az: 109:52:54.6
Alt: -11:19:11.4
RA: 13:27:08.32
Dec: -9:56:54.6
Skyfield comet 88P/Howell
Lat: -33deg 00' 00.0"
Lon: 151deg 00' 00.0"
Az: 109deg 53' 03.4"
Alt: -11deg 19' 23.4"
RA: 13h 26m 04.51s
Dec: -09deg 50' 37.9"
If you try printing out the PyEphem observer location:
print(observer.lat)
print(observer.lon)
you will find that your observer's coordinates are otherwise than you expect:
-1890:45:38.6
8651:39:45.7
That's because PyEphem interprets floating-point angles as radians, which is how your numbers -33 and 151 are being interpreted. To have PyEphem interpret them as degrees, send strings instead:
observer.lat = str(latitude)
observer.lon = str(longitude)
The script should then show very close agreement between the two libraries.
I now consider it a mistake for PyEphem to have tried making it so “convenient” to switch between degrees and radians, but it cannot be fixed without breaking all existing PyEphem scripts. So I wrote Skyfield, which never makes odd decisions based on type, and always tries to label values degrees or radians so that it's clear to the programmer what units are expected.