pyephem sidereal time gives unexpected result

Question

I'm using ephem for the first time, and having trouble understanding the output of oberver.sidereal_time()

I've written a couple scripts to determine solar time from hour angle. The first one uses ephem to compute right ascension and a formula from Meeus' Astronomical Algorithms to get Greenwich mean sidereal time, which can be converted to local mean sidereal with the longitude.

import sys
from datetime import datetime, time, timedelta
import ephem

def hour_angle(dt, longit, latit, elev):
    obs = ephem.Observer()
    obs.date = dt.strftime('%Y/%m/%d %H:%M:%S')
    obs.lon = longit
    obs.lat = latit
    obs.elevation = elev
    sun = ephem.Sun()
    sun.compute(obs)
    # get right ascention
    ra = ephem.degrees(sun.g_ra)

    # get sidereal time at greenwich (AA ch12)
    jd = ephem.julian_date(dt)
    t = (jd - 2451545.0) / 36525
    theta = 280.46061837 + 360.98564736629 * (jd - 2451545) \
            + .000387933 * t**2 - t**3 / 38710000

    # hour angle (AA ch13)
    ha = (theta + longit - ra * 180 / ephem.pi) % 360
    return ha

def main():
    if len(sys.argv) != 6:
        print 'Usage: hour_angle.py [YYYY/MM/DD] [HH:MM:SS] [longitude] [latitude] [elev]'
        sys.exit()
    else:
        dt = datetime.strptime(sys.argv[1] + ' ' + sys.argv[2], '%Y/%m/%d %H:%M:%S')
        longit = float(sys.argv[3])
        latit = float(sys.argv[4])
        elev = float(sys.argv[5])

    # get hour angle
    ha = hour_angle(dt, longit, latit, elev)

    # convert hour angle to timedelta from noon
    days = ha / 360
    if days > 0.5:
        days -= 0.5
    td = timedelta(days=days)

    # make solar time
    solar_time = datetime.combine(dt.date(), time(12)) + td
    print solar_time

if __name__ == '__main__':
    main()

This gives output I would expect when I plug in some data:

> python hour_angle_ephem.py 2012/11/16 20:34:56 -122.2697 37.8044 3.0
2012-11-16 12:40:54.697115

The second script I wrote calculates right ascension the same way, but uses ephem's sidereal_time() to get the local apparent sidereal time.

import sys
from datetime import datetime, time, timedelta
import math
import ephem

def solartime(observer, sun=ephem.Sun()):
    sun.compute(observer)
    # sidereal time == ra (right ascension) is the highest point (noon)
    t = observer.sidereal_time() - sun.ra
    return ephem.hours(t + ephem.hours('12:00')).norm  # .norm for 0..24

def main():
    if len(sys.argv) != 6:
        print 'Usage: hour_angle.py [YYYY/MM/DD] [HH:MM:SS] [longitude] [latitude] [elev]'
        sys.exit()
    else:
        dt = datetime.strptime(sys.argv[1] + ' ' + sys.argv[2], '%Y/%m/%d %H:%M:%S')
        longit = float(sys.argv[3])
        latit = float(sys.argv[4])
        elev = float(sys.argv[5])

    obs = ephem.Observer()
    obs.date = dt.strftime('%Y/%m/%d %H:%M:%S')
    obs.lon = longit
    obs.lat = latit
    obs.elevation = elev
    solar_time = solartime(obs)
    print solar_time

if __name__ == '__main__':
    main()

This does not get me the output I would expect.

python hour_angle_ephem2.py 2012/11/16 20:34:56 -122.2697 37.8044 3.0
9:47:50.83

AFAIK, the only difference between the two scripts is that the first bases hour angle on local mean sidereal time, while the second bases hour angle on local apparent sidereal time, which takes into account the nutation of the earth, which I think should be a very small factor. Instead I'm seeing a difference of about three hours. Can anyone explain to me what is going on?

Answer 1

When you provide PyEphem with a raw floating-point number where it expects an angle, then it trusts that you have converted the angle to radians first — since it always treats floating point angles as radians, to keep things consistent. But in your second script, you are getting longitudes and latitudes that are expressed in degrees and providing them to PyEphem as though they are in radians. You can see the result if you add a print statement or two to see what the .lon and .lat attributes of your Observer look like:

print observer.lon  #--> -7005:32:16.0
print observer.lat  #-->  2166:01:57.2

I think that what you instead want to do is simply provide your raw longitude and latitude strings to PyEphem so that it interprets them as human-readable degrees instead of machine-readable radians, by removing the float() calls around argv[3] and argv[4] in your second script. You should then find that it returns a value closer to what you are expecting:

$ python tmp11.py 2012/11/16 20:34:56 -122.2697 37.8044 3.0
12:40:55.59