Calcolo degli equinozi e solstizi
(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009 - Cap 27)

y = Anno / 1000

Y = (anno - 2000) / 1000

Nelle formule che seguono fate attenzione a non confondere y (minuscola) con Y (maiuscola)
in pratica per gli anni inferiori a 1000 si usa y (minuscola)
mentre per gli anni superiori a 1000 si usa Y (maiuscola)

Per l'equinozio di Marzo
Giorno Giuliano iniziale vicino all' equinozio di Marzo per l'anno considerato
se l'anno è inferiore a 1000
JDE = 1721139.29189 + 365242.13740 * y + 0.06134 * y² + 0.00111 * y³ - 0.00071 * y^4
se l'anno è superiore a 1000
JDE = 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
JDE = 1721233.25401 + 365241.72562 * y - 0.05323 * y² + 0.00907 * y³ + 0.00025 * y^4
se l'anno è superiore a 1000
JDE = 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
JDE = 1721325.70455 + 365242.49558 * y - 0.11677 * y² - 0.00297 * y³ + 0.00074 * y^4
se l'anno è superiore a 1000
JDE = 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
JDE = 1721414.39987 + 365242.88257 * y - 0.00769 * y² - 0.00933 * y³ - 0.00006 * y^4
se l'anno è superiore a 1000
JDE = 2451900.05952 + 365242.74049 * Y - 0.06223 * Y² - 0.00823 * Y³ + 0.00032 * Y^4

Una volta trovato il JDE 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 ricava:
Longitudine dell'eclittica (L_VSOP87) espresso in °decimali
Raggio vettore (R_VSOP87) espresso in UA

Si calcola ora 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_UT - 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)

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

dopo aver trovato la nutazione e sapendo che la correzione per il sistema FK5 = - 0.09033" arcosecondi
(Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009 - Cap 25)
quindi CORR_FK5° = - 0.09033 / 3600

ricaviamo 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

ora non ci resta che trovare 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 calcola la correzione (c) che fornisce la quantità in giorni giuliani da aggiungere al JDE iniziale per avvicinarsi all'orario dell' equinozio o solstizio preso in considerazione

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

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

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

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

Nuovo Giorno Giuliano
JDN = JDE + c

Con il nuovo Giorno Giuliano si ripetono i calcoli dall'inizio ponendo JDN al posto di JDE 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 si converte il giorno giuliano in data Gregoriana con le seguenti formule: (Astronomical Algorithms Second Edition - Jean Meeus 1998 corr. 2009 - Cap 7)

Si aggiunga 0.5 al Giorno Giuliano e siano Z la sua Parte Intera e F la sua Parte frazionaria
Se Z < 2299161 si prenda A = Z
Se Z è maggiore o uguale a 2299161 si calcoli
α = 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

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

ET = UT + ΔT

UT = ET - ΔT

Polynomial Expressions for Delta T (ΔT)

Adapted from "Five Millennium Canon of Solar Eclipses" [Espenak and Meeus]

Using the ΔT values derived from the historical record and from direct observations.

We define the decimal year "y" as follows:

y = year + (month - 0.5)/12

Before the year -500, calculate:

Between years -500 and +500 calculate:

Between years +500 and +1600 calculate:

Between years +1600 and +1700, calculate:

Between years +1700 and +1800, calculate:

Between years +1800 and +1860, calculate:

Between years 1860 and 1900, calculate:

Between years 1900 and 1920, calculate:

Between years 1920 and 1941, calculate:

Between years 1941 and 1961, calculate:

Between years 1961 and 1986, calculate:

Between years 1986 and 2005, calculate:

Between years 2005 and 2050, calculate:

This expression is derived from estimated values of ΔT in the years 2010 and 2050. The value for 2010 (66.9 seconds) is based on a linearly extrapolation from 2005 using 0.39 seconds/year (average from 1995 to 2005). The value for 2050 (93 seconds) is linearly extrapolated from 2010 using 0.66 seconds/year (average rate from 1901 to 2000).

Between years 2050 and 2150, calculate:

The last term is introduced to eliminate the discontinuity at 2050.

After 2150, calculate:

All values of ΔT based on Morrison and Stephenson [2004] assume a value for the Moon's secular acceleration of -26 arcsec/cy^2. However, the ELP-2000/82 lunar ephemeris employed in the Canon uses a slightly different value of -25.858 arcsec/cy^2. Thus, a small correction "c" must be added to the values derived from the polynomial expressions for ΔT before they can be used in the Canon

Since the values of ΔT for the interval 1955 to 2005 were derived independent of any lunar ephemeris, no correction is needed for this period.

The uncertainty in ΔT over this period can be estimated from scatter in the measurements.