Algoritmi

Menù Algoritmi

Calcolo del Giorno Giuliano

Unita la data nel formato AAAA.MMGGdd
Esempio 21 Marzo 2012 alle ore 12:30 = 2012.032152
Se MM > 2 si è preso a = AAAA e m = MM
Se MM = 1 oppure 2 si è preso a = AAAA - 1 e m = MM + 12
Se il numero AAAA.MMGGdd è >= 1582.1015
Calcolato
A = Parte Intera di (a / 100)
B = 2 - A + Parte Intera di (A / 4)
Se invece AAAA.MMGGdd è < 1582.1015
allora A = B = 0

JD = Parte Intera di [365.25 * (a + 4716)] + Parte Intera di [30.6001 * (m + 1)] + GG.dd + B - 1524.5

(Astronomia con il computer - Jean Meeus 1990)

Torna a inizio pagina

Calcolo del Tempo Sidereo apparente

Calcolato il JD corrispondente alla data alle ore 0

Secoli T a partire dal 1/1/2000
T = (JD0 - 2451545.0) / 36525

Tempo siderale medio a Greenwich a ore 0 UT
θº = 100.46061837 + 36000.770053608 * T + 0.000387933 * T² - T³ / 38710000
il risultato è espresso in gradi ed è stato convertito in ore decimali
θ in ore dec. = θº / 15

Tempo siderale medio locale
Ts medio = θ in ore dec. + 1.00273790935 * orario UT (in ore decimali)

Calcoli per la nutazione e l'obliquità dell'ecclittica

Calcolato il JD corrispondente alla data ed alle ore UT

Secoli T a partire dal 1/1/2000
T = (JD - 2451545.0) / 36525

Longitudine media del sole
L = 280.4665° + 36000.7698° * T + 0.000303° * T²
Longitudine media della luna
L1 = 218.3164° + 481267.8812° * T - 0.001599° * T²
Anomalia media del sole
M = 357.5291° + 35999.0503° * T - 0.000154° * T²
Anomalia media della luna
M1 = 134.9634° + 477198.8675° * T + 0.008721° * T²
Longitudine media del nodo ascendente dell'orbita lunare sull'ecclittica misurata a partire dall'equinozio medio della data
Ω = 125.0443° - 1934.1363° * T + 0.002075° * T²

Nutazione in longitudine il risultato è espresso in arcosecondi
Δψ = - (17.1996” + 0.01742” * T) * sen Ω
- (1.3187” + 0.00016” * T) * sen (2L)
- 0.2274” * sen (2L1)
+ 0.2062” * sen (2Ω)
+ (0.1426” - 0.00034” * T) * sen M
+ 0.0712” * sen M1
- (0.0517” - 0.00012” * T) * sen (2L + M)
- 0.0386” * sen (2L1 - Ω)
- 0.0301” * sen (2L1 + M1)
+ 0.0217” * sen (2L - M)
- 0.0158” * sen (2L - 2L1 + M1)
+ 0.0129” * sen (2L - Ω)
+ 0.0123” * sen (2L1 - M1)

Nutazione in obliquità il risultato è espresso in arcosecondi
Δε = + (9.2025 + 0.00089 * T) * cos Ω
+ (0.5736 - 0.00031 * T) * cos 2L
+ 0.0977 * cos 2L1
- 0.0895 cos 2Ω
+ 0.0224 * cos (2L + M)
+ 0.0200 * cos (2L1 - Ω)
+ 0.0129 * cos (2L1 + M)
- 0.0095 * cos (2L - M)
- 0.0070 * cos (2L - Ω)

Obliquità
ε0 = 23°26'21".448 - 46".8150 * T - 0".00059 * T² + 0".001813 * T³
ε = ε0 + Δε

Correzione in secondi = Δψ * cos ε / 15

Ts apparente = Ts medio + (Correzione / 3600)

(Calcules Astronomiques - Jean Meeus 2014)

Torna a inizio pagina

Calcolo delle coordinate solari
metodo con l'ausilio della teoria Bretagnon Francou VSOP-87

Il calcolo della posizione del sole è stato effettuato con l'estratto dei più importanti termini periodici della teoria Bretagnon Francou VSOP-87 (vedi procedura calcolo coordinate planetarie relative al pianeta terra)

ottenendo:
λT - Longitudine dell'eclittica terrestre
βT - Latitudine dell'eclittica terrestre
R - Raggio vettore sole-terra

λ - Longitudine geocentrica del sole = λT + 180
β - Latitudine geocentrica del sole = -βT

Conversione delle coordinate geocentriche eclittiche nel sistema FK5

Secoli Giuliani dall'epoca standard J2000 alla data inserita
T = (JD - 2451545.0) ÷ 36525

λ' = λ - 1°.397 * T - 0°.00031 * T²

Δλ = - 0".09033
Δβ = + 0".03916 * (cosλ' - senλ')
(entrambi i 2 precedenti risultati sono espressi in arcosecondi)

λ Sole convertita in FK5 = λ + Δλ/3600
β Sole convertita in FK5 = β + Δβ/3600

Calcolo dell'aberrazione (A)

in questo programma l'aberrazione è stata calcolata con la seguente formula:

Millenni Giuliani dall'epoca standard J2000 alla data inserita
ro = (JD - 2451545.0) ÷ 365250

ABB_Δλ = 3548.193
+ 118.568 * sin( 87.5287 + 359993.7286 * ro)
+ 2.476 * sin( 85.0561 + 719987.4571 * ro)
+ 1.376 * sin( 27.8502 + 4452671.1152 * ro)
+ 0.119 * sin( 73.1375 + 450368.8564 * ro)
+ 0.114 * sin(337.2264 + 329644.6718 * ro)
+ 0.086 * sin(222.5400 + 659289.3436 * ro)
+ 0.078 * sin(162.8136 + 9224659.7915 * ro)
+ 0.054 * sin( 82.5823 + 1079981.1857 * ro)
+ 0.052 * sin(171.5189 + 225184.4282 * ro)
+ 0.034 * sin( 30.3214 + 4092677.3866 * ro)
+ 0.033 * sin(119.8105 + 337181.4711 * ro)
+ 0.023 * sin(247.5418 + 299295.6151 * ro)
+ 0.023 * sin(325.1526 + 315559.5560 * ro)
+ 0.021 * sin(155.1241 + 675553.2846 * ro)
+ 7.311 * ro * sin(333.4515 + 359993.7286 * ro)
+ 0.305 * ro * sin(330.9814 + 719987.4571 * ro)
+ 0.010 * ro * sin(328.5170 + 1079981.1857 * ro)
+ 0.309 * ro² * sin(241.4518 + 359993.7286 * ro)
+ 0.021 * ro² * sin(205.0482 + 719987.4571 * ro)
+ 0.004 * ro² * sin(297.8610 + 4452671.1152 * ro)
+ 0.010 * ro² * sin(154.7066 + 359993.7286 * ro)
(gli argomenti del seno sono espressi in gradi decimali)

