pyephem 恒星时间给出了意想不到的结果

Question

我是第一次使用ephem，无法理解 oberver.sidereal_time() 的输出

我编写了几个脚本来从小时角度确定太阳时。第一个使用 ephe 计算赤经，并使用 Meeus 的天文算法中的公式来获得格林威治平均恒星时，可以将其转换为经度的局部平均恒星时。

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()

当我插入一些数据时，这给出了我期望的输出：

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

我编写的第二个脚本以相同的方式计算赤经，但使用 ephem 的 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()

这并没有得到我期望的输出。

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

AFAIK，两个脚本之间的唯一区别是第一个基于当地平均恒星时的小时角，而第二个基于当地视恒星时的小时角，这考虑到地球的章动，我认为应该是一个很小的因素。相反，我看到了大约三个小时的差异。谁能给我解释一下这是怎么回事？

Answer 1

当您为 PyEphem 提供一个原始浮点数时，它需要一个角度，然后它相信您首先将角度转换为弧度 - 因为它总是将浮点角度视为弧度，以保持一致。但是在您的第二个脚本中，您将获得以度数表示的经度和纬度，并将它们提供给 PyEphem，就好像它们以弧度为单位一样。如果您添加一个或两个print 语句以查看Observer 的.lon 和.lat 属性是什么样的，您可以看到结果：

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

我认为您想要做的只是将原始经度和纬度字符串提供给 PyEphem，以便它将它们解释为人类可读的度数而不是机器可读的弧度，方法是删除 float() 周围的调用 @987654328 @ 和 argv[4] 在您的第二个脚本中。然后您应该会发现它返回的值更接近您的预期：

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