Calcul des phases lunaires pour la date du jour

[SPH]

ha, je savais pas.

[PhM]

Oui, je sais, je n'avance pas rapidement. Voilà ou j'en suis :



Mais, il y a une bonne raison, je me suis aperçu dans mes calculs et, notamment, pour l'illumination de la Lune, d'une lente dérive dans le temps. Pour étudier cela de plus près, j'ai relevé les valeurs sur Stellarium ainsi que sur mon programme sur 5 années (2025/2029), tous les jours de l'année à 0h00.

Je me suis donc remis au PureBasic, voici les résultats qui, sur 5 années révèlent une précision meilleure que +/- 2,5% :



Le programme journalier mis à jour :

Code : Tout sélectionner

; ----------- Calcul des phases lunaires pour la date du jour (Version optimisée) ---------------
;                         Programme de contrôle : Stellarium (v25.2)
;                     Version finale PhM - octobre 2025 - PB 6.21 (x64)
; -----------------------------------------------------------------------------------------------

EnableExplicit

; ------------------ Constantes
#DegToRad = #PI / 180
#RadToDeg = 180 / #PI
Global Lunaison.f = 29.530588853 ; Longueur moyenne du mois synodique

; ------------------ Structures de données
Structure LunarData
  elongation.f
  distance.f
  anomaly.f
EndStructure

; ------------------ Déclaration des procédures

Declare.f JourJulien(annee, mois, jour, heure, minute, seconde)
Declare.f AgeLune(JD.f, *data.LunarData, annee, mois, jour)
Declare.f IlluminationLune(JD.f, *data.LunarData, annee, mois, jour)
Declare.s PhaseLune(age.f, illumination.f)
Declare.f CalculUTCOffset(annee, mois, jour)
Declare.s SaisonHoraire(offset.f)
Declare.f PositifMod360(f.f)
; Nouvelle procédure pour calculer Soleil et Lune ensemble
Declare CalculateSolarLunar(JD.f, *lambdaS.Float, *lambdaL.Float, *M.Float, *Mp.Float, *L.Float, *D.Float)

; ------------------ Date et heure locale avec correction automatique heure d'été/hiver

Global an = Year(Date())
Global mo = Month(Date())
Global jo = Day(Date())

;Global an = 2025
;Global mo = 1
;Global jo = 1

Global UTCOffset.f = CalculUTCOffset(an, mo, jo)

Global heure = Hour(Date()) - UTCOffset
Global minute = Minute(Date())
Global seconde = Second(Date())

;Global heure = 0 - UTCOffset
;Global minute = 0
;Global seconde = 0

; ------------------ Boucle principale du programme ------------------

Define JD.f = JourJulien(an, mo, jo, heure, minute, seconde)

Define Lunar.LunarData
Define age.f = AgeLune(JD, @Lunar, an, mo, jo)
Define illumination.f = IlluminationLune(JD, @Lunar, an, mo, jo)
Define phase$ = PhaseLune(age, illumination)

; Affichage des résultats
MessageRequester("Observation lunaire du " + Str(jo) + "/" + Str(mo) + "/" + Str(an) + " (" + SaisonHoraire(UTCOffset) + ")",
  "Heure locale : " + Str(heure + UTCOffset) + "h" + Str(minute) + "mn" + #LF$ +
  "Âge lunaire : " + StrF(age, 2) + " jours" + #LF$ +
  "Illumination : " + StrF(illumination, 2) + " %" + #LF$ +
  "Phase : " + phase$)

End

; ------------------ Fonction utilitaire pour modulo 360 correct
Procedure.f PositifMod360(f.f)
  f = f - Int(f / 360) * 360
  If f >= 360 : f - 360 : EndIf
  If f < 0 : f + 360 : EndIf
  ProcedureReturn f
EndProcedure

; ------------------ Procédures de calculs