A = -0.005775518 * R * ABB_Δλ
(Aberrazione espressa in arcosecondi)
dove R è il raggio vettore sole - terra estratto dai termini periodici del Vsop87

L'aberrazione può essere calcolata anche più semplicemente con la formula che segue
ma anche se le differenze tra i due sistemi sono minime si ritiene più precisa la precedente

A = - 20".4898 / R
(Aberrazione espressa in arcosecondi)

Δψ - Nutazione in longitudine (Vedi calcolo tempo sidereo)
(espressa in arcosecondi)

λ apparente del sole corretta da aberrazione e nutazione = λ Sole convertita in FK5 + A / 3600 + Δψ / 3600
β apparente del sole = β Sole convertita in FK5 (rimane invariata)

Obliquità
ε0 = 23°26'21".448 - 46".8150 * T - 0".00059 * T² + 0".001813 * T³
ε = ε0 + Δε
per il Δε vedi "Calcolo del tempo sidereo apparente"

Declinazione del sole
δ = arcsen (sen β * cos ε + cos β * sen ε * sen λ)

Ascensione Retta del sole
N = sen λ * cos ε - tan β * sen ε
D = cos λ
tan α = N / D
Se N è positivo e D è positivo ==> α = α
Se N è positivo e D è negativo ==> α = α +180
Se N è negativo e D è positivo ==> α = α +360
Se N è negativo e D è negativo ==> α = α +180

Calcolo coordinate altazimutali
(vedi calcolo delle coordinate altazimutali)

Calcolo del sorgere, passagio al meridiano e tramonto del Sole
Vedi calcolo di sorgere, trasito e tramonto

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo delle coordinate solari
metodo senza l'utilizzo della teoria Bretagnon Francou VSOP-87
questo programma non usa questo sistema perché ha un accuratezza inferiore al sistema Vsop87
le formule che seguono sono solo informative
tuttavia il solo risultato Equazione del Tempo è calcolato con questo sistema

Secoli Giuliani dall'epoca standard J2000 alla data inserita
T = (JD - 2451545.0) ÷ 36525

Longitudine geometrica media del Sole
L = 280.46646° + 36000.76983° × (T) + 0.0003032° × (T)²

Anomalia media del Sole
M = 357.52911° + 35999.05029° × (T) - 0.0001537° × (T)²

Equazione del centro del Sole
C = [1.914602° - 0.004817° × (T) - 0.000014° × (T)²] × sen M + [0.019993° - 0.000101° × (T)] × sen (2 × M) + 0.000289° × sen (3 × M)

Longitudine vera del Sole
LV = L + C

Longitudine apparente del Sole
LA = LV - 0.00569° - 0.00478° × sen (125.04° - 1934.136° × T)

Calcolo dell' obliquità dell'ecclittica

Longitudine media del sole
L = 280.4665° + 36000.7698° + 0.000303° * T²
Longitudine media della luna
L1 = 218.3164° + 481267.8812° * T - 0.001599° * T²
Anomalia media del sole
M = 357.5291° + 35999.0503° * T - 0,000154° * T²
Anomalia media della luna
M1 = 134.9634° + 477198.8675° * T + 0,008721° * T²
Longitudine media del nodo ascendente dell'orbita lunare sull'ecclittica misurata a partire dall'equinozio medio della data
Ω = 125.0443° - 1934.1363° * T + 0,002075° * T²

Obliquità
Convertita Δε da arcosecondi a gradi decimali
Δε in ° = Δε / 3600
ε0 = 23°26'21".448 - 46".8150 * T - 0".00059 * T² + 0".001813 * T³
ε vera = ε0 + Δε
ε = ε vera + 0.00256 * cos Ω

Declinazione del sole
δ = arcsen (sen ε × sen La)

Ascensione Retta del sole
N = cos ε × sen LA
D = cos LA
tan α = N / D
Se N è positivo e D è positivo ==> α = α
Se N è positivo e D è negativo ==> α = α +180
Se N è negativo e D è positivo ==> α = α +360
Se N è negativo e D è negativo ==> α = α +180

Eccentricità dell'orbita
Ec = 0.016708634 - 0.000042037 × (T) - 0.0000001267 × (T)²

Equazione del tempo
ET hms = {[tan (ε ÷ 2)]² × sen (2 × L) - 2 × Ec × sen M + (4 × Ec) × [tan (ε ÷ 2)]² × sen M × cos (2 × L) - (1 ÷ 2) × [tan (ε ÷ 2)]ˆ4 × sen (4 × L) - (5 ÷ 4) × (Ec)² × sen (2 × M)} × 180° ÷ 3.14159265359 ÷ 15

Calcolo del sorgere, passagio al meridiano e tramonto del Sole
Vedi calcolo di sorgere, trasito e tramonto

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo delle coordinate lunari

Secoli Giuliani dall'epoca standard J2000 alla data inserita
T = (JD - 2451545.0) ÷ 36525

Longitudine media della luna incluso effetto tempo-luce
L1 = 218.3164477° + 481267.88123421° * T - 0.0015786° * T² + T³ / 538841 - T^4 / 65194000

Elongazione media della luna
D = 297.8501921 + 445267.1114034 * T - 0.0018819 * T² + T³ / 545868 - T^4 / 113065000

Anomalia media del sole
M = 357.5291092° + 35999.0502909° * T - 0.0001536° * T² + T³ / 24490000

Anomalia media della luna
M1 = 134.9633964° + 477198.8675055° * T + 0.0087414° * T² + T³ / 69699 - T^4 / 14712000

Distanza media della luna dal suo nodo ascendente
F = 93.2720950 + 483202.0175233 * T - 0.0036539 * T² - T³ / 3526000 + T^4 / 863310000

Solitamente questi angoli sono maggiori di 360° conviene quindi ridurli con la formula:
X = X - 360 * INT (X / 360)
dove INT è la parte intera della frazione (X / 360)

Coefficienti di calcolo
A1 = 119.75 + 131.849 * T
A2 = 53.09 + 479264.290 * T
A3 = 313.45 + 481266.484 * T
E = 1 - 0.002516 * T - 0.0000074 * T²

Calcolata la sommatoria di Σl e Σr della tabella Chapront ELP-2000/82

Moltiplicati i termini contenenti M oppure -M per il coefficiente E
Moltiplicati i termini contenenti 2M oppure -2M per il coefficiente E²

Esempio per i primi termini del coefficiente seno dell'argomento in tabella
Σl =
+ 6288774 * sen (M1)
+1274027 * sen (2D-M1)
+658314 * sen (2D)
+213618 * sen (2M1)
-185116 * E * (M)
-114332 * sen (2F)
omissis
omissis
omissis
+2236 * E² * sen (2D-2M)
ecc, ecc

