Search code examples
pythonastronomypyephem

How to calculate longitude using PyEphem


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.


Solution

  • 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