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.

Lian

Just guessing, but could it have something to do with the fact, that you entered '2011-05-04' without specific time information? – Cristinecristiona

it is fine, as long as the variation is marginal. but 43 vs 274, is huge. something must be missing. – Lian

Call me crazy, but I thought the convention was -180 <= longitude_in_degrees <= 180 ... – Underdone

@John Machin - it doesn't have to. Some systems mod 360 and the Western hemisphere is 180 .. 359. – Cristinecristiona

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

Commercialism

Thank you. It is perfect answer for me. It solves my problem. However, is Ecliptic translation required only for Sun? or to all planets (Planets also not matching in my case, so I will try to translate them). Thanks. – Lian

Yes, ecliptic longitude seems to be required for all bodies. See code for printing the tables. Can you explain what you are using this for? – Commercialism

yes, I used ecliptic for all bodies, and results are good. I'm trying to use it for my personal research project in stock market, to see the correlation or reality of Gann methods. – Lian

Excellent explanation Brandon. Btw, in the example above, (Ecliptic(sun).lon) returns 0.7513583544261719. To return 43:02:58.8, i think the "print" statement needs to be used – Ninnetta

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. :)

Sisely