Σr =
come sopra ma con il coefficiente coseno dell'argomento in tabella

Calcolata la sommatoria di Σβ della tabella Chapront ELP-2000/82

Moltiplicati i termini contenenti M oppure -M per il coefficiente E
Moltiplicati i termini contenenti 2M oppure -2M per il coefficiente E²

Σβ =
come per la tabella precedente

Aggiunto a Σl
+3958 * sen (A1)
+1962 * sen (L1 - F)
+318 * sen (A2)

Aggiunto a Σβ
-2235 * sen (L1)
+382 * sen (A3)
+175 * sen (A1 - F)
+175 * sen (A1 + F)
+127 * sen (L1 - M1)
-115 * sen (L1 + M1)

Coordinate eclittiche vere della luna
Λ = L1 + Σl / 1000000
β = Σβ / 1000000

Distanza della luna in km
Δ = 385000.56 + Σr / 1000

Parallasse lunare
π = arcsen (6378.14 / Δ)

Coordinate ecclittiche apparenti della luna
Λ app = Λ + Δψ (vedi caclolo tempo sidereo apparente)
β app = β

Coordinate equatoriali della luna (geocentriche)
sen δ = sen β * cos ε + cos β * sen ε * sen Λ

N = sen Λ * cos ε - tan β * sen ε
D = cos Λ

tan α = N / D

Se N è positivo e D è positivo ==> α = α
Se N è positivo e D è negativo ==> α = α +180
Se N è negativo e D è positivo ==> α = α +360
Se N è negativo e D è negativo ==> α = α +180

Convertite le coordinate equatoriali da geocentriche a topocentriche

a = 6378.14 (Raggio terrestre)
f = 1 / 298.257
b = a * (1 - f)
ba = 1 - f
u = tan(ba) * tan(latitudine)
h = metri di elevazione sul livello del mare
ro_sin_fi = ba * sin(u) + h / 6378140 * sin(latitudine)
ro_cos_fi = cos(u) + h / 6378140 * cos(latitudine)

ts = tempo siderale , ar = ascensione, H = angolo orario
H = ts - ar - longitudine
la longitudine va inserita negativa se di segno est e positiva se di segno ovest

Δα = tan(-ro_cos_fi * sin(parallasse) * sin(H) / cos(δ geocentrica) - ro_cos_fi * sin(parallasse) * cos(H))

α (topocentrica) = Δα + α (geocentrica)

δ topocentrica = tan(sin(δ geocentrica) - ro_sin_fi * sin(parallasse) * cos(Δα) / cos(δ geocentrica) - ro_cos_fi * sin(parallasse) * cos(H))

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo della frazione illuminata della luna

cos Φ = sin(δ sole) * sin(δ luna) + cos(δ sole) * cos(δ luna) * cos(α sole - α luna)
tan i = dist terra-sole * sin(Φ) / dist terra-luna - dist terra-sole * cos(Φ)
k = (1 + cos(i)) / 2
percentuale = parte intera arrotondata di (k * 100)

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Calcolo del centro (X) della parte illuminata della luna

Num = cos (dec sole) * sin (ar sole - ar luna)
Den = sin(dec sole) * cos(dec luna) - cos(dec sole) * sin(dec luna) * cos(ar sole - ar luna))
tan X = Num / Den
Se (X) > 0 la luna ha la gobba a levante (fase calante)
Se (X) < 0 la luna ha la gobba a ponente (fase crescente)

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo delle fasi lunari

Vengono calcolati i giorni occorsi dall'inizio dell'anno
INT = parte intera del calcolo
B = 1 per gli anno bisestili
B = 2 per gli anni NON bisestili
(in questo programma si è utilizzato B = 2)

Numero Giorni = INT(275 * mese / 9) - B * INT((mese + 9)/12) + giorno - 30

quindi l'anno decimale
anno.dec = Numero Giorni / 365 + anno

K = (anno_dec - 2000) * 12.3685

Una volta ottenuto k, quest'ultimo è preso ed inserito nelle successive formule come numero intero.
Es: se k = -35.34, viene preso solo il numero -35 oppure
se k = +16.07 viene preso solo il numero +16

Calcolo della fase luna nuova

Tempo in secoli giuliani dall'epoca 2000
T = K / 1236.85

Correzione per l'eccentricità
E = 1 - 0.002516 * T - 0.0000074 * T²

Giorno giuliano aprossimato
jd_1 = 2451550.09766 + 29.530588861 * K
+ 0.00015437 * T^2
- 0.000000150 * T^3
+ 0.00000000073 * T^4

Anomalia media del sole
M = 2.5534 + 29.10535670 * K
- 0.0000014 * T^2
- 0.00000011 * T^3

Anomalia media della luna
M1 = 201.5643 + 385.81693528 * K
+ 0.0107582 * T^2
+ 0.00001238 * T^3
- 0.000000058 * T^4

Argomento della latitudine della luna
F = 160.7108 + 390.67050284 * K
- 0.0016118 * T^2
- 0.00000227 * T^3
+ 0.000000011 * T^4

Longitudine del nodo ascendente dell'orbita lunare
O = 124.7746 - 1.56375588 * K
+ 0.0020672 * T^2
+ 0.00000215 * T^3

Prima correzione valida solo per luna nuova
C1 = - 0.40720 * sin(M1)
+ 0.17241 * E * sin(M)
+ 0.01608 * sin(2 * M1)
+ 0.01039 * sin(2 * F)
+ 0.00739 * E * sin(M1 - M)
- 0.00514 * E * sin(M1 + M)
+ 0.00208 * E^2 * sin(2 * M)
- 0.00111 * sin(M1 - 2 * F)
- 0.00057 * sin(M1 + 2 * F)
+ 0.00056 * E * sin(2 * M1 + M)
- 0.00042 * sin(3 * M1)
+ 0.00042 * E * sin(M + 2 * F)
+ 0.00038 * E * sin(M - 2 * F)
- 0.00024 * E * sin(2 * M1 - M)
- 0.00017 * sin(O)
- 0.00007 * sin(M1 + 2 * M)
+ 0.00004 * sin(2 * M1 - 2 * F)
+ 0.00004 * sin(3 * M)
+ 0.00003 * sin(M1 + M - 2 * F)
+ 0.00003 * sin(2 * M1 + 2 * F)
- 0.00003 * sin(M1 + M + 2 * F)
+ 0.00003 * sin(M1 - M + 2 * F)
- 0.00002 * sin(M1 - M - 2 * F)
- 0.00002 * sin(3 * M1 + M)
+ 0.00002 * sin(4 * M1)

