Alors que j'attends encore des pièces pour terminer mes horloges, j'ai travaillé, assez intensément, sur mon prochain projet.
Comme il y a une partie très importante des calculs développée en PB, je peux en parler ici.
Il s'agira d'un lustre / planétarium centrée sur la Terre qui indiquera la position horizontale et la hauteur de tous les astres du système solaire, y compris, la Lune et le Soleil vues depuis votre position exacte sur Terre. Il faudra donc les indiquer (lignes 65 à 67) dans le programme de calcul que je vais vous fournir.
Global MyLat.d = xx.xxxxx ; en décimales
Global MyLon.d = xx.xxxx ; en décimales
Global MyAlt.d = xx.x ; en m
Commencer par cela car, sinon, le programme ne fonctionnera pas.
Pour le contrôle, j'ai fait un script pour Stellarium qui calcule les mêmes données que le programme permettant ainsi de vérifier les calculs nombreux et variés.
Ce programme sera convertit en une bibliothèque pour Arduino pour gérer les indications azimutales et de hauteur en degrés.
Voici ce programme :
Code : Tout sélectionner
; ---------------- Calculs lunaires, solaires et planétaires pour déterminer ----------------
; l'azimut et la hauteur de chaque astre du système solaire
; Programme de contrôle : Stellarium (v25.4)
; Version 1.00 PhM - février 2026 - PB 6.30 (x64)
; -------------------------------------------------------------------------------------------
EnableExplicit
;----- Structures
Structure SkyCoord
RA.d
Dec.d
Azimuth.d
Altitude.d
Dist.d
Phase.d
Age.d
PhaseName.s
UTC_Time.s
Local_Time.s
SunRA.d
SunDec.d
SunAz.d
SunAlt.d
SunEclX.d
SunEclY.d
; Planètes
MercRA.d
MercDec.d
MercAz.d
MercAlt.d
VenRA.d
VenDec.d
VenAz.d
VenAlt.d
MarsRA.d
MarsDec.d
MarsAz.d
MarsAlt.d
JupRA.d
JupDec.d
JupAz.d
JupAlt.d
SatRA.d
SatDec.d
SatAz.d
SatAlt.d
UraRA.d
UraDec.d
UraAz.d
UraAlt.d
NepRA.d
NepDec.d
NepAz.d
NepAlt.d
EndStructure
Structure Nutation
dPsi.d
dEps.d
Eps.d
EndStructure
;- --- Constantes et Globales
Global MyLat.d = 0
Global MyLon.d = 0
Global MyAlt.d = 0
#RAD = #PI / 180.0
#DEG = 180.0 / #PI
Global Sky.SkyCoord
;- --- Fonctions Utilitaires
Procedure.d Wrap360(a.d)
Protected res.d = a - 360.0 * Int(a / 360.0)
If res < 0
res + 360.0
EndIf
ProcedureReturn res
EndProcedure
Procedure.d AddRefraction(Alt.d)
If Alt < -0.575
ProcedureReturn Alt
EndIf
Protected R.d = 1.02 / Tan((Alt + 10.3 / (Alt + 5.11)) * #RAD)
ProcedureReturn Alt + (R / 60.0)
EndProcedure
Procedure GetHorizontal(RA.d, Dec.d, LST.d, *Az.Double, *Alt.Double)
Protected H.d = Wrap360(LST - RA) * #RAD
Protected Lat.d = MyLat * #RAD
Protected d.d = Dec * #RAD
Protected SinAlt.d = Sin(Lat) * Sin(d) + Cos(Lat) * Cos(d) * Cos(H)
*Alt\d = ASin(SinAlt) * #DEG
Protected Y.d = Sin(H)
Protected X.d = Cos(H) * Sin(Lat) - Tan(d) * Cos(Lat)
*Az\d = Wrap360(ATan2(X, Y) * #DEG + 180.0)
EndProcedure
Procedure.d CorrectParallax(Alt.d, DistanceUA.d)
If DistanceUA < 0.00001
ProcedureReturn Alt
EndIf
Protected HP.d = ASin(6378.14 / (DistanceUA * 149597870.7)) * #DEG
ProcedureReturn Alt - (HP * Cos(Alt * #RAD))
EndProcedure
Procedure GetNutation(T.d, *N.Nutation)
Protected Omega.d = (125.04452 - 1934.136261 * T) * #RAD
Protected L.d = (280.4665 + 36000.7698 * T) * #RAD
*N\dPsi = (-17.20 / 3600.0) * Sin(Omega) - (1.32 / 3600.0) * Sin(2 * L)
*N\dEps = (9.20 / 3600.0) * Cos(Omega) + (0.57 / 3600.0) * Cos(2 * L)
Protected Eps0.d = 23.43929111 - (46.8150 / 3600.0) * T
*N\Eps = (Eps0 + *N\dEps) * #RAD
EndProcedure
;- --- Moteurs de Calcul
Procedure GetSunPosition(T_Cent.d, LST_Deg.d, *Out.SkyCoord, TrueObl.d)
Protected L.d = Wrap360(280.46607 + 36000.7698 * T_Cent)
Protected M.d = Wrap360(357.52911 + 35999.0503 * T_Cent) * #RAD
Protected Lambda.d = Wrap360(L + 1.9146 * Sin(M) + 0.02 * Sin(2 * M)) * #RAD
*Out\SunEclX = Cos(Lambda)
*Out\SunEclY = Sin(Lambda)
*Out\SunRA = Wrap360(ATan2(Cos(Lambda), Sin(Lambda) * Cos(TrueObl)) * #DEG)
*Out\SunDec = ASin(Sin(Lambda) * Sin(TrueObl)) * #DEG
Protected Az.Double, Alt.Double
GetHorizontal(*Out\SunRA, *Out\SunDec, LST_Deg, @Az, @Alt)
*Out\SunAz = Az\d
*Out\SunAlt = Alt\d
EndProcedure
Procedure GetPlanet(T_Cent.d, LST_Deg.d, Planete.s, *Out.SkyCoord, TrueObl.d)
Protected a.d, e.d, i.d, Node.d, peri.d, M.d, E_anom.d, v.d, r.d, xG.d, yG.d, zG.d, RA.d, Dec.d, n.i
Protected DistGeo.d, Az.Double, Alt.Double, AltTopo.d
Select Planete
Case "Mercure"
a = 0.387098
e = 0.205630
i = 7.0047 * #RAD
Node = Wrap360(48.331 + 1.186 * T_Cent) * #RAD
peri = Wrap360(77.456 + 1.556 * T_Cent) * #RAD
M = Wrap360(174.796 + 149472.674 * T_Cent) * #RAD
Case "Venus"
a = 0.723330
e = 0.006773
i = 3.3946 * #RAD
Node = Wrap360(76.680 + 0.586 * T_Cent) * #RAD
peri = Wrap360(131.533 + 0.659 * T_Cent) * #RAD
M = Wrap360(50.447 + 58517.815 * T_Cent) * #RAD
Case "Mars"
a = 1.523688
e = 0.093405
i = 1.8497 * #RAD
Node = Wrap360(49.557 + 0.772 * T_Cent) * #RAD
peri = Wrap360(336.060 + 1.060 * T_Cent) * #RAD
M = Wrap360(19.387 + 19140.303 * T_Cent) * #RAD
Case "Jupiter"
a = 5.202597
e = 0.048464
i = 1.3031 * #RAD
Node = Wrap360(100.464 + 1.021 * T_Cent) * #RAD
peri = Wrap360(14.331 + 1.613 * T_Cent) * #RAD
M = Wrap360(20.020 + 3034.906 * T_Cent) * #RAD
Case "Saturne"
a = 9.541498
e = 0.054444
i = 2.4844 * #RAD
Node = Wrap360(113.665 + 1.395 * T_Cent) * #RAD
peri = Wrap360(92.598 + 1.964 * T_Cent) * #RAD
M = Wrap360(317.373 + 1222.113 * T_Cent) * #RAD
Case "Uranus"
a = 19.187979
e = 0.047193
i = 0.7733 * #RAD
Node = Wrap360(74.006 + 0.498 * T_Cent) * #RAD
peri = Wrap360(170.964 + 1.485 * T_Cent) * #RAD
M = Wrap360(142.951 + 428.482 * T_Cent) * #RAD
Case "Neptune"
a = 30.069033
e = 0.008590
i = 1.7691 * #RAD
Node = Wrap360(131.784 + 1.102 * T_Cent) * #RAD
peri = Wrap360(44.971 + 0.198 * T_Cent) * #RAD
M = Wrap360(260.247 + 218.459 * T_Cent) * #RAD
EndSelect
E_anom = M
For n = 1 To 5
E_anom = M + e * Sin(E_anom)
Next
v = 2.0 * ATan(Sqr((1.0+e)/(1.0-e)) * Tan(E_anom/2.0))
r = a * (1.0 - e * Cos(E_anom))
Protected u.d = v + (peri - Node)
Protected xH.d = r * (Cos(Node) * Cos(u) - Sin(Node) * Sin(u) * Cos(i))
Protected yH.d = r * (Sin(Node) * Cos(u) + Cos(Node) * Sin(u) * Cos(i))
Protected zH.d = r * (Sin(u) * Sin(i))
xG = xH + *Out\SunEclX
yG = yH + *Out\SunEclY
zG = zH
Protected LS.d = Wrap360(280.466 + 36000.77 * T_Cent) * #RAD
Protected ConstAb.d = (20.495 / 3600.0)
xG - ConstAb * Cos(LS)
yG - ConstAb * Sin(LS)
DistGeo = Sqr(xG*xG + yG*yG + zG*zG)
RA = Wrap360(ATan2(xG, yG * Cos(TrueObl) - zG * Sin(TrueObl)) * #DEG)
Dec = ASin((yG * Sin(TrueObl) + zG * Cos(TrueObl)) / DistGeo) * #DEG
GetHorizontal(RA, Dec, LST_Deg, @Az, @Alt)
AltTopo = CorrectParallax(Alt\d, DistGeo)
Select Planete
Case "Mercure"
*Out\MercRA = RA
*Out\MercDec = Dec
*Out\MercAz = Az\d
*Out\MercAlt = AddRefraction(AltTopo)
Case "Venus"
*Out\VenRA = RA
*Out\VenDec = Dec
*Out\VenAz = Az\d
*Out\VenAlt = AddRefraction(AltTopo)
Case "Mars"
*Out\MarsRA = RA
*Out\MarsDec = Dec
*Out\MarsAz = Az\d
*Out\MarsAlt = AddRefraction(AltTopo)
Case "Jupiter"
*Out\JupRA = RA
*Out\JupDec = Dec
*Out\JupAz = Az\d
*Out\JupAlt = AddRefraction(AltTopo)
Case "Saturne"
*Out\SatRA = RA
*Out\SatDec = Dec
*Out\SatAz = Az\d
*Out\SatAlt = AddRefraction(AltTopo)
Case "Uranus"
*Out\UraRA = RA
*Out\UraDec = Dec
*Out\UraAz = Az\d
*Out\UraAlt = AddRefraction(AltTopo)
Case "Neptune"
*Out\NepRA = RA
*Out\NepDec = Dec
*Out\NepAz = Az\d
*Out\NepAlt = AddRefraction(AltTopo)
EndSelect
EndProcedure
Procedure.d GetCurrentJD(*Out.SkyCoord = 0)
Protected T_Now = Date()
Protected vYear = Year(T_Now)
Protected vMonth = Month(T_Now)
Protected vDay = Day(T_Now)
Protected vHour = Hour(T_Now)
Protected Offset = 1
Protected LastSunMarch = 31 - DayOfWeek(Date(vYear, 3, 31, 2, 0, 0))
Protected LastSunOct = 31 - DayOfWeek(Date(vYear, 10, 31, 3, 0, 0))
If (vMonth > 3 And vMonth < 10) Or (vMonth = 3 And (vDay > LastSunMarch Or (vDay = LastSunMarch And vHour >= 2))) Or (vMonth = 10 And (vDay < LastSunOct Or (vDay = LastSunOct And vHour < 3)))
Offset = 2
EndIf
If *Out
*Out\Local_Time = FormatDate("%dd/%mm/%yyyy %hh:%ii:%ss", T_Now)
*Out\UTC_Time = FormatDate("%dd/%mm/%yyyy %hh:%ii:%ss", AddDate(T_Now, #PB_Date_Hour, -Offset))
EndIf
Protected h_dec.d = (vHour - Offset) + Minute(T_Now) / 60.0 + Second(T_Now) / 3600.0
Protected y_calc = vYear
Protected m_calc = vMonth
If m_calc <= 2
y_calc - 1
m_calc + 12
EndIf
Protected B_calc = 2 - Int(y_calc / 100) + Int(Int(y_calc / 100) / 4)
ProcedureReturn Int(365.25 * (y_calc + 4716)) + Int(30.6001 * (m_calc + 1)) + vDay + B_calc - 1524.5 + (h_dec + 69.2 / 3600.0) / 24.0
EndProcedure
Procedure.s PhaseLune(elongation.d)
If elongation < 5.0 Or elongation > 355.0
ProcedureReturn "Nouvelle lune"
EndIf
If Abs(elongation - 180.0) < 8.0
ProcedureReturn "Pleine lune"
EndIf
If Abs(elongation - 90.0) < 2.0
ProcedureReturn "Premier quartier"
EndIf
If Abs(elongation - 270.0) < 2.0
ProcedureReturn "Dernier quartier"
EndIf
If elongation < 180.0
If elongation < 90.0
ProcedureReturn "Premier croissant"
Else
ProcedureReturn "Gibbeuse croissante"
EndIf
Else
If elongation < 270.0
ProcedureReturn "Gibbeuse décroissante"
Else
ProcedureReturn "Dernier croissant"
EndIf
EndIf
EndProcedure
Procedure Update()
Shared Sky
Protected Az.Double
Protected Alt.Double
Protected JD_Now.d = GetCurrentJD(@Sky)
Protected T.d = (JD_Now - 2451545.0) / 36525.0
Protected LST.d = Wrap360(280.460618 + 360.985647 * (JD_Now - 2451545.0) + MyLon)
Protected Nut.Nutation
GetNutation(T, @Nut)
; --- SOLEIL
GetSunPosition(T, LST, @Sky, Nut\Eps)
; --- ARGUMENTS LUNAIRES
Protected L_Lune.d = Wrap360(218.316 + 481267.881 * T)
Protected M_Lune.d = Wrap360(134.963 + 477198.867 * T) * #RAD
Protected F_Lune.d = Wrap360(93.272 + 483202.017 * T) * #RAD
Protected D_Lune.d = Wrap360(297.850 + 445267.111 * T) * #RAD
Protected M_Sol.d = Wrap360(357.529 + 35999.050 * T) * #RAD
; --- LONGITUDE / LATITUDE LUNE
Protected Long.d = L_Lune + 6.289 * Sin(M_Lune) - 1.274 * Sin(M_Lune - 2 * D_Lune) + 0.658 * Sin(2 * D_Lune) - 0.185 * Sin(M_Sol)
Protected Latit.d = (5.128 * Sin(F_Lune)) * #RAD
Protected LongNut.d = (Long + Nut\dPsi) * #RAD
; --- COORDONNÉES ÉQUATORIALES LUNE
Sky\RA = Wrap360(ATan2(Cos(LongNut), Sin(LongNut) * Cos(Nut\Eps) - Tan(Latit) * Sin(Nut\Eps)) * #DEG)
Sky\Dec = ASin(Sin(Latit) * Cos(Nut\Eps) + Cos(Latit) * Sin(Nut\Eps) * Sin(LongNut)) * #DEG
; --- DISTANCE LUNE PRÉCISE (Termes de Brown/Chapront)
Protected DistKm.d
DistKm = 385001
DistKm - 20905 * Cos(M_Lune)
DistKm - 3699 * Cos(2 * D_Lune - M_Lune)
DistKm - 2956 * Cos(2 * D_Lune)
DistKm - 570 * Cos(M_Lune + 2 * D_Lune)
Sky\Dist = DistKm / 149597870.7
; --- COORDONNÉES HORIZONTALES LUNE
GetHorizontal(Sky\RA, Sky\Dec, LST, @Az, @Alt)
Sky\Azimuth = Az\d
Sky\Altitude = AddRefraction(CorrectParallax(Alt\d, Sky\Dist))
; --- PHASE ET AGE
Protected Elong.d = Wrap360(Long - (280.466 + 36000.77 * T))
Sky\Phase = (1 - Cos(Elong * #RAD)) / 2 * 100
Sky\Age = Wrap360(Elong) * 29.53 / 360.0
Sky\PhaseName = PhaseLune(Elong)
; --- CALCUL DES PLANÈTES
GetPlanet(T, LST, "Mercure", @Sky, Nut\Eps)
GetPlanet(T, LST, "Venus", @Sky, Nut\Eps)
GetPlanet(T, LST, "Mars", @Sky, Nut\Eps)
GetPlanet(T, LST, "Jupiter", @Sky, Nut\Eps)
GetPlanet(T, LST, "Saturne", @Sky, Nut\Eps)
GetPlanet(T, LST, "Uranus", @Sky, Nut\Eps)
GetPlanet(T, LST, "Neptune", @Sky, Nut\Eps)
; --- MISE À JOUR DES GADGETS
SetGadgetText(1, StrD(Sky\RA, 4) + "° / " + StrD(Sky\Dec, 4) + "°")
SetGadgetText(3, StrD(Sky\Azimuth, 4) + "° / " + StrD(Sky\Altitude, 4) + "°")
; Affichage Distance UA et Km (Gadget 5)
SetGadgetText(5, StrD(Sky\Dist, 8) + " UA / " + StrD(DistKm, 0) + " Km")
SetGadgetText(6, StrD(Sky\Phase, 2) + " %")
SetGadgetText(7, Sky\PhaseName + " (" + StrD(Sky\Age, 2) + " jours)")
SetGadgetText(8, Sky\Local_Time + " / " + Sky\UTC_Time)
SetGadgetText(9, StrD(Sky\SunAz, 4) + "° / " + StrD(Sky\SunAlt, 4) + "°")
; Mercure
SetGadgetText(20, StrD(Sky\MercRA, 4) + "° / " + StrD(Sky\MercDec, 4) + "°")
SetGadgetText(23, StrD(Sky\MercAz, 4) + "° / " + StrD(Sky\MercAlt, 4) + "°")
; Vénus
SetGadgetText(25, StrD(Sky\VenRA, 4) + "° / " + StrD(Sky\VenDec, 4) + "°")
SetGadgetText(22, StrD(Sky\VenAz, 4) + "° / " + StrD(Sky\VenAlt, 4) + "°")
; Mars
SetGadgetText(26, StrD(Sky\MarsRA, 4) + "° / " + StrD(Sky\MarsDec, 4) + "°")
SetGadgetText(27, StrD(Sky\MarsAz, 4) + "° / " + StrD(Sky\MarsAlt, 4) + "°")
; Jupiter
SetGadgetText(28, StrD(Sky\JupRA, 4) + "° / " + StrD(Sky\JupDec, 4) + "°")
SetGadgetText(29, StrD(Sky\JupAz, 4) + "° / " + StrD(Sky\JupAlt, 4) + "°")
; Saturne
SetGadgetText(35, StrD(Sky\SatRA, 4) + "° / " + StrD(Sky\SatDec, 4) + "°")
SetGadgetText(36, StrD(Sky\SatAz, 4) + "° / " + StrD(Sky\SatAlt, 4) + "°")
; Uranus
SetGadgetText(39, StrD(Sky\UraRA, 4) + "° / " + StrD(Sky\UraDec, 4) + "°")
SetGadgetText(40, StrD(Sky\UraAz, 4) + "° / " + StrD(Sky\UraAlt, 4) + "°")
; Neptune
SetGadgetText(43, StrD(Sky\NepRA, 4) + "° / " + StrD(Sky\NepDec, 4) + "°")
SetGadgetText(44, StrD(Sky\NepAz, 4) + "° / " + StrD(Sky\NepAlt, 4) + "°")
EndProcedure
;- --- Programme Principal
If OpenWindow(0, 0, 0, 440, 620, "Calculs Astronomiques du système solaire v1.00 - PhM", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
Define y = 20
TextGadget(10, 20, y, 160, 20, "Lune RA / Dec :")
TextGadget(1, 190, y, 240, 20, "")
y + 20
TextGadget(12, 20, y, 160, 20, "Lune Az / Alt :")
TextGadget(3, 190, y, 240, 20, "")
y + 20
TextGadget(14, 20, y, 160, 20, "Lune Distance :")
TextGadget(5, 190, y, 240, 20, "")
y + 20
TextGadget(15, 20, y, 160, 20, "Lune illumination :")
TextGadget(6, 190, y, 240, 20, "")
y + 20
TextGadget(16, 20, y, 160, 20, "Phase Lunaire :")
TextGadget(7, 190, y, 240, 20, "")
y + 40
TextGadget(18, 20, y, 160, 20, "Soleil Az / Alt :")
TextGadget(9, 190, y, 240, 20, "")
y + 40
TextGadget(24, 20, y, 160, 20, "Mercure RA / Dec :")
TextGadget(20, 190, y, 240, 20, "")
y + 20
TextGadget(19, 20, y, 160, 20, "Mercure Az / Alt :")
TextGadget(23, 190, y, 240, 20, "")
y + 35
TextGadget(30, 20, y, 160, 20, "Vénus RA / Dec :")
TextGadget(25, 190, y, 240, 20, "")
y + 20
TextGadget(21, 20, y, 160, 20, "Vénus Az / Alt :")
TextGadget(22, 190, y, 240, 20, "")
y + 35
TextGadget(31, 20, y, 160, 20, "Mars RA / Dec :")
TextGadget(26, 190, y, 240, 20, "")
y + 20
TextGadget(32, 20, y, 160, 20, "Mars Az / Alt :")
TextGadget(27, 190, y, 240, 20, "")
y + 35
TextGadget(33, 20, y, 160, 20, "Jupiter RA / Dec :")
TextGadget(28, 190, y, 240, 20, "")
y + 20
TextGadget(34, 20, y, 160, 20, "Jupiter Az / Alt :")
TextGadget(29, 190, y, 240, 20, "")
y + 35
TextGadget(37, 20, y, 160, 20, "Saturne RA / Dec :")
TextGadget(35, 190, y, 240, 20, "")
y + 20
TextGadget(38, 20, y, 160, 20, "Saturne Az / Alt :")
TextGadget(36, 190, y, 240, 20, "")
y + 35
TextGadget(41, 20, y, 160, 20, "Uranus RA / Dec :")
TextGadget(39, 190, y, 240, 20, "")
y + 20
TextGadget(42, 20, y, 160, 20, "Uranus Az / Alt :")
TextGadget(40, 190, y, 240, 20, "")
y + 35
TextGadget(45, 20, y, 160, 20, "Neptune RA / Dec :")
TextGadget(43, 190, y, 240, 20, "")
y + 20
TextGadget(46, 20, y, 160, 20, "Neptune Az / Alt :")
TextGadget(44, 190, y, 240, 20, "")
y + 35
TextGadget(17, 20, y, 160, 20, "Temps (Local / UTC) :")
TextGadget(8, 190, y, 240, 20, "")
Update()
AddWindowTimer(0, 1, 1000)
Repeat
Select WaitWindowEvent()
Case #PB_Event_Timer
Update()
Case #PB_Event_CloseWindow
End
EndSelect
ForEver
EndIf
