# vim: fileencoding=cp1251

from __future__ import division
import math, time

# Аргумент -- дата, достаточно близкая к новолунию, в годах
# например, 2008.25 -- для новолуния, ближайшего к 1 апреля

# Результат -- момент новолуния в виде unixtime с примерной
# поправкой на високосные секунды (чтобы после его перевода
# в годы..секунды получалось близко к тому, что надо)

# Для расчёта полнолуния величина k в программе должна быть
# полуцелой (число месяцев с половиной). Величины, отличные
# от целых и полуцелых, подставлять в эти формулы нельзя!

def new_moon(guess):

    # Юлианский день 1 января 1970
    start_of_epoch = 2440587.5

    # Округление и синус в градусах
    round = lambda x: int(x+0.5)
    sin = lambda x: math.sin(math.radians(x))

    # Если эту функцию подставить вместо round, будут полнолуния
    unround = lambda x: int(x)+0.5

    # Примерные месяц и столетие
    k = round((guess-1900) * 12.3685)
    T = k / 1236.85

    # Момент среднего новолуния
    JD = 2415020.75933 + 29.53058868 * k \
                       + 0.0001178 * T**2 \
                       - 0.000000155 * T**3 \
                       + 0.00033 * sin(166.56 + 132.87 * T - 0.009173 * T**2) 

    # Средняя аномалия Солнца
    M = 359.2242 + 29.10535608 * k \
                 - 0.0000333 * T**2 \
                 - 0.00000347 * T**3

    # Средняя аномалия Луны
    M1 = 306.0253 + 385.81691806 * k \
                  + 0.0107306 * T**2 \
                  + 0.00001236 * T**3

    # Аргумент широты Луны
    F = 21.2964 + 390.67050646 * k \
                - 0.0016528 * T**2 \
                - 0.00000239 * T**3

    # Загоняем углы в 0..360
    M, M1, F = [x % 360.0 for x in M, M1, F]

    # Поправки
    adjusments = [
            +(0.1734 - 0.000393*T) * sin(M),
            + 0.0021 * sin(2*M),
            - 0.4068 * sin(M1),
            + 0.0161 * sin(2*M1),
            - 0.0004 * sin(3*M1),
            + 0.0104 * sin(2*F),
            - 0.0051 * sin(M+M1),
            - 0.0074 * sin(M-M1),
            + 0.0004 * sin(2*F+M),
            - 0.0004 * sin(2*F-M),
            - 0.0006 * sin(2*F+M1),
            + 0.0010 * sin(2*F-M1),
            + 0.0005 * sin(M+2*M1),
    ]

    # Новолуние -- юлианский день
    new_moon = JD + sum(adjusments)

    # Новолуние -- секунды с 1 января 1970
    seconds = (new_moon - start_of_epoch) * 86400

    # Високосные секунды -- приблизительно
    leap_seconds = (guess - 1972) / 1.5 + 42

    # Значение для пересчёта в местное время
    return seconds - leap_seconds

seconds_in_year = 31556952 # Григорианское приближение
seconds_in_month = 2551443 # Синодический месяц

about = 0
while 1:
    guess = about / seconds_in_year + 1970
    if guess > 2038: break
    unixtime = new_moon(guess)
    print time.localtime(unixtime)
    about += seconds_in_month