Paranetri per la seconda correzione
a_1 = 299.77 + 0.107408 * K - 0.009173 * T^2
a_2 = 251.88 + 0.016321 * K
a_3 = 251.83 + 26.651886 * K
a_4 = 349.42 + 36.412478 * K
a_5 = 84.66 + 18.206239 * K
a_6 = 141.74 + 53.303771 * K
a_7 = 207.14 + 2.453732 * K
a_8 = 154.84 + 7.306860 * K
a_9 = 34.52 + 27.261239 * K
a_10 = 207.19 + 0.121824 * K
a_11 = 291.34 + 1.844379 * K
a_12 = 161.72 + 24.198154 * K
a_13 = 239.56 + 25.513099 * K
a_14 = 331.55 + 3.592518 * K

Seconda correzione
C2 = + 0.000325 * sin(a_1)
+ 0.000165 * sin(a_2)
+ 0.000164 * sin(a_3)
+ 0.000126 * sin(a_4)
+ 0.000110 * sin(a_5)
+ 0.000062 * sin(a_6)
+ 0.000060 * sin(a_7)
+ 0.000056 * sin(a_8)
+ 0.000047 * sin(a_9)
+ 0.000042 * sin(a_10)
+ 0.000040 * sin(a_11)
+ 0.000037 * sin(a_12)
+ 0.000035 * sin(a_13)
+ 0.000023 * sin(a_14)

Giorno Giuliano luna nuova
jd = jd_1 + C1 + C2

Il Giorno Giuliano viene poi convertito in data civile con le formule descritte nel calcolo degli equinozi e solstizi

Calcolo della fase luna piena

Al K calcolato per la luna nuova va aggiunto il valore 0.5
K = K + 0.5

Riprtuti i calcoli come sopra ma stavolta K, T e tutti gli altri parametri sono diversi

Tempo in secoli giuliani dall'epoca 2000
T = K / 1236.85

Correzione per l'eccentricità
E = 1 - 0.002516 * T - 0.0000074 * T²

Giorno giuliano aprossimato
jd_1 = 2451550.09766 + 29.530588861 * K
+ 0.00015437 * T^2
- 0.000000150 * T^3
+ 0.00000000073 * T^4

Anomalia media del sole
M = 2.5534 + 29.10535670 * K
- 0.0000014 * T^2
- 0.00000011 * T^3

Anomalia media della luna
M1 = 201.5643 + 385.81693528 * K
+ 0.0107582 * T^2
+ 0.00001238 * T^3
- 0.000000058 * T^4

Argomento della latitudine della luna
F = 160.7108 + 390.67050284 * K
- 0.0016118 * T^2
- 0.00000227 * T^3
+ 0.000000011 * T^4

Longitudine del nodo ascendente dell'orbita lunare
O = 124.7746 - 1.56375588 * K
+ 0.0020672 * T^2
+ 0.00000215 * T^3

Prima correzione solo per luna piena
C1 = - 0.40614 * sin(M1)
+ 0.17302 * E * sin(M)
+ 0.01614 * sin(2 * M1)
+ 0.01043 * sin(2 * F)
+ 0.00734 * E * sin(M1 - M)
- 0.00515 * E * sin(M1 + M)
+ 0.00209 * E^2 * sin(2 * M)
- 0.00111 * sin(M1 - 2 * F)
- 0.00057 * sin(M1 + 2 * F)
+ 0.00056 * E * sin(2 * M1 + M)
- 0.00042 * sin(3 * M1)
+ 0.00042 * E * sin(M + 2 * F)
+ 0.00038 * E * sin(M - 2 * F)
- 0.00024 * E * sin(2 * M1 - M)
- 0.00017 * sin(O)
- 0.00007 * sin(M1 + 2 * M)
+ 0.00004 * sin(2 * M1 - 2 * F)
+ 0.00004 * sin(3 * M)
+ 0.00003 * sin(M1 + M - 2 * F)
+ 0.00003 * sin(2 * M1 + 2 * F)
- 0.00003 * sin(M1 + M + 2 * F)
+ 0.00003 * sin(M1 - M + 2 * F)
- 0.00002 * sin(M1 - M - 2 * F)
- 0.00002 * sin(3 * M1 + M)
+ 0.00002 * sin(4 * M1)

Giorno giuliano luna piena
jd = jd_1 + C1 + C2

Il Giorno Giuliano viene poi convertito in data civile con le formule descritte nel calcolo degli equinozi e solstizi

Calcolo della fase primo quarto ed ultimo quarto

Al K calcolato per la luna nuova va aggiunto il valore 0.25 per il primo quarto e 0.75 per l'ultimo quarto
K = K + 0.25 (primo quarto)
K = K + 0.75 (ultimo quarto)

Ripetuti i calcoli come sopra ma stavolta K, T e tutti gli altri parametri sono diversi

Tempo in secoli giuliani dall'epoca 2000
T = K / 1236.85

Correzione per l'eccentricità
E = 1 - 0.002516 * T - 0.0000074 * T²

Giorno giuliano aprossimato
jd_1 = 2451550.09766 + 29.530588861 * K
+ 0.00015437 * T^2
- 0.000000150 * T^3
+ 0.00000000073 * T^4

Anomalia media del sole
M = 2.5534 + 29.10535670 * K
- 0.0000014 * T^2
- 0.00000011 * T^3

Anomalia media della luna
M1 = 201.5643 + 385.81693528 * K
+ 0.0107582 * T^2
+ 0.00001238 * T^3
- 0.000000058 * T^4

Argomento della latitudine della luna
F = 160.7108 + 390.67050284 * K
- 0.0016118 * T^2
- 0.00000227 * T^3
+ 0.000000011 * T^4

Longitudine del nodo ascendente dell'orbita lunare
O = 124.7746 - 1.56375588 * K
+ 0.0020672 * T^2
+ 0.00000215 * T^3

Prima correzione valida per luna al primo ed all'ultimo quarto
C1 = - 0.62801 * sin(M1)
+ 0.17172 * E * sin(M)
- 0.01183 * E * sin(M1 + M)
+ 0.00862 * sin(2 * M1)
+ 0.00804 * sin(2 * F)
+ 0.00454 * E * sin(M1 - M)
+ 0.00204 * E^2 * sin(2 * M)
- 0.00180 * sin(M1 - 2 * F)
- 0.00070 * sin(M1 + 2 * F)
- 0.00040 * sin(3 * M1)
- 0.00034 * E * sin(2 * M1 - M)
+ 0.00032 * E * sin(M + 2 * F)
+ 0.00032 * E * sin(M - 2 * F)
- 0.00028 * E^2 * sin(M1 + 2 * M)
+ 0.00027 * E * sin(2 * M1 + M)
- 0.00017 * sin(O)
- 0.00005 * sin(M1 - M - 2 * F)
+ 0.00004 * sin(2 * M1 + 2 * F)
- 0.00004 * sin(M1 + M + 2 * F)
+ 0.00004 * sin(M1 - 2 * M)
+ 0.00003 * sin(M1 + M - 2 * F)
+ 0.00003 * sin(3 * M)
+ 0.00002 * sin(2 * M1 - 2 * F)
+ 0.00002 * sin(M1 - M + 2 * F)
- 0.00002 * sin(3 * M1 + M)

