Calculer les dates du calendrier correspondant aux phases de la lune:

• nouvelle lune
• premier quartier
• pleine lune
• dernier quartier

J'ai respecté le plus fidèlement possible l'excellent livre “Calculs astronomiques à l'usage des amateurs” de Jean Meeus, édité par la Société Astronomique de France (3 rue Beethoven 75016 PARIS). Ce petit livre peut être acheté par correspondance auprès de cet éditeur à un prix très raisonnable (env. 13 Euros).

Merci aussi à Christian Guine qui m'a mis le pied à l'étrier en me communicant son propre code.

Voici le code complet.

Ce code est “autoporteur”: vous pouvez le copier dans un fichier et le lancer dans une console ou dans IDLE. Et, bien sûr, changer les appels aux fonctions pour obtenir d'autres résultats.

Il est par ailleurs assez “auto-documenté”.

#!/usr/bin/python
# -*- coding: utf-8 -*-
 
from __future__ import division
 
from math import *
 
######################################################################
def bissextile(an):
    """dit si l'année donnée est bissextile ou non (True=oui, False=non)"""
    if (an % 4)==0:
        if ((an % 100)==0) and ((an % 400)<>0):
            return False
        else:
            return True
    else:
        return False
 
######################################################################
def postdate(D1,D2):
    """dit si une date D2 'j/m/a' est postérieure ou égale à une autre date D1 'j/m/a'"""
    (J1, M1, A1) = D1.split('/')
    (J2, M2, A2) = D2.split('/')
    if int(A2)>int(A1):
        return True
    elif int(A2)==int(A1):
        if int(M2)>int(M1):
            return True
        elif int(M2)==int(M1):
            if int(J2)>=int(J1):
                return True
    return False
 
######################################################################
def jj2date(JJ):
    """calcul d'une date (J,M,A) à partir du jour julien des éphémérides"""
    JJ += 0.5
    Z = int(JJ)
    F = JJ - Z
    if Z < 2299161:
        A = Z
    else:
        alpha = int((Z - 1867216.25) / 36524.25)
        A = Z + 1 + alpha - int(alpha / 4)
    B = A + 1524
    C = int((B - 122.1) / 365.25)
    D = int(365.25 * C)
    E = int((B - D) /30.6001)
 
    JD = B - D - int(30.6001 * E) + F          # calcul du jour décimal
 
    J = int(JD)                                # calcul du jour
 
    if (E < 13.5):                             # calcul du mois
        M = E - 1
    else:
        M = E - 13
 
    if (M > 2.5):                              # calcul de l'année
        A = C - 4716
    else:
        A = C - 4715
 
    return "%02d/%02d/%04d" % (J, M, A)
 
######################################################################
def calculphaseslune(k):
    """calcul de la date de la phase de la lune correspondant à la valeur k"""
 
    # calcul de T en fonction de k (formule 26.3 du livre de Jean MEEUS)
    T = k / 1236.85  
 
    # calcul de JJ = phase moyenne de la lune (formule 26.1 du livre de Jean MEEUS)
    JJ = 2415020.75933 + 29.53058868*k  + 0.0001178*T*T - 0.000000155*T*T*T   \
        + 0.00033 * sin(radians(166.56) + radians(132.87)*T - radians(0.009173)*T*T)
 
    # anomalie moyenne du soleil à l'instant JJ
    M = 359.2242 + 29.10535608*k - 0.0000333*T*T - 0.00000347*T*T*T
 
    # anomalie moyenne de la lune  à l'instant JJ
    MP = 306.0253 + 385.81691806*k + 0.0107306*T*T + 0.00001236*T*T*T
 
    # argument de la latitude de la lune  à l'instant JJ
    F = 21.2964 + 390.67050646*k - 0.0016528*T*T - 0.00000239*T*T*T
 
    # application de corrections selon les décimales de k
    kp = int(round((k - int(k))*100,0))
    if (kp==0) or (kp==50):
        # => nouvelle lune ou pleine lune
        JJ += (0.1734 - 0.000393*T) * sin(radians(M))  \
                + 0.0021 * sin(radians(2*M)) \
                - 0.4068 * sin(radians(MP))  \
                + 0.0161 * sin(radians(2*MP))  \
                - 0.0004 * sin(radians(3*MP))  \
                + 0.0104 * sin(radians(2*F))  \
                - 0.0051 * sin(radians(M+MP))  \
                - 0.0074 * sin(radians(M-MP))  \
                + 0.0004 * sin(radians(2*F+M))  \
                - 0.0004 * sin(radians(2*F-M))  \
                - 0.0006 * sin(radians(2*F+MP))  \
                + 0.0010 * sin(radians(2*F-MP))  \
                + 0.0005 * sin(radians(M+2*MP))
    else:
        JJ += (0.1721 - 0.0004*T) * sin(radians(M))  \
                + 0.0021 * sin(radians(2*M))  \
                - 0.6280 * sin(radians(MP))  \
                + 0.0089 * sin(radians(2*MP))  \
                - 0.0004 * sin(radians(3*MP))  \
                + 0.0079 * sin(radians(2*F))  \
                - 0.0119 * sin(radians(M+MP))  \
                - 0.0047 * sin(radians(M-MP))  \
                + 0.0003 * sin(radians(2*F+M))  \
                - 0.0004 * sin(radians(2*F-M))  \
                - 0.0006 * sin(radians(2*F+MP))  \
                + 0.0021 * sin(radians(2*F-MP))  \
                + 0.0003 * sin(radians(M+2*MP))  \
                + 0.0004 * sin(radians(M-2*MP))  \
                - 0.0003 * sin(radians(2*M+MP))
        if kp==25:
            # => premier quartier
            JJ += 0.0028 - 0.0004*cos(radians(M)) + 0.0003*cos(radians(MP))
        else:
            # => dernier quartier
            JJ += 0.0028 + 0.0004*cos(radians(M)) - 0.0003*cos(radians(MP))
 
    # retour de la date du calendrier grégorien à partir du jour Julien des éphémérides
    return jj2date(JJ)
 
