pbp_com/apps/almanac/utils/moon_phase.py

""" Function for generating arbitrary moon jds based on a given date.

Copyright (c) 2009, Colin Powell
All rights reserved.

    Basic key
    0 = 'New Moon'
    1 = 'Waxing Crescent Moon'
    2 = 'Quarter Moon'
    3 = 'Waxing Gibbous Moon'
    4 = 'Full Moon'
    5 = 'Waning Gibbous Moon'
    6 = 'Last Quarter Moon'
    7 = 'Waning Crescent'

    Trig key
    0 = New
    4 = Waxing crecent
    8 = Quarter
    15 = Full
    20 = Wanning gibbous
    24 = Last Quarter
    29 = New
"""

import math
from datetime import datetime

def moon_phase(date):
    """ Returns a percentage representing the state of the moon on a given date
    
    """
    date = datetime.strptime(date, "%Y-%m-%d")
    if date.month < 3:
        year=date.year-1
        month=date.month+12
    else:
        year=date.year
        month=date.month
    month+=1
    jd=(365.25*year)+(30.6*month)+(date.day-694039.09)
    jd/=29.53058868
    b=int(jd)
    jd-=b
    #if jd < 4:
    #    b=4-jd
    #else:
    #    b=8-jd
    b=jd*8+0.5
    if b==8:
        b=0
    
    return b

def jdn(date): 
    a = (14 - date.month) // 12
    y = date.year + 4800 - a
    m = date.month + (12 * a) - 3
    p = date.day + (((153 * m) + 2) // 5) + (365 * y)
    q = (y // 4) - (y // 100) + (y // 400) - 32045
    return p + q

def moon_phase_accurate_but_broken(date):
    n=math.floor(12.37 * (date.year-1900 + ((1.0 + date.month-0.5)/12.0)))
    RAD=3.14159265/180.0
    t=n/1236.85
    t2=t*t
    ase=359.2242 + 29.105356*n
    am=306.0253+385.816918*n+0.010730*t2
    xtra=0.75933 + 1.53058868 * n + ((1.178e-4) - (1.55e-7) * t) * t2
    xtra+=(0.1734-3.93e-4*t)*math.sin(RAD*ase)-0.4068*math.sin(RAD*am)
    if xtra > 0.0:
        i=math.floor(xtra)
    else:
        i=math.ceil(xtra - 1.0)
    j1=jdn(date)
    jd=(2415020+28*n)+i
    p = (j1-jd + 30)%30
    if p < 14.5:
        p=14.5-p
    else:
        p=29-p
    return (p/14.5)