C1_BIS = 0.00306 - 0.00038 * E * cos(M) + 0.00026 * cos(M1) - 0.00002 * cos(M1 - M) + 0.00002 * cos(M1 - M) + 0.00002 * cos(2 * F)

Giorno giuliano luna al primo quarto
jd = jd_1 + C1 + C1_BIS + C2

Giorno giuliano luna all'ultimo quarto
jd = jd_1 + C1 - C1_BIS + C2

Il Giorno Giuliano viene poi convertito in data civile con le formule descritte nel calcolo degli equinozi e solstizi

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo FK5 J2000 della posizione apparente di una stella

Nr. di secoli a partire dal 1/1/2000
T = (JD - 2451545.0) / 36525

Correzione dei moti propri nello spazio

α0 e δ0 coordinate equatoriali all'epoca di riferimento J2000
r = distanza della stella in parsec
Δr = velocita' radiale in parsec / anno
Δα e Δδ = moti propri in arcosecondi / annno

x = r cos δ0 cos α0
y = r cos δ0 sin α0
z = r sin δ0

Δx = (x/r) Δr - zΔδ cos α0 - yΔα
Δy = (y/r)Δr - zΔδ sin α0 + xΔα
Δz = (z/r)Δr + rΔδ cosδ0

x1 = x + tΔx
y1 = y + tΔy
z1 = z + tΔz

tan α = y1 / x1 (sin α con lo stesso segno di Y1)
tan δ = z1 / Radice quadrata di x1² + y1²

Correzione di α0 e δ0 per la precessione

ζ = 2306.2181" * T + 0.30188" * T² - 0.017998" * T³
z = 2306.2181" * T + 1.09468" * T² + 0.018203" * T³
θ = 2004.3109" * T - 0.42665" * T² - 0.041833" * T³
A = cos δ0 * sin (α0 + ζ)
B = cos θ * cos δ0 * cos (α0 + ζ) - sin θ * sin δ0
C = sin θ * cos δ0 * cos (α0 + ζ) + cos θ * sin δ0
tan (α1 - z) = (A / B)
sin δ1 = C
α1 = tan(A/B) + z
δ1 = asin(C)

Calcolo della nutazione e obliquità dell'ecclittica

Longitudine media del sole
L = 280.4665° + 36000.7698° + 0.000303° * T²
Longitudine media della luna
L1 = 218.3164° + 481267.8812° * T - 0.001599° * T²
Anomalia media del sole
M = 357.5291° + 35999.0503° * T - 0.000154° * T²
Anomalia media della luna
M1 = 134.9634° + 477198.8675° * T + 0.008721° * T²
Longitudine media del nodo ascendente dell'orbita lunare sull'ecclittica misurata a partire dall'equinozio medio della data
Ω = 125.0443° - 1934.1363° * T + 0.002075° * T²

Obliquità
Convertire Δε da arcosecondi a gradi decimali
Δε in ° = Δε / 3600
ε0 = 23°26'21".448 - 46".8150 * T - 0".00059 * T² + 0".001813 * T³
ε = ε0 + Δε

Correzione Δ da apportare ad α1 e δ1 per la nutazione e obliquità dell'ecclittica

Convertita Δψ da arcosecondi a gradi decimali
Δψ in ° = Δψ / 3600
Δα2 = (cos ε + sen ε * sen α1 * tan δ1) * Δψ - (cos α1 * tan δ1) * Δε
Δδ2 = (sen ε * cos α1) * Δψ + (sen α1) * Δε

Calcolo dell’aberrazione annua della luce

Longitudine media del sole (stessa formula per la nutazione vedi sopra)
L = 280.4665° + 36000.7698° + 0.000303° * T²

Anomalia media del sole (stessa formula per la nutazione vedi sopra)
M = 357.5291° + 35999.0503° * T - 0,000154° * T²

Equazione del centro del sole
C = + (1.914602° - 0.004817° * T - 0.000014° * T²) * sin M + (0.019993° - 0.000101° * T) * sin 2M + 0.000289° * sin 3M

Longitudine vera del sole
Lv = L + C

Eccentricità dell'orbita
e = 0.016708634 - 0.000042037 * T - 0.0000001267 * T²

Longitudine del perielio
Π = 102.93735 + 1.71946 * T + 0.00046 * T²

Costante dell'aberrazione
k = 20.49552"

Correzione Δ da apportare ad α1 e δ1 per l' aberrazione

Convertita k da arcosecondi a gradi decimali
Δα3= - k * [(cos α1 * cos Lv * cos ε + sen α1 * sen Lv) / cos δ1]
+ e * k * [(cos α1 * cos Π * cos ε + sen α1 * sen Π) / cos δ1]

Δδ3 = - k * [cos Lv * cos ε * (tan ε * cos δ1 - sen α1 * sen δ1) + cos α1 * sen δ1*sen L1]
+ e * k [cos Π * cos ε * (tan ε* cos δ1 - sin α1 * sin δ1) + cos α1 * sin δ1 * sin Π]

Coordinate della posizione apparente dell' astro
α app = α1 + Δα2 + Δα3
δ app = δ1 + Δδ2 + Δδ3

(Calcules Astronomiques - Jean Meeus 2014)

Torna a inizio pagina

Calcolo delle coordinate planetarie

Il calcolo della posizione dei pianeti è stato effettuato con l'estratto dei più importanti termini periodici della teoria Bretagnon Francou VSOP-87, questa fornisce 3 parametri orbitali dei pianeti:
L (Longitudine eclittica)
B (Latitudine eclittica)
R (distanza in Unità astronomiche)

Calcolati i millenni dall'epoca J2000 all'istante dato con la seguente formula:
ro = (JD - 2451545.0) / 365250

una volta ottenuto ro, il valore di ogni singolo termine è stato calcolato con:
A * cos(B + C*ro)
dove A,B e C sono i termini periodici descritti nelle tabelle dei pianeti (pulsanti "T. Periodici")
le quantità B e C sono espresse in radianti, il coefficiente A è espresso in unità di 10^-8 radianti

Mercurio , Venere , Terra , Marte ,

Giove , Saturno , Urano , Nettuno