######################################################################
def _phaseslune(D):
    """utilitaire de calcul de k pour la date D 'j/m/a' (phases de la lune)"""
 
    # extraction de la date
    lst=D.split('/')
    J = int(lst[0])
    M = int(lst[1])
    A = int(lst[2])
 
    # calcul de l'année décimale
    if bissextile(A):
        A += ((0,31,60,91,121,152,182,213,244,274,305,335,366)[M-1] + J) / 366
    else:
        A += ((0,31,59,90,120,151,181,212,243,273,304,334,365)[M-1] + J) / 365
 
    # renvoi de k
    return (A-1900)*12.3685
 
######################################################################
def phaseslune(D, n=1):
    """calcul des n dates de phase de la lune qui commence(nt) à la date D 'j/m/a' """
 
    # calcul du k initial basé sur D
    k=_phaseslune(D)
 
    # calcul de la 1ère date (éventuellement antérieure à D)
    k=int(k)
    p=0
    d = [p, calculphaseslune(k)]
 
    # trouver la 1ère date de phase lunaire postérieure ou égale à D
    while not postdate(D, d[1]):
        k+=0.25
        p = (p+1) % 4
        d = [p, calculphaseslune(k)]
    L=[d]
 
    # si n>1, calculer autant de dates de phase lunaire que demandé
    for i in xrange(0,n):
        k+=0.25
        p = (p+1) % 4
        d = [p, calculphaseslune(k)]
        L.append(d)
 
    # retour de la liste des dates cherchées avec le type de phase pour chaque date
    return L
 
######################################################################
def phaseslune2(D1, D2):
    """calcul des dates de phase de la lune à partir de D1 'j/m/a' et jusqu'à D2 'j/m/a' exclue"""
 
    # calcul du k initial basé sur D
    k=_phaseslune(D1)
 
    # calcul de la 1ère date (éventuellement antérieure à D)
    k=int(k)
    p=0
    d = [p, calculphaseslune(k)]
 
    # trouver la 1ère date de phase lunaire postérieure ou égale à D
    while not postdate(D1, d[1]):
        k+=0.25
        p = (p+1) % 4
        d = [p, calculphaseslune(k)]
    L=[d]
 
    # calculer autant de dates de phase lunaire qu'il y a jusqu'à D2 exclue
    while True:
        k+=0.25
        p = (p+1) % 4
        d = [p, calculphaseslune(k)]
        if postdate(D2, d[1]):
            break
        else:
            L.append(d)
 
    # retour de la liste des dates cherchées avec le type de phase pour chaque date
    return L
 
######################################################################
 
# exemples d'utilisation
 
# calculer la date à partir du jour Julien des Ephémérides
print jj2date(2436116.31)  # affiche: 04/10/1957
print jj2date(1842713.0)  # affiche: 27/01/0333
print jj2date(2443259.9)  # affiche: 26/04/1977
print "----------------------------------------------------------------------"
 
# calculer la date de la 1ère phase lunaire postérieure ou égale à la date donnée
print phaseslune('1/8/2007')  # affiche: [[3, '05/08/2007']]
print phaseslune('15/11/2007')  # affiche: [[1, '17/11/2007']]
print phaseslune('25/7/2008')  # affiche: [[3, '25/07/2008']]
print "----------------------------------------------------------------------"
 
# calculer les n=10 dates des phases lumaires commençant à la date D donnée
d='1/1/2008'
n=10
print d, n
L=phaseslune(d,n)
for i in xrange(0,n):
    print L[i]
# affiche:
# [0, '08/01/2008']
# [1, '15/01/2008']
# [2, '22/01/2008']
# [3, '30/01/2008']
# [0, '07/02/2008']
# [1, '14/02/2008']
# [2, '21/02/2008']
# [3, '29/02/2008']
# [0, '07/03/2008']
# [1, '14/03/2008']
print "----------------------------------------------------------------------"
 
# calculer les dates des phases lumaires commençant à la date D1 donnée et antérieures à D2
D1='1/7/2008'
D2='1/8/2008'
print D1,D2
L=phaseslune2(D1,D2)
for i in xrange(0,len(L)):
    print L[i]
# affiche: 
# [0, '03/07/2008']
# [1, '10/07/2008']
# [2, '18/07/2008']
# [3, '25/07/2008']
print
# cas d'un intervalle trop faible:
print phaseslune2('1/7/2009', '2/7/2009') # affiche: [[]] car aucune date de convient
print "----------------------------------------------------------------------"
 
# affichage plus facile à lire:
p=["nouvelle lune","premier quartier","pleine lune","dernier quartier"]
 
x=phaseslune('25/7/2008')
print x[0][1], ' : ', p[x[0][0]]  # affiche: "25/07/2008  :  dernier quartier"
 
x=phaseslune('5/3/2009')
print x[0][1], ' : ', p[x[0][0]]  # affiche: 11/03/2009  :  pleine lune
 
x=phaseslune('1/7/2009')
print x[0][1], ' : ', p[x[0][0]]  # affiche: 07/07/2009  :  pleine lune