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?

note: [the `solartime()` function is from my answer](http://stackoverflow.com/a/13425515/4279). It is ok to use it but you must provide proper reference (a link to the answer would be enough) as [the license says](http://creativecommons.org/licenses/by-sa/3.0/) — jfs, Nov 20 '12 at 06:52
btw, [the answer](http://stackoverflow.com/a/13425515/4279) explicitly says that you should use radians for the longitude if you use float. — jfs, Nov 20 '12 at 07:02

score 2 · Accepted Answer · answered Nov 20 '12 at 05:27

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

offtop: could you comment on [the answer](http://stackoverflow.com/a/13425515/4279). Is it correct to use `sun.ra` with `observer.sidereal_time()` there, should it be `sun.a_ra` instead? — jfs, Nov 20 '12 at 07:00
J.F., it depends upon whether you are interested in when it *looks* like the sun is exactly overhead, or whether — behind the distortion of the atmosphere — you are interested in when it *really is* overhead. The plain `.ra` should give you its apparent position once all of the distortions that light undergoes are accounted for. The fancier `.a_ra` gives you the position before applying any of the effects that affect real-world observing. — Brandon Rhodes, Nov 20 '12 at 15:59
Thank you. My intent is too emulate a simple sundial i.e., 12:00 should correspond to the shortest (or no) shadow so it seems `.ra` is correct. — jfs, Nov 20 '12 at 17:36

pyephem sidereal time gives unexpected result

1 Answers1