La sommatoria dei termini così ottenuti fornisce i riultati parziali di L, B ed R
L0 = Σ [A * cos(B + C*ro)]
L1 = Σ [A * cos(B + C*ro)]
L2 = Σ [A * cos(B + C*ro)]
L3 = Σ [A * cos(B + C*ro)]
L4 = Σ [A * cos(B + C*ro)]
--- ecc ecc ---
B0 = Σ [A * cos(B + C*ro)]
B1 = Σ [A * cos(B + C*ro)]
B2 = Σ [A * cos(B + C*ro)]
B3 = Σ [A * cos(B + C*ro)]
B4 = Σ [A * cos(B + C*ro)]
--- ecc ecc ---
R0 = Σ [A * cos(B + C*ro)]
R1 = Σ [A * cos(B + C*ro)]
R2 = Σ [A * cos(B + C*ro)]
R3 = Σ [A * cos(B + C*ro)]
R4 = Σ [A * cos(B + C*ro)]
--- ecc ecc ----

Dopo aver calcolato tutti i termini periodici si sono ottenuti i risultati eliocentrici finali di L,B ed R
L = (L0 + L1*ro + L2*ro^2 + L3*ro^3 + L4*ro^4 + L5*ro^5)/10^8
B = (B0 + B1*ro + B2*ro^2 + B3*ro^3 + B4*ro^4 + B5*ro^5)10^8
R = (R0 + R1*ro + R2*ro^2 + R3*ro^3 + R4*ro^4 + R5*ro^5)10^8

I risultati sono riferiti all'equinozio della data definita dalla teoria planetaria VSOP
che differisce di poco rispetto allo standard FK5
per convertire questi dati allo standard FK5 sono state utilizzate le seguenti formule:
T = (JD - 2451545.0) / 36525
L1 = L - 1°.397 * T - 0°.00031 * T²
ΔL = -0".09033 + 0".03916 * (cos L1 + sin L1) * tan B
ΔB = +0".03916 * (cos L1 - sin L1)
L (Sistema FK5) = L + ΔL
B (Sistema FK5) = B + ΔB

Le coordinate eclittiche geocentriche sono state ottenute mettendo nel sistema anche il pianeta terra
ovvero calcolati L,B ed R del pianeta interessato, poi LE, BE ed RE della terra NON convertite nel sistema FK5
quindi sono state calcolate le coordinate rettangolari X,Y,Z del pianeta
X = R * cos B * cos L - RE * cos BE * cos LE
Y = R * cos B * sin L - RE * cos BE * sin LE
Z = R * sin B - RE * sin BE

Quindi le coordinate geocentriche sono state ricavate da:
tan L(Geo) = Y / X
tan B(Geo) = Z / √ (X² + Y²)

Queste sono però le coordinate geocentriche del pianeta all'istante considerato
nella realtà invece la luce del pianeta arriva a noi diversi minuti dopo, per ottenere
le coordinate geocentriche reali è stata tenuta in considerazione la distanza che la luce
impiega per arrivare fino a noi.
Calcolata la distanza dal pianeta considerato alla terra con la seguente formula:
Δ = √ (X² + Y² + Z²)
ro1 = 0.0057755183 * Δ

Con la nuova ro1 sono stati ripetuti i calcoli come descritto sopra e sono state ottenute nuovamente
L,B ed R che sono state convertite nel sistena FK5 come descritto precedentemente.

Al solo valore di L ottenuto è stata aggiunta la correzione per la nutazione Δψ
descritta nel "Calcolo del tempo sidereo apparente"

Con le coordinate eclittiche finali sono state calcolate le coordinate equatoriali come descritto nel calcolo delle coordinate solari.

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo della coordinate Altazimutali

In questa formula λ Est ha segno negativo mentre λ Ovest ha segno positivo
H è l'angolo orario locale, misurato in direzione Ovest, da Sud
H = Ts - α - λ
P = Angolo al polo
Se H > 180 ==> P = 360 - H (ed è di segno positivo est)
Se H < 180 ==> P = H (ed è di segno negativo ovest)

sen h = sen φ * sen δ + cos φ * cos δ * cos P

N = sen P
D = tan δ * cos φ - sen φ * cos P

tan az = N / D

Se N è positivo e D è positivo ==> az = az
Se N è positivo e D è negativo ==> az = az +180
Se N è negativo e D è positivo ==> az = az +360
Se N è negativo e D è negativo ==> az = az +180

Torna a inizio pagina

Calcolo della coordinate Ecclittiche

sen β = sen δ * cos ε - cos δ * sen ε * sen α

N = sen α * cos ε + tan δ * sen ε
D = cos α

tan λ = N / D

Se N è positivo e D è positivo ==> λ = λ
Se N è positivo e D è negativo ==> λ = λ +180
Se N è negativo e D è positivo ==> λ = λ +360
Se N è negativo e D è negativo ==> λ = λ +180

Torna a inizio pagina

Calcolo del sorgere, passaggio al meridiano e tramonto

φ = latitudine , λ = longitudine, L = λ * -1

hο = -0°34' = -0.5667°dec (per stelle e pianeti)
hο = -0°50' = -0.8333°dec (per il sole)
hο = -0.7275°dec * π - 0.5667 π = parallasse lunare (per la luna)
π = parallasse lunare
se non si ha la possibilità di calcolare la parallasse lunare un valore approssimato di ho per la luna è +0.125

Questo programma al posto della rifrazione media approssimata pari a 0°34' calcola la rifrazione esatta con le formule descritte successivamente.

Sono state calcolate le ascensioni rette e declinazioni relative all'astro del giorno prima, del giorno in corso e del giorno dopo
alle ore 0 di tempo dinamico
(α1 α2 α3 δ1 δ2 δ3)

E' stato calcolato il tempo sidereo del giorno in corso alle ore 0 di UT (Φo)

Cos H = (sen hο - sen φ sen δ2) / (cos φ cos δ2)
se il secondo membro di questa formula è maggiore di 1 oppure minore di -1 allora l'astro è circumpolare

E' stato calolato mo (transito), m1 (sorgere), m2 (tramonto) con le seguenti formule:
mo = (α2 + L - Φo) / 360
m1 = mo - H / 360
m2 = mo + H / 360

i valori di m devono essere compresi tra 0 e + 1, ai valori fuori da questo range va aggiunto o sottratto il valore 1
Es: -0.1709 diventa +0.8291 mentre +1.1853 diventa 0.1853

Sono stati calcolati i tempi siderei (θo, θ1, θ2) per i valori di mo, m1, m2 con le seguenti formule:
θo = Φo + 360.985647 * mo
θ1 = Φo + 360.985647 * m1
θ2 = Φo + 360.985647 * m2

Sono stati interpolati i valori di Ascensione retta e declinazione del giorno prima, del giorno in corso e del giorno dopo
utilizzando i seguenti fattori di interpolazione:
n0 = mo + ΔT/86400
n1 = m1 + ΔT/86400
n2 = m2 + ΔT/86400
ΔT = Delta T (differenza tra Tempo Dinamico e Tempo Universale)

Formule di interpolazione generiche:
conoscendo y1,y2,y3
a = y2 - y1
b = y3 - y2
c = b - a
y = y2 + n/2 * (a + b + n*c)