Procedure.f JourJulien(annee, mois, jour, heure, minute, seconde)
  Protected a, b, y, m
  y = annee
  m = mois

  If m <= 2
    y = annee - 1
    m = mois + 12
  EndIf

  a = y / 100
  b = 2 - a + a / 4
  Protected JD.f = Int(365.25 * (y + 4716)) + Int(30.6001 * (m + 1)) + jour + b - 1524.5
  JD = JD + (heure + minute / 60.0 + seconde / 3600.0) / 24.0

  ProcedureReturn JD
EndProcedure

; Procédure unique pour calculer Soleil et Lune
Procedure CalculateSolarLunar(JD.f, *lambdaS.Float, *lambdaL.Float, *M.Float, *Mp.Float, *L.Float, *D.Float)
  Protected T.f = (JD - 2451545.0) / 36525.0

  ; --- Soleil ---
  Protected M_val.f = 357.5291092 + 35999.0502909 * T - 0.0001536 * T * T + Pow(T, 3) / 24490000
  M_val = PositifMod360(M_val)
  *M\f = M_val

  Protected Ls.f = 280.46646 + 36000.76983 * T + 0.0003032 * T * T
  Ls = PositifMod360(Ls)

  *lambdaS\f = Ls + 1.914602 * Sin(M_val * #DegToRad) - 0.004817 * Sin(2 * M_val * #DegToRad) - 0.019993 * Sin(3 * M_val * #DegToRad)
  *lambdaS\f = PositifMod360(*lambdaS\f)

  ; --- Lune ---
  Protected L_val.f = 218.3164591 + 481267.88134236 * T - 0.0013268 * T * T + Pow(T, 3) / 538841 - Pow(T, 4) / 65194000
  Protected T3.f = T * T * T
  Protected T4.f = T3 * T
  L_val = L_val + T3 / 538841 - T4 / 65194000
  L_val = PositifMod360(L_val)
  *L\f = L_val

  Protected Mp_val.f = 134.9633964 + 477198.8675055 * T + 0.0087414 * T * T + Pow(T, 3) / 69699 - Pow(T, 4) / 147120000
  Mp_val = PositifMod360(Mp_val)
  *Mp\f = Mp_val

  Protected D_val.f = L_val - *lambdaS\f
  D_val = PositifMod360(D_val)
  *D\f = D_val

  ; --- Correction de longitude lunaire (ELP2000/82)
  Protected correctionLune.f = 6.289 * Sin(Mp_val * #DegToRad) - 1.274 * Sin((2*L_val - *lambdaS\f) * #DegToRad) + 0.658 * Sin(2*(L_val - *lambdaS\f) * #DegToRad) + 0.214 * Sin(2 * D_val * #DegToRad)
  Protected lambdaL_val.f = L_val + correctionLune
  lambdaL_val = PositifMod360(lambdaL_val)
  *lambdaL\f = lambdaL_val

  ; --- Longitude vraie de la Lune (ajout des termes périodiques)
  Protected term1.f = 0.00014 * Sin((2*D_val - Mp_val) * #DegToRad)
  Protected term2.f = 0.00012 * Sin((2*D_val + Mp_val) * #DegToRad)
  Protected term3.f = 0.00006 * Sin((2*D_val - 2*Mp_val) * #DegToRad)
  Protected term4.f = 0.00006 * Sin((2*D_val + 2*Mp_val) * #DegToRad)
  Protected term5.f = 0.00004 * Sin((2*D_val - 3*Mp_val) * #DegToRad)
  Protected term6.f = 0.00004 * Sin((2*D_val + 3*Mp_val) * #DegToRad)
  Protected term7.f = 0.00003 * Sin((2*D_val - 4*Mp_val) * #DegToRad)
  Protected term8.f = 0.00003 * Sin((2*D_val + 4*Mp_val) * #DegToRad)
  Protected term9.f = 0.00002 * Sin((2*D_val - 5*Mp_val) * #DegToRad)
  Protected term10.f = 0.00002 * Sin((2*D_val + 5*Mp_val) * #DegToRad)

  lambdaL_val = lambdaL_val + term1 + term2 + term3 + term4 + term5 + term6 + term7 + term8 + term9 + term10
  lambdaL_val = PositifMod360(lambdaL_val)
  *lambdaL\f = lambdaL_val
EndProcedure

Procedure.f AgeLune(JD.f, *data.LunarData, annee, mois, jour)
  Protected lambdaS.f, lambdaL.f, M.f, Mp.f, L.f, D.f
  CalculateSolarLunar(JD, @lambdaS, @lambdaL, @M, @Mp, @L, @D)

  Protected diff.f = lambdaL - lambdaS
  If diff < 0 : diff + 360 : EndIf

  Protected age.f = diff / 360.0 * Lunaison

  ; --- CORRECTION DYNAMIQUE ---
  Protected v_corr.f = *data\anomaly
  If v_corr > 180 : v_corr = 360 - v_corr : EndIf
  v_corr = (v_corr - 90) / 90.0

  Protected correction_dynamique.f = 0.05 * v_corr
  age = age + correction_dynamique

  ; --- Correction empirique linéaire ---
  Protected correction_empirique.f = 0.0

  ; Exemple : ajouter 0.01 jour par jour depuis le 1er octobre
  Protected jours_depuis_debut.f = (mois - 10) * 30 + (jour - 1)
  correction_empirique = 0.01 * jours_depuis_debut

  age = age + correction_empirique

  age = age - Int(age / Lunaison) * Lunaison
  If age < 0 : age + Lunaison : EndIf
  ProcedureReturn age
EndProcedure

Procedure.f IlluminationLune(JD.f, *data.LunarData, annee, mois, jour)
  Protected lambdaS.f, lambdaL.f, M.f, Mp.f, L.f, D.f
  CalculateSolarLunar(JD, @lambdaS, @lambdaL, @M, @Mp, @L, @D)

  Protected T.f = (JD - 2451545.0) / 36525.0

  ; --- Calcul des positions du Soleil ---
  Protected M_val.f = 357.5291092 + 35999.0502909 * T - 0.0001536 * T * T + Pow(T, 3) / 24490000
  M_val = PositifMod360(M_val)
  Protected Ls.f = 280.46646 + 36000.76983 * T + 0.0003032 * T * T
  Ls = PositifMod360(Ls)
  Protected lambdaS_val.f = Ls + 1.914602 * Sin(M_val * #DegToRad) - 0.004817 * Sin(2 * M_val * #DegToRad) - 0.019993 * Sin(3 * M_val * #DegToRad)
  lambdaS_val = PositifMod360(lambdaS_val)

  ; --- Calcul des positions de la Lune ---
  Protected L_val.f = 218.3164591 + 481267.88134236 * T - 0.0013268 * T * T + Pow(T, 3) / 538841 - Pow(T, 4) / 65194000
  Protected Mp_val.f = 134.9633964 + 477198.8675055 * T + 0.0087414 * T * T + Pow(T, 3) / 69699 - Pow(T, 4) / 147120000
  Mp_val = PositifMod360(Mp_val)
  Protected D_val.f = L_val - lambdaS_val
  D_val = PositifMod360(D_val)

  ; --- Correction de longitude lunaire ---
  Protected correctionLune.f = 6.289 * Sin(Mp_val * #DegToRad) - 1.274 * Sin((2*L_val - lambdaS_val) * #DegToRad) + 0.658 * Sin(2*(L_val - lambdaS_val) * #DegToRad) + 0.214 * Sin(2 * D_val * #DegToRad)
  Protected lambdaL_val.f = L_val + correctionLune
  lambdaL_val = PositifMod360(lambdaL_val)

  ; --- Argument de latitude F et latitude lunaire beta ---
  Protected F.f = 93.2720950 + 483202.0175233 * T - 0.0036539 * T * T - Pow(T, 3) / 3526000 + Pow(T, 4) / 863310000
  F = PositifMod360(F)
  Protected beta.f = 5.128 * Sin(F * #DegToRad)

  ; --- Anomalie vraie de la Lune ---
  Protected ecc.f = 0.0549006
  Protected M_rad.f = Mp_val * #DegToRad
  Protected E.f = M_rad + ecc * Sin(M_rad) + (ecc * ecc / 2.0) * Sin(2.0 * M_rad)
  Protected tanHalfE.f = Tan(E / 2.0)
  Protected v.f = 2.0 * ATan(Sqr((1.0 + ecc) / (1.0 - ecc)) * tanHalfE)
  If v < 0 : v + 2.0 * #PI : EndIf
  *data\anomaly = v * #RadToDeg

  ; --- Distance Terre-Lune ---
  Protected rL.f = 382700.0 * (1.0 - ecc * Cos(v))
  *data\distance = rL

  ; --- Coordonnées cartésiennes ---
  Protected lambdaS_rad.f = lambdaS_val * #DegToRad
  Protected betaS_rad.f = 0
  Protected rS.f = 149597870 * (1.00014 - 0.01671 * Cos(M_val * #DegToRad) - 0.00014 * Cos(2 * M_val * #DegToRad))
  Protected Xs.f = rS * Cos(lambdaS_rad) * Cos(betaS_rad)
  Protected Ys.f = rS * Sin(lambdaS_rad) * Cos(betaS_rad)
  Protected Zs.f = rS * Sin(betaS_rad)

  Protected lambdaL_rad.f = lambdaL_val * #DegToRad
  Protected betaL_rad.f = beta * #DegToRad
  Protected Xl.f = rL * Cos(lambdaL_rad) * Cos(betaL_rad)
  Protected Yl.f = rL * Sin(lambdaL_rad) * Cos(betaL_rad)
  Protected Zl.f = rL * Sin(betaL_rad)

  ; --- Vecteurs ---
  Protected Xls.f = Xs - Xl
  Protected Yls.f = Ys - Yl
  Protected Zls.f = Zs - Zl
  Protected Xle.f = -Xl
  Protected Yle.f = -Yl
  Protected Zle.f = -Zl

  ; --- Produits scalaires et angle de phase ---
  Protected dot_ls_le.f = Xls * Xle + Yls * Yle + Zls * Zle
  Protected mag_ls.f = Sqr(Xls * Xls + Yls * Yls + Zls * Zls)
  Protected mag_le.f = Sqr(Xle * Xle + Yle * Yle + Zle * Zle)
  Protected cos_i.f = dot_ls_le / (mag_ls * mag_le)
  If cos_i > 1 : cos_i = 1 : ElseIf cos_i < -1 : cos_i = -1 : EndIf
  Protected i.f = ACos(cos_i) * #RadToDeg

  ; --- Illumination de base ---
  Protected k.f = (1 + Cos(i * #DegToRad)) / 2 * 100

  ; --- Correction de parallaxe ---
  Protected parallax_correction.f = 0.0025 * Sin((lambdaL_val - lambdaS_val) * #DegToRad) * (384400 / *data\distance)
  i = i + parallax_correction * #RadToDeg

  ; --- Recalcul final de k après corrections ---
  k = (1 + Cos(i * #DegToRad)) / 2 * 100

  ; --- Compensation empirique adaptative de la courbe de la tendance pluri-annuelle ---
  Protected age_calc.f = (lambdaL_val - lambdaS_val) / 360.0 * Lunaison
  If age_calc < 0 : age_calc + Lunaison : EndIf
  k = k - (Pow(((an-2000)/ 90),3) + Sqr(age_calc/55) - Sqr(jour/500))

  If k < 0 : k = 0 : EndIf
  If k > 100 : k = 100 : EndIf

  ProcedureReturn k
EndProcedure

Procedure.s PhaseLune(age.f, illumination.f)

  ;Phase lunaire Icône 🌙 % d'illumination approximatif Description rapide
  ;Nouvelle lune 🌑 ~0 % Face visible entièrement dans l'ombre
  ;Premier croissant 🌒 ~1 % à 49 % Fine lumière croissante à droite
  ;Premier quartier 🌓 ~50 % Moitié droite éclairée
  ;Gibbeuse croissante 🌔 ~51 % à 99 % Presque pleine, en croissance
  ;Pleine lune 🌕 ~100 % Face totalement éclairée
  ;Gibbeuse décroissante 🌖 ~99 % à 51 % Décroît après la pleine lune
  ;Dernier quartier 🌗 ~50 % Moitié gauche éclairée
  ;Dernier croissant 🌘 ~49 % à 1 % Fine lumière décroissante à gauche

  ; --- Nouvelle lune ---
  If illumination < 0.5
    ProcedureReturn "🌑 Nouvelle lune"
  EndIf

  ; --- Premier croissant ---
  If illumination >= 0.5 And illumination < 50.0 And age < 7.4
    ProcedureReturn "🌒 Premier croissant"
  EndIf

  ; --- Premier quartier ---
  If illumination >= 48.0 And illumination <= 52.0 And age >= 7.4 And age <= 8.4
    ProcedureReturn "🌓 Premier quartier"
  EndIf

  ; --- Gibbeuse croissante ---
  If illumination > 52.0 And illumination < 99.5 And age > 8.4 And age < 15.8
    ProcedureReturn "🌔 Gibbeuse croissante"
  EndIf

  ; --- Pleine lune ---
  If illumination >= 99.5 And age >= 14.8 And age <= 15.8
    ProcedureReturn "🌕 Pleine lune"
  EndIf

  ; --- Gibbeuse décroissante ---
  If illumination > 52.0 And illumination < 99.5 And age > 15.8 And age < 22.2
    ProcedureReturn "🌖 Gibbeuse décroissante"
  EndIf

  ; --- Dernier quartier ---
  If illumination >= 48.0 And illumination <= 52.0 And age >= 22.2 And age <= 23.2
    ProcedureReturn "🌗 Dernier quartier"
  EndIf

  ; --- Dernier croissant ---
  If illumination < 50.0 And age > 23.2
    ProcedureReturn "🌘 Dernier croissant"
  EndIf

  ; --- Fallback basé sur l'âge si illumination ambigüe ---
  If age < 7.4
    ProcedureReturn "🌒 Premier croissant"
  ElseIf age < 8.4
    ProcedureReturn "🌓 Premier quartier"
  ElseIf age < 14.8
    ProcedureReturn "🌔 Gibbeuse croissante"
  ElseIf age < 15.8
    ProcedureReturn "🌕 Pleine lune"
  ElseIf age < 22.2
    ProcedureReturn "🌖 Gibbeuse décroissante"
  ElseIf age < 23.2
    ProcedureReturn "🌗 Dernier quartier"
  Else
    ProcedureReturn "🌘 Dernier croissant"
  EndIf

EndProcedure

Procedure.f CalculUTCOffset(annee, mois, jour)
  Protected dimancheMars, dimancheOctobre, i
  For i = 31 To 25 Step -1
    If DayOfWeek(Date(annee, 3, i, 0, 0, 0)) = 0
      dimancheMars = i
      Break
    EndIf
  Next
  For i = 31 To 25 Step -1
    If DayOfWeek(Date(annee, 10, i, 0, 0, 0)) = 0
      dimancheOctobre = i
      Break
    EndIf
  Next
  If (mois > 3 And mois < 10) Or (mois = 3 And jour >= dimancheMars) Or (mois = 10 And jour < dimancheOctobre)
    ProcedureReturn 2.0
  Else
    ProcedureReturn 1.0
  EndIf
EndProcedure

Procedure.s SaisonHoraire(offset.f)
  If offset = 2.0
    ProcedureReturn "Heure d'été : UTC+2"
  Else
    ProcedureReturn "Heure d'hiver : UTC+1"
  EndIf
EndProcedure

J'ai également refait la version sur x années avec les mêmes critères (sur demande, si vous souhaitiez l'obtenir).

J'y retourne...

[PhM]

Ce matin, grosse séance de câble avec 79/370 leds câblées, 22% de la tache !



A suivre...

[Jacobus]

Ahah tout ton bazar me rappelle ma jeunesse d'étudiant en électronique (nostalgie ).
Il fallait tout faire; tôlerie, câblage, peinture... même les circuits imprimés en cuivre sur bakélite où on devait tracer notre circuit avec des lignes et des pastilles en plastique noir avant de les tremper dans l'acide... tout un univers de geek du premier stade, de l'époque où on ne pouvait imprimer que des longues ramettes de données texte et ou les ordis n'affichaient eux aussi que du texte ou du code, en vert ou en orange suivant le type d'appareil... Heureusement que tout cela a changé. Ma première conception je dois l'avoir encore dans un coin, une alimentation en 12V de la taille d'une batterie de bagnole
En tout cas j'aime bien ton projet, très sympa.

[PhM]

Salut Jacobus,

Oui, c'est du bazar à première vue, mais tout est utile, tu t'en doutes. Le câblage des 370 leds n'étant pas la partie la plus excitante de la réalisation mais, il faut bien passé par là avant de reprendre le code de commande des leds pour en faire des mots via l'Arduino.

Pour la nostalgie d'un passé technologique, il suffit que je me retourne vers mes étagères pour y retourner directement (avec du très moderne aussi).





Tant mieux si tu aimes ce projet ainsi que quelques autres, au regard du compteur. Cela m'a beaucoup étonné étant donné que je m'éloigne de PureBasic en grande partie. A quand une version de PureBasic compatible avec le monde Arduino...

Je vais donc continuer, jusqu'à l'accrochage au mur de cette horloge, à vous tenir au courant de l'avancement de ce projet.J'ai également l'intention, en fin de projet, de faire un dossier de fabrication (zip complet : manuel d'instruction, pdf annexes, plans mécaniques, schéma électrique, codes PB et Arduino, bibliothèque LunaPhases) pour le mettre à disposition pour ceux qui souhaiteraient en savoir plus.

Voilà, cette semaine j'espère bien avancé sur le câblage des leds, A suivre donc...

[Jacobus]

Ne soit pas étonné que ça plaise, et c'est parfaitement lié à PB puisque il te sert à créer ton projet.
Ce qui plaît c'est que ton code débouche sur du concret, une chose réelle que tu peux toucher contrairement aux codes en général qui restent des données qui n'existent pas matériellement. Le fait de créer quelque chose de matériel reste une des grandes joies de la vie
Tu devrais le proposer à un industriel qui en ferait un objet utile et décoratif distribuable au plus grand nombre. Ça peut, par exemple, faire un chouette cadeau de Noël
Amuses toi bien

[PhM]

Progression en cours avec 146 leds câblées avec déjà, environ 4,5m de câble en série.



Voici le schéma de principe de l'ensemble de l'horloge. Il n'est pas très complexe car, c'est le code Arduino qui fait tout le travail.



A suivre...

[venom]

Ce qui plaît c'est que ton code débouche sur du concret, une chose réelle que tu peux toucher contrairement aux codes en général qui restent des données qui n'existent pas matériellement. Le fait de créer quelque chose de matériel reste une des grandes joies de la vie

Je suis bien d'accord avec toi jacobus. Quand j'ai fait mon premier projet/programme et que j'ai vu les composants se mettent en route et que tout fonctionnait comme écrit dans mon code....

Quel satisfaction que j'ai ressenti Du bonheur que je pouvais toucher

Bravo phM ça avance bien. Courage pour les LEDs un moment dit a passer






@++