如何使用PyEphem计算经度

Question

我试着用PyEphem来计算太阳的纬度和经度，但结果和天文历书不一致。
太阳：2011年5月4日 星期三 04:14:08 13TA12 = 大约43度（根据网站www.findyourfate.com）

a = Sun()
a.compute('2011-05-04')
>>> a.hlon
274:18:49.1
>>> a.hlat
0:00:00.1

这可能出什么问题了？怎么计算行星或太阳的经度呢？是以太阳为中心的还是以地球为中心的？

Answer 1

这里有一个多余的字符需要去掉，才能让这个代码正常工作。

在最后一个答案的倒数第二行：

Ecliptic(b).long 需要改成 Ecliptic(b).lon

至少在我这台64位的Windows 7机器上运行32位的2.7版本时是这样的。

另外，如果能把表格中的黄道经度显示成度数，那就更好了。 :)

Answer 2

这是个有趣的问题，值得详细解答。

第一个问题是，PyEphem 这个库要求日期的格式是 YYYY/mm/dd，而不是 YYYY-mm-dd。

>>> from ephem import *
>>> Date('2011-05-04')
2010/6/26 00:00:00
>>> Date('2011/05/04')
2011/5/4 00:00:00

（这种行为看起来非常不方便。我已经向 Brandon Craig Rhodes 报告了这个问题， 从 PyEphem 版本 3.7.5.1 开始，这个问题已经得到修正。）

第二个问题是，在 PyEphem 中，hlon 通常指的是一个天体的 日心经度（也就是以太阳为中心的坐标系统中的经度）。对于太阳来说，这个概念是没有意义的。因此，作为一个特别的、没有文档说明的例外，当你查询太阳的 hlon 和 hlat 时，你得到的是 地球 的日心经度和纬度。

（如果这个能在文档中说明就好了。我也向他们报告了这个问题，PyEphem 3.7.5.1 现在已经包含了我建议的文档。）

我相信你真正想要的是太阳的 黄道经度。你可以使用 PyEphem 的 Ecliptic 函数来找到一个天体的黄道坐标：

>>> sun = Sun()
>>> sun.compute('2011/05/04')
>>> Ecliptic(sun).lon
43:02:58.9

不过，findyourfate.com 报告的结果是 "13TA12"（也就是金牛座的 13°12′，对应 PyEphem 的 43:12:00）。那缺失的 0:09 到哪里去了呢？我认为这与选择的历元有关（也就是考虑多少岁差）。默认情况下，PyEphem 使用标准的天文学历元 J2000.0。但 findyourfate.com 似乎使用的是 日期历元：

>>> sun.compute('2011/05/04', '2011/05/04')
>>> Ecliptic(sun).lon
43:12:29.0

希望这些都能让你明白，如果不明白请随时问我。

---

如果你想在 Python 中生成整个表格，可以参考下面的代码。我不知道用 PyEphem 计算月球交点的简单方法，所以我没有实现那部分。（我想你可以通过迭代搜索来做到这一点，但我还没有尝试过。）

from ephem import *
import datetime
import itertools
import math

zodiac = 'AR TA GE CN LE VI LI SC SG CP AQ PI'.split()

def format_zodiacal_longitude(longitude):
    "Format longitude in zodiacal form (like '00AR00') and return as a string."
    l = math.degrees(longitude.norm)
    degrees = int(l % 30)
    sign = zodiac[int(l / 30)]
    minutes = int(round((l % 1) * 60))
    return '{0:02}{1}{2:02}'.format(degrees, sign, minutes)

def format_angle_as_time(a):
    """Format angle as hours:minutes:seconds and return it as a string."""
    a = math.degrees(a) / 15.0
    hours = int(a)
    minutes = int((a % 1) * 60)
    seconds = int(((a * 60) % 1) * 60)
    return '{0:02}:{1:02}:{2:02}'.format(hours, minutes, seconds)

def print_ephemeris_for_date(date, bodies):
    date = Date(date)
    print datetime.datetime(*date.tuple()[:3]).strftime('%A')[:2],
    print '{0:02}'.format(date.tuple()[2]),
    greenwich = Observer()
    greenwich.date = date
    print format_angle_as_time(greenwich.sidereal_time()),
    for b in bodies:
        b.compute(date, date)
        print format_zodiacal_longitude(Ecliptic(b).long),
    print

def print_ephemeris_for_month(year, month, bodies):
    print
    print (datetime.date(year, month, 1).strftime('%B %Y').upper()
           .center(14 + len(bodies) * 7))
    print
    print 'DATE  SID.TIME',
    for b in bodies:
        print '{0: <6}'.format(b.name[:6].upper()),
    print
    for day in itertools.count(1):
        try:
            datetuple = (year, month, day)
            datetime.date(*datetuple)
            print_ephemeris_for_date(datetuple, bodies)
        except ValueError:
            break

def print_ephemeris_for_year(year):
    bodies = [Sun(), Moon(), Mercury(), Venus(), Mars(), Jupiter(), 
              Saturn(), Uranus(), Neptune(), Pluto()]
    for month in xrange(1, 13):
        print_ephemeris_for_month(year, month, bodies)
        print