Sostituendo i valori di y1, y2, y3 prima con le 3 ascensioni rette e poi con le 3 declinazioni
è stata trovata una singola ascensione retta ed una singola declinazione interpolata per il fattore di interpolazione n0
(α0 e δ0 relativa al fattore di interpolazione n0)

ripetuta l'interpolazione per il fattore n1
(α1 e δ1 relativa al fattore di interpolazione n0)

ripetuta l'interpolazione per il fattore n2
(α2 e δ2 relativa al fattore di interpolazione n2)

Con i valori di (θo θ1 θ2) sono stati calcolati i relativi angoli orari con le seguenti formule:
Ho = θo - L - α0 (meridiano)
H1 = θ1 - L - α1 (sorgere)
H2 = θ2 - L - α2 (tramonto)

Sono state calcolate le altezze degli astri per gli istanti del sorgere e del tramonto con la seguente formula:
sen h1 = sen φ * sen δ1 + cos φ * cos δ1 * cos H1 (sorgere)
sen h2 = sen φ * sen δ2 + cos φ * cos δ2 * cos H2 (tramonto)

E' stata calcolate la differenza di m da aggiungere ai valori di m trovati in precedenza
Δm0 = - H0 / 360 (per il transito)
Δm1 = (h1 - ho) / 360 * cos δ1 * cos φ * cos H1 (per il sorgere)
Δm2 = (h2 - ho) / 360 * cos δ2 * cos φ * cos H2 (per il tramonto)

Aggiunti i (Δm) ai valori di (m) calcolati in precedenza e convertiti in 24 esimi
per ottenere gli orari UT

Passaggio in meridiano = (Δm0 + mo) * 24
Sorgere = (Δm1 + m1) * 24
Tramonto = (Δm2 + m2) * 24

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo della rifrazione

Calcolo della rifrazione astronomica con pressione atmosferica a 1010 millibar ed una temperatura di 10° Celsius
(formula di Bennett)

1
R = --------------------------------------
7.31
tan ( ho + ------------ )
ho + 4.4

L' altezza osservata per il calcolo del sorgere e tramonto è pari a 0

Affinamento della formula

A = -0.06 * sin(14.7 * (rifrazione/60) + 13)

Correzione della rifrazione astronomica in funzione della pressione atmosferica e della temperatura esterna

correzione = (P / 1010) x (283 / (273 + T ))

dove P è la pressione atmosferica espressa in millibar e T è la temperatura esterna espressa in gradi celsius

Rifrazione corretta = (R+A) * correzione

Il risultato è espresso in primi d'arco ed è stato poi convertito in gradi dividendolo per 60

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo degli elementi orbitali dei pianeti

Descrizione degli elementi:
L = Longitudine media del pianeta
a = Semiassemaggiore dell'orbita
e = Eccentricità dell'orbita
i = inclinazione sul piano dell'eclittica
Ω = Longitudine del nodo ascendente
π = Longitudine del perielio

Gli elementi orbitali dei pianeti sono stati calcolati tramite polinomi nella seguente forma:

a0 + a1*T + a2*T² + a3*T³

dove T = (JD - 2451545.0) / 36525

i valori di a0,a1,a2,a3 sono descritti nella tabella sottostante

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo perielio ed afelio dei pianeti

Calcolo del perielio

Vengono calcolati i giorni occorsi dall'inizio dell'anno
INT = parte intera del calcolo
B = 1 per gli anno bisestili
B = 2 per gli anni NON bisestili
(in questo programma si è utilizzato B = 2)

Numero Giorni = INT(275 * mese / 9) - B * INT((mese + 9)/12) + giorno - 30

quindi l'anno decimale
anno.dec = Numero Giorni / 365 + anno

Vengono calcolate le costanti k per ogni pianeta
(Mercurio) k = 4.15201 * (anno.dec - 2000.12)
(Venere) k = 1.62549 * (anno.dec - 2000.53)
(Terra) k = 0.99997 * (anno.dec - 2000.01)
(Marte) k = 0.53166 * (anno.dec - 2001.78)
(Giove) k = 0.08430 * (anno.dec - 2011.20)
(Saturno) k = 0.03393 * (anno.dec - 2003.52)
(Urano) k = 0.01190 * (anno.dec - 2051.1)
(Nettuno) k = 0.00607 * (anno.dec - 2047.5)

i relativi Giorni Giuliani del perielio derivano dalle seguenti formule
(Mercurio) JD = 2451590.257 + 87.96934963 * k - 0.0000000000 * k²
(Venere) JD = 2451738.233 + 224.7008188 * k - 0.0000000327 * k²
(Terra) JD = 2451547.507 + 365.2596358 * k + 0.0000000156 * k²
(Marte) JD = 2452195.026 + 686.9957857 * k - 0.0000001187 * k²
(Giove) JD = 2455636.936 + 4332.897065 * k + 0.0001367 * k²
(Saturno) JD = 2452830.12 + 10764.21676 * k + 0.000827 * k²
(Urano) JD = 2470213.5 + 30694.8767 * k - 0.00541 * k²
(Nettuno) JD = 2468895.1 + 60190.33 * k + 0.03429 * k²

Ottenuti i Giorni Giuliani vengono calcolate le date gregoriane con le formule espresse
negli algoritmi "Calcolo degli Equinozi e Solstizi

Nei risultati di questo programma gli orari sono espressi in Tempo dinamico
per convertirli in Tempo universale basta sottrarre il ΔT.
UT = ET - ΔT

Calcolo dell' afelio

Stessa procedura come sopra ponendo
k = k + 0.5

Per il solo pianeta Terra i giorni giuliani del perielio ed afelio
sono stati corretti per ottenere un risultato più preciso con le seguenti formule:

Perielio
A1 = 328.41 + 132.788585 * k
A2 = 316.13 + 584.903153 * k
A3 = 346.20 + 450.380738 * k
A4 = 136.95 + 659.306737 * k
A5 = 249.52 + 329.653368 * k

correzione = +1.278 * sin(A1) - 0.055 * sin(A2) - 0.091 * sin(A3) - 0.056 * sin(A4) - 0.045 * sin(A5)

JD corretto = JD perielio + correzione

Afelio
ponendo k = k + 0.5
A1 = 328.41 + 132.788585 * k
A2 = 316.13 + 584.903153 * k
A3 = 346.20 + 450.380738 * k
A4 = 136.95 + 659.306737 * k
A5 = 249.52 + 329.653368 * k

correzione = -1.352 * sin(A1) + 0.061 * sin(A2) + 0.062 * sin(A3) + 0.029 * sin(A4) + 0.031 * sin(A5)

JD corretto = JD afelio + correzione

Il Giorno Giuliano viene poi convertito in data civile con le formule descritte nel calcolo degli equinozi e solstizi

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

Torna a inizio pagina

Calcolo di congiunzioni ed elongazioni

