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