How to calculate longitude using PyEphem

Question

tried to calculate sun lat and long using PyEphem but not matching with ephemeris
SUN: 2011 MAY 04 We 04 14:46:08 13TA12 = 43 degrees approx (As per website www.findyourfate.com)

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

What could be wrong? How to calculate longitude of planet/sun. Helio/Geocentric.

Answer 1

An interesting question, deserving a detailed answer.

First problem. PyEphem takes its dates in the form YYYY/mm/dd, not 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

(This behaviour seems extremely unhelpful. I reported it to Brandon Craig Rhodes as a bug and starting with PyEphem version 3.7.5.1 this behavior has been corrected.)

Second problem. In PyEphem, hlon is normally the heliocentric longitude of a body (the longitude in a sun-centred co-ordinate system). This makes no sense for the sun. So, as a special undocumented exception, when you look up the hlon and hlat of the Sun you get the heliocentric longitude and latitude of the Earth.

(It would be nice if this were documented. I reported this too and PyEphem 3.7.5.1 now carries documentation that follows my recommendation.)

What you want instead, I believe, is the ecliptic longitude of the sun. You can find the ecliptic coordinates of a body using Pyephem's Ecliptic function:

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

However, findyourfate.com reports "13TA12" (that is, 13°12′ of Taurus, corresponding to PyEphem's 43:12:00). What happened to the missing 0:09? I think this comes down to choice of epoch (that is, how much precession to take into account). By default, PyEphem uses the standard astronomical epoch J2000.0. But findyourfate.com seems to be using epoch-of-date:

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

I hope this is all clear: do ask if not.

---

If you want to generate the whole table in Python, you can do it as in the code below. I don't know an easy way to compute the lunar node using PyEphem, so I haven't implemented that bit. (I expect you can do it by iterative search, but I haven't tried.)

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

Answer 2

There is an extra character that has got to go to get this to work. On the next to the last line in the last answer:

Ecliptic(b).long must change to Ecliptic(b).lon

At least in 2.7 running 32 bit on my 64 bit Windows 7 machine.

Also, it would be nice if the table were printed out in degrees for Ecliptic Longitude. :)