A program to calculate exact time of the four special moon phases.

moon.py

#!/usr/bin/env python3

import sys
from skyfield.api import load, utc
from datetime import datetime


TS_ZERO = datetime(1970, 1, 1)
NAMES = ['朔', '上弦', '望', '下弦']


time_scale = load.timescale()
planets = load('de421.bsp')


def date2ts(date):
    return int((date - TS_ZERO).total_seconds())


def diff(func):
    def get_value(point):
        return (func(point+1)-func(point-1))/2
    return get_value


def find_root(f, initial, value):
    def g(t):
        y = f(t) - value
        if y > 270:
            return y - 360
        else:
            return y
    d = diff(g)
    x0 = initial
    x1 = x0 - g(x0)/d(x0)
    n = 0
    while abs(x0-x1) > 1:
        x0 = x1
        x1 = x1 - g(x1)/d(x1)
        n = n + 1
        if n > 100000:
            raise Exception(
                'unable to find root when value = %d' % value
            )
    return round(x1)


def delta_l(t):
    T = time_scale.utc(datetime.utcfromtimestamp(t).replace(tzinfo=utc))
    earth = planets['earth']
    sun = planets['sun']
    moon = planets['moon']
    _, l1, __ = earth.at(T).observe(sun).ecliptic_latlon()
    _, l2, __ = earth.at(T).observe(moon).ecliptic_latlon()
    return (l2.degrees - l1.degrees) % 360


def main():
    initial_year = int(sys.argv[1])
    t = date2ts(datetime(initial_year-1, 12, 31))
    begin = datetime(initial_year, 1, 1)
    end = datetime(initial_year+1, 1, 1)
    angle = 0
    while datetime.fromtimestamp(t) < end:
        t = find_root(delta_l, t, angle)
        if datetime.fromtimestamp(t) > begin:
            print(NAMES[angle//90] + '\t' + str(datetime.fromtimestamp(t)) )
        t += 7*86400
        angle = (angle + 90) % 360
        if angle == 0:
            print()
    #print(delta_l(ts.now()))


if __name__ == '__main__':
    main()