Per il calcolo delle congiunzioni Superiori ed inferiori dei pianeti
e per il calcolo delle elongazioni massime di Mercurio e Venere
si è fatto riferimento al Capitolo 36 del libro:
Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009
Qui troverete i termini periodici nelle tabelle 36A, 36B, 36C necessarie al calcolo.

Torna a inizio pagina

Calcolo degli equinozi e solstizi

y = Anno / 1000

Y = (anno - 2000) / 1000

Per l'equinozio di Marzo
Giorno Giuliano iniziale vicino all' equinozio di Marzo per l'anno considerato
se l'anno è inferiore a 1000
JD = 1721139.29189 + 365242.13740 * y + 0.06134 * y² + 0.00111 * y³ - 0.00071 * y^4
se l'anno è superiore a 1000
JD = 2451623.80984 + 365242.37404 * Y + 0.05169 * Y² - 0.00411 * Y³ - 0.00057 * Y^4

Per il solstizio di Giugno
Giorno Giuliano iniziale vicino al solstizio di Giugno per l'anno considerato
se l'anno è inferiore a 1000
JD = 1721233.25401 + 365241.72562 * y - 0.05323 * y² + 0.00907 * y³ + 0.00025 * y^4
se l'anno è superiore a 1000
JD = 2451716.56767 + 365241.62603 * Y + 0.00325 * Y² + 0.00888 * Y³ - 0.00030 * Y^4

Per l'equinozio di Settembre
Giorno Giuliano iniziale vicino all' equinozio di Settembre per l'anno considerato
se l'anno è inferiore a 1000
JD = 1721325.70455 + 365242.49558 * y - 0.11677 * y² - 0.00297 * y³ + 0.00074 * y^4
se l'anno è superiore a 1000
JD = 2451810.21715 + 365242.01767 * Y - 0.11575 * Y² + 0.00337 * Y³ + 0.00078 * Y^4

Per il solstizio di Dicembre
Giorno Giuliano iniziale vicino al solstizio di Dicembre per l'anno considerato
se l'anno è inferiore a 1000
JD = 1721414.39987 + 365242.88257 * y - 0.00769 * y² - 0.00933 * y³ - 0.00006 * y^4
se l'anno è superiore a 1000
JD = 2451900.05952 + 365242.74049 * Y - 0.06223 * Y² - 0.00823 * Y³ + 0.00032 * Y^4

Una volta trovato il JD interessato, tramite il calcolo della posizione dei pianeti relativo alla terra con l'estratto dei più importanti termini periodici della teoria Bretagnon Francou VSOP-87
(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009 - Cap 33 + Appendix III)
si sono ricavate:
Longitudine dell'eclittica (L_VSOP87) espresso in °decimali
Raggio vettore (R_VSOP87) espresso in UA

E' stata calcolata la nutazione in longitudine con le seguenti formule:
(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009 - Cap 22)

Nr. di secoli a partire dal 1/1/2000
T = (JD - 2451545.0) / 36525

Longitudine media del sole
L = 280.4665° + 36000.7698° + 0.000303° * T²
Longitudine media della luna
L1 = 218.3164° + 481267.8812° * T - 0.001599° * T²
Anomalia media del sole
M = 357.5291° + 35999.0503° * T - 0.000154° * T²
Anomalia media della luna
M1 = 134.9634° + 477198.8675° * T + 0.008721° * T²
Longitudine media del nodo ascendente dell'orbita lunare sull'ecclittica misurata a partire dall'equinozio medio della data
Ω = 125.0443° - 1934.1363° * T + 0.002075° * T²

Convertita Δψ da arcosecondi a gradi decimali Δψ° = Δψ / 3600

dopo aver trovato la nutazione e sapendo che la correzione per il sistema FK5 = - 0.09033" arcosecondi
e stata calcolata CORR_FK5° = - 0.09033 / 3600

ricavata l'aberrazione tramite la formula:
(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009 - Cap 25)
AB = - 20.4898" / R_VSOP87
quindi AB° = AB / 3600

Calcolata la longitudine apparente del sole ovvero quella opposta al pianeta terra tramite la formula:
LA_sole = L_VOSOP87 - 180 + Δψ° + CORR_FK5° + AB°

con questa longitudine appena ricavata si è calcolata la correzione (c) che fornisce la quantità in giorni giuliani da aggiungere al JD iniziale per avvicinarsi all'orario dell' equinozio o solstizio preso in considerazione

Per l' equinozio di Marzo
Correzione in giorni giuliani da apportare al JD
c = 58 * sen (0 * 90 - LA_sole)

Per il solstizio di Giugno
Correzione in giorni giuliani da apportare al JD
c = 58 * sen (1 * 90 - LA_sole)

Per l' equinozio di Settembre
Correzione in giorni giuliani da apportare al JD
c = 58 * sen (2 * 90 - LA_sole)

Per il solstizio di Dicembre
Correzione in giorni giuliani da apportare al JD
c = 58 * sen (3 * 90 - LA_sole)

Nuovo Giorno Giuliano
JDN = JD + c

Con il nuovo Giorno Giuliano si sono ripetuti i calcoli dall'inizio ponendo JDN al posto di JD fino ad avere una longitudine apparente del sole prossima a 0° per Marzo, 90° per Giugno, 180 per Settembre, 270 per Dicembre; questa condizione si ha quando il valore di c è prossino a 0.000005 (ovviamente questo è possibile solo tramite un ciclo di programmazione per computer)

con l'ultimo JDN ottenuto, ed avendo c prossimo a 0.000005 è stato convertito il giorno giuliano in data Gregoriana con le seguenti formule:

Aggiunto 0.5 al Giorno Giuliano, siano Z la sua Parte Intera e F la sua Parte frazionaria
Se Z < 2299161 si è preso A = Z
Se Z è maggiore o uguale a 2299161 si è calcolato
α = Parte Intera ((Z - 1867216.25) / 36524.25)
A = Z + 1 + α - Parte Intera (α / 4)
B = A + 1524
C = Parte Intera ((B - 122.1) / 365.25)
D = Parte Intera (365.25 * C)
E = Parte Intera ((B - D) / 30.6001)

Il giorno del mese (con decimali) = B - D - Parte Intera (30.6001 * E) + F

Il mese = E - 1 se E < 13.5
Il mese = E - 13 se E > 13.5

L'anno = C - 4716 se m > 2.5
L'anno = C - 4715 se m < 2.5

Gli orari ottenuti sono in Tempo Effemeridi per convertirli in Tempo universale basta sottrarre il ΔT.

UT = ET - ΔT

(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009)

In questo programma il ΔT è stato calcolato con le formule fornite dal sito NASA
POLYNOMIAL EXPRESSIONS FOR DELTA T (ΔT)
Five Millennium Canon of Solar Eclipses [Espenak and Meeus]

https://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html