ARCHEOASTRONOMIA LIGUSTICA
Pubblicato
in: Atti del XIII seminario A.L.S.S.A. di Archeoastronomia, Genova 09-10 aprile
2011. pp. 43-57.
Il calcolo FK4 B1950.0 della
precessione delle stelle
Mario Codebò
Si descrive qui l’algoritmo per il
calcolo della posizione apparente di una stella ad una determinata data
- in Archeoastronomia quasi sempre nel passato più o meno remoto - partendo
dalle sue coordinate per un’epoca standard (per es. il 1900, il
1950, il 2000 o, prossimamente, il 2050)[1].
Tali coordinate sono:
1) l’ascensione retta α, comunemente
espressa in unità sessagesimali del tempo (ore h, minuti m e
secondi s);
2) la declinazione δ,
espressa in unità di gradi sessagesimali di circonferenza (arcogradi °,
arcoprimi ’ e arcosecondi ”);
3) il moto proprio in ascensione retta
mpα espresso in secondi di tempo s per anno;
4) il moto proprio in declinazione mpδ
in secondi di circonferenza ” per anno.
Come vedremo fra poco, a queste quattro coordinate
fondamentali si possono aggiungere, per ottenere una maggiore precisione di
calcolo dei moti propri, le seguenti:
5) la distanza in parsec[2]
r;
6) la velocità radiale Δr in parsec per
anno col segno proprio + o -.
Queste coordinate si possono reperire sui cataloghi
stellari, la maggior parte dei quali sono oggi rintracciabili su Internet. Uno,
con note storiche, è visibile sul sito:
http://www.daviddarling.info/encyclopedia/A/astronomical_catalogues.html.
Il migliore servizio di banca-dati astronomici è sul sito
http://cdsweb.u-strasbg.fr/
del Centre de Donées Astronomiques de Strasbourg. Da qui si può accedere ai tre
servizi Sinbad, VizieR e Aladin, nonché ad un ricchissimo
elenco di cataloghi e banche dati (http://vizier.u-strasbg.fr/cats/cats.html).
Si tratta però di un sito non facilmente gestibile.
Due siti utili sono: quello dello United States Naval
Observatory U.S.N.O. (http://ad.usno.navy.mil/star/)
e quello dello Smithsonian Astrophysical Observatory S.A.O. (http://heasarc.gsfc.nasa.gov/W3Browse/star-catalog/sao.html,
poi link a sinistra in alto: browse this table).
Un altro sito di facile consultazione è http://www.alcyone.de/SIT/bsc/index.html
dove, nel Bright Star Catalogue 5th ed. 1991, si troverà l’elenco
di tutte le 88 costellazioni in cui la International Astronomical Union I.A.U.
suddivise il cielo nel 1933. Entrando in quella che interessa, si ottiene un
elenco delle stelle, numerate secondo vari cataloghi, che ne fanno parte e
cliccando sul numero all’estrema sinistra si apre una pagina con i link di
alcuni di essi. Cliccando su quello dello S.A.O. si otterranno i parametri di
quella particolare stella sia per il 2000 che per il 1950, mentre cliccando
sugli altri si otterranno i parametri del 2000 e del 1900.
Nel calcolo della posizione di una stella si prendono
in considerazione tre date (espresse in giorni giuliani delle effemeridi[3]
JD):
1) la data di partenza, ossia la data a partire dalla
quale si vogliono effettuare i calcoli (per es. la data odierna);
2) la data dell’epoca standard
prescelta (B1900.0, B1950.0, J2000.0, J2050.0, ecc.);
3) la data di arrivo, ossia la data, più o meno lontana
nel tempo, per la quale si vogliono calcolare le coordinate della stella.
La lettera B o J premessa all’anno
indica che quest’ultimo è, rispettivamente, l’anno
besseliano o l’anno giuliano. L’anno besseliano,
utilizzato fino al 1984 e di lunghezza uguale all’anno tropico
(pari a 365,2421988 giorni nel 1900 d. C., secondo Newcomb), incominciava nell’istante
in cui la longitudine media del Sole, affetta dall’aberrazione
(pari a -0°00’20,5”) e misurata dall’equinozio medio
della data, era esattamente 280°00’00,00”. Tale istante è
sempre molto vicino all’inizio dell’anno civile gregoriano.
I cataloghi stellari relativi agli anni B1900.0 e B1950.0 appartengono al
Fundamental Katalog 4 o FK4.
L’anno giuliano cominciò ad essere
utilizzato a partire dal 1984 per decisione dell’International Astronomical
Union UAI. Dura 365,25 giorni esatti e comincia sempre alle ore UT 12:00:00 del
01/01 di ciascun anno. I cataloghi stellari relativi all’anno
J2000.0 appartengono al Fundamental Katalog 5 o FK5.
Lo 0 dopo il punto nella numerazione dell’anno
indica l’inizio dell’anno stesso, sia che si
tratti di anno besseliano che di anno giuliano. Espresse come sopra ed in
giorni giuliani, le date d’inizio degli anni B1900.0, B1950.0,
J2000.0 e J2050.0 sono, rispettivamente: JD 2415020,3135 (= 01/01/1900 UT 19:31:26);
JD 2433282,4235 (= 01/01/1950 UT 22:09:50); JD 2451545,00 (= 01/01/2000 UT
12:00:00); JD 2469807,00 (= 01/01/2050 UT 12:00:00).
Oltre all’uso dell’anno
besseliano o giuliano esistono altre differenze tra FK4 e FK5 ed è possibile
utilizzare l’algoritmo di calcolo previsto per il sistema FK5 con
le coordinate del sistema FK4 apportando tre correzioni (Meeus 2005, pp.
139-140).
In questo articolo prenderemo in considerazione le
coordinate della sola epoca standard dell’anno B1950.0, che
appartengono al Fundamental Katalog 4 o FK4.
Il metodo di seguito descritto è quello classico di
Newcomb, preciso anche per distanze di tempo relativamente lunghe, ma comunque
non utilizzabile per più di alcuni millenni dal presente (per es., dà risultati
assolutamente errati per 30000 anni dal presente!).
Precessione degli equinozi e moto proprio danno la posizione
media delle stelle.
Poiché invece in Archeoastronomia generalmente occorre
conoscerne la posizione apparente, è necessario calcolare i seguenti
altri effetti:
1) la nutazione;
2) l’aberrazione annua della luce;
3) la parallasse annua.
La nutazione è un piccolo movimento periodico dell’equatore
celeste dovuto alla forza gravitazionale della Luna.
L’aberrazione annua è un fenomeno
relativistico della luce che sposta l’immagine della stella.
La parallasse annua, sempre inferiore a 0°00’00,8”
per tutte le stelle visibili ad occhio nudo tranne tredici (Meeus 1988 p. 72 e
2005, p. 150), può essere trascurata nei calcoli archeoastronomici.
Per la dettagliata descrizione di questi tre fenomeni
si rinvia a Smart 1977, capp. VIII-X, ed a Zagar 1984, capp. VII-IX.
In conclusione, la procedura per calcolare le
coordinate di una stella in un tempo diverso dall’attuale sono,
nell’ordine di esecuzione, le seguenti:
1) calcolo del moto proprio;
2) calcolo della precessione;
3) calcolo della nutazione;
4) calcolo dell’aberrazione annua.
Di seguito si descrive l’algoritmo
completo per l’epoca standard B1950.0 ed il sistema FK4 secondo Meeus
1988 e 1990, capp. 14, 15, 16 e 18 con alcune modifiche, rinviando ad altri
prossimi lavori la descrizione degli algoritmi FK4 B1900.0 e FK5 J2000.0.
Presi dagli almanacchi cartacei o sul web α
e δ della stella per il 1950, si trasformano in un’unica
unità di misura, in genere in gradi sessagesimali moltiplicando per 15 la
misura in tempo di α.
Poi si calcola la differenza di tempo T in
secoli tropici tra la data dell’epoca standard B1950.0 e la data di
arrivo JD:
T = (2433282,4235[4]
- JD) / 36524,2199[5].
Poi si calcolano le variazioni di α e δ
per effetto dei moti propri. Per fare ciò innanzitutto si moltiplica la
differenza di tempo T in secoli tropici per 100 e si ottiene la
differenza di tempo in anni tropici, si moltiplica tale differenza di tempo in
anni tropici per i moti propri della stella in α e δ e
si sommano a questi prodotti i valori iniziali di α e δ
dati dagli almanacchi per l’epoca standard del B1950.0:
α0 = α + (T * 100) * mpα
δ0 = δ + (T * 100) * mpδ.
Se gli almanacchi consultati riportano anche la
distanza in parsec r e la velocità radiale Δr in parsec per
anno[6]
col segno proprio + o -, si possono calcolare con precisione maggiore,
soprattutto per notevoli distanze di tempo, gli effetti dei moti propri della
stella, procedendo nel seguente modo[7]:
α = ascensione retta all’epoca iniziale
δ = declinazione all’epoca iniziale
r = distanza in parsec
Δr = velocità radiale in parsec
all’anno
mpαr = moto
proprio α radianti all’anno
mpδr = moto
proprio δ radianti all’anno
t = numero di anni
dall’epoca di partenza a quella di arrivo, negativo nel passato e positivo nel
futuro[8]
x = r * cos δ *
cos a
y = r * cos δ *
sen a
z = r * sen δ
Δx = (x / r) * Δr -
z * mpδr * cos α -
y * mpαr
Δy = (y / r) * Δr -
z * mpδr * sen α +
x * mpαr
Δz = (z / r) * Δr +
r * mpδr * cos δ
x1 = x + t * Δx
y1 = y + t * Δy
z1 = z + t * Δz
tan α0 = y1
/ x1
tan δ0 = z1
/ Ö(x12
+ y12).
Si procede poi a calcolare gli effetti della precessione,
della nutazione e dell’aberrazione annua della luce come già descritto sopra.
Poi si calcola l’effetto della precessione
degli equinozi con le formule rigorose di Newcomb:
ζ = 0°00’2304,948” * T + 0°00’00,302”
* T² + 0°00’00,018” * T³;
z = 0°00’2304,948” * T + 0°00’01,093”
* T² + 0°00’00,019” * T³;
θ = 0°00’2004,255” * T - 0°00’00,426”
* T² - 0°00’00,042” * T³.
Poi si calcolano:
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.
Poiché la tangente di un angolo è la stessa
rispettivamente nei quadranti I e III nonché II e IV, per calcolare tan (α1
- z) e collocarla nel quadrante esatto si può procede in due modi:
1) se il denominatore di A / B è minore di 0 (B <
0) al risultato di A / B si aggiungono 180°; se invece esso è maggiore di 0 (B
> 0) non si aggiunge nulla e il risultato è già la tangente nell’angolo
corretto;
2) si calcolano le coordinate polari di B / A (occorre
invertire i fattori rispetto alla formula originaria A / B), ottenendone due
risultati, il secondo dei quali è quello cercato. Le calcolatrici scientifiche
hanno un apposito tasto che trasforma le coordinate rettangolari in polari: pol
(B; A) e due specifiche memorie dove sono immagazzinati i due risultati[9].
Se tale secondo risultato fosse negativo, lo si somma algebricamente a 360°,
ottenendo così il corretto valore positivo.
Es.:
1a) tan (α1 - z) = 1,0405017
/ -0,4315299 = -67,47462466 + 180° = 112,5253753° = 112°31’31,3”;
1b) tan (α1 - z) = 1,0405017 / -0,4315299 = pol (-0,4315299;
1,0405017) = 112,5253753° = 112°31’31,3”;
2a) tan (α1 - z) = 2,456 / 1,852 =
52,98108896° = 52°58’51,92”;
2b) tan (α1 - z) = pol (1,852; 2,456)
= 52,98108896° = 52°58’51,92”.
Il risultato di tan (α1 - z) si
aggiunge a z che è già noto. Il risultato è l’ascensione
retta cercata.
Invece il calcolo sin δ1 = C non
richiede alcuna trasformazione, essendo il risultato già la declinazione
cercata.
Si ottengono così l’ascensione retta α1
e la declinazione δ1 corrette per i moti propri e
per la precessione degli equinozi.
Ora si calcola l’effetto della nutazione.
Ottenuta la differenza di tempo in secoli giuliani dal
1950 con la formula
T = (JD - 2415020,0) / 36525[10]
si calcolano i seguenti parametri:
1) longitudine
media del Sole L = 279,6967° + 36000,7689° * T + 0,000303° * T²;
2) longitudine
media della Luna L1 = 270,4342° + 481267,8831° * T - 0,001133° * T²;
3) anomalia
media del Sole M = 358,4758° + 35999,0498° * T - 0,000150° * T²;
4) anomalia
media della Luna M1 = 296,1046° + 477198,8491° * T + 0,009192° * T²;
5) longitudine
del nodo ascendente della Luna Ω = 259,1833° - 1934,1420° * T + 0,002078°
* T².
Ora si calcolano i valori della nutazione in
longitudine Δψ e della
nutazione in obliquità Δε, i cui coefficienti sono qui scritti, per risparmio
di spazio, in secondi sessagesimali con decimali (e nel calcolo devono essere
scritti in forma completa; per es.: 0,01737” = 0°00’00,01737”):
Δψ = - (17,2327” + 0,01737” * T) * sen Ω
- (1,2729”
+ 0,00013”
* T) * sen (2L) + 0,2088” * sen (2Ω) - 0,2037” * sen (2L1) +
(0,1261” - 0,00031” * T) * sen M +
0,0675” * sen M1
- (0,0497”
- 0,00012”
* T) * sen (2L + M) - 0,0342” * sen (2L1 - Ω) - 0,0261” * sen (2L1
+ M1) + 0,0214” * sen (2L - M) - 0,0149” * sen (2L - 2L1 +
M1) + 0,0124” * sen (2L - Ω) + 0,0114” * sen (2L1 - M1)
Δε = +
(9,21” + 0,00091” * T) * cos Ω + (0,5522” - 0,00029” * T) * cos (2L) -
0,0904” * cos (2Ω) + 0,0884” * cos 2L1 + 0,0216” * cos (2L + M)
+ 0,0183” * cos (2L1 - Ω) + 0,0113” * cos (2L1 + M1)
- 0,0093” * cos (2L - M) - 0,0066” * cos (2L - Ω).
Ora si calcolano le variazioni in ascensione retta α2
e in declinazione δ2 per effetto della nutazione. Per
fare ciò è necessario prima calcolare l’obliquità dell’eclittica
ε con la formula di Laskar:
U = T / 100
ε = 23°26’21,448” - 0°00’4680,93” * (T / 100) - 0°00’01,55” * (T / 100)² + 0°00’1999,25” * (T / 100)³ - 0°00’51,38” * (T / 100)4 - 0°00’249,67” * (T / 100)5 - 0°00’39,05” * (T / 100)6 + 0°00’07,12” * (T / 100)7 + 0°00’27,87” * (T / 100)8 + 0°00’05,79” * (T / 100)9 + 0°00’02,45” * (T / 100)10
poi si risolvono le seguenti formule:
α2 = (cos ε + sen
ε * sen α1 * tan δ1) * Δψ - (cos
α1 * tan δ1) * Δε
δ2 =
(sen ε *
cos α1) *
Δψ +
(sen α1) *
Δε.
Ora si calcola longitudine vera del Sole (dal
01/01/1900 UT 12:00:00):
T = (JD - 2415020,0) / 36525[11]
1) longitudine media geometrica del Sole L = 279,69668° +
36000,76892° * T + 0,0003025° * T²;
2) anomalia media del Sole M = 358,47583° + 35999,04975° * T -
0,00015 * T² - 0,0000033° * T³;[12]
3) equazione del centro del Sole C = + (1,91946° - 0,004789° *
T - 0,000014° * T²) * sen M + (0,020094° - 0,0001° * T) * sen (2M) + 0,000293°
* sen (3M);
4) longitudine vera del Sole Lv = L + C.
Ora si calcolano le variazioni in ascensione retta Δα3
e declinazione Δδ3 per effetto dell’aberrazione
annua della luce (ε è ancora l’obliquità dell’eclittica
calcolata con la formula di Laskar):
α3=
-0°00’20,49” * [(cos α1 *
cos L1 * cos ε +
sen α1 *
sen L1) / cos δ1]
δ3 =
-0°00’20,49” * [cos L1 * cos ε * (tan ε * cos δ1 -
sen α1 *
sen δ1) +
cos α1 *
sen δ1 *
sen L1].
Ora si hanno tutti i parametri necessari per calcolare
la posizione apparente della stella α4 e δ4 all’epoca
voluta sommando algebricamente le correzioni, rispettivamente:
α4 = α1 +
α2 + α3
δ4 = δ1 + δ2 + δ3.
Esempio numerico
Si vuole calcolare la posizione apparente di
Spica (α Virginis) all’equinozio di primavera
del 350 d. C..
Con uno dei metodi descritti in Meeus 1990, cap. 20,
od in Meeus 2005, cap. 27 si calcola la data dell’equinozio di
primavera dell’anno 350 d. C., corrispondente al 20/03/350 d. C. ore
13h 00m 17s UT.
Dal sito http://www.alcyone.de/SIT/bsc/index.html
si ricava che le coordinate di Spica nel 1950.0, secondo il catalogo S.A.O.,
sono:
α = 13h 22m 33,301s
δ = -10°54’03,36”
mpα = -0h 00m
00,0029s
mpδ = -0°00’00,033”.
Poiché l’ascensione retta α
è data in unità di tempo, la si riduce a gradi sessagesimali moltiplicandola
per 15:
α 13h 22m 33,301s * 15 = 200°38’19,51”.
Col metodo descritto in Meeus 1990, cap. 3, Meeus
2005, cap. 7, e Codebò 2010, si calcola il giorno giuliano JD del
20/03/350 d. C.: ore 13h 00,28m ottenendo JD 1848974,04186.
Poi si calcola:
T = (1848974,04186 - 2433282,4235) / 36524,2199 =
-15,9978333073
α0 = α + [(T *
100) * mpα] = 200,658084882
δ0
= δ + [(T * 100) * mpδ] = -10,8862686528
ζ
= 0°00’2304,948” * T + 0°00’00,302” * T² +
0°00’00,018” * T³ = -10,2418280209
z
= 0°00’2304,948” * T + 0°00’01,093” * T² +
0°00’00,019” * T³ = -10,1867316808
θ
= 0°00’2004,255” * T - 0°00’00,426” * T² -
0°00’00,042” * T³ = -8,88911159455
A
= cos δ0 * sin (α0 + ζ) = -0,177544568207
B
= cos θ * cos δ0 * cos (α0 + ζ) - sin θ *
sin δ0 = -0,983403711782
C
= sin θ * cos δ0 * cos (α0 + ζ) + cos
θ * sin δ0 = -0,0373505816924
tan (α1- z) = 180,047256593
sin δ1 = C = -2,14052858662
α1 = 180,047256593
δ1 = -2,14052858662
T = (JD - 2415020,0) / 36525 =
-15,4974937205
L = 279,6967° + 36000,7689° * T + 0,000303° * T² =
-557641,920487
L1 = 270,4342° + 481267,8831° * T -
0,001133° * T² = -7458175,83412
M = 358,4758° + 35999,0498° * T - 0,000150° * T² =
-557536,608444
M1 = 296,1046° + 477198,8491° * T +
0,009192° * T² = -7395087,85508
Ω = 259,1833° - 1934,1420°
* T + 0,002078° * T² = 30234,0358776
Δψ = - (17,2327” + 0,01737”
* T) * sen Ω - (1,2729” + 0,00013” * T) * sen
(2L) + 0,2088” * sen (2Ω) - 0,2037” * sen (2L1)
+ (0,1261” - 0,00031” * T) * sen M + 0,0675” * sen M1 - (0,0497” -
0,00012” * T) * sen (2L + M) - 0,0342” * sen (2L1 - Ω) - 0,0261” * sen (2L1 + M1) +
0,0214” * sen (2L - M) - 0,0149” * sen (2L - 2L1 + M1) +
0,0124” * sen (2L - Ω) + 0,0114” * sen (2L1
- M1) = 0,0005968330633
Δε = + (9,21” + 0,00091” *
T) * cos Ω + (0,5522” - 0,00029” * T) * cos
(2L) - 0,0904” * cos (2Ω) + 0,0884” * cos 2L1
+ 0,0216” * cos (2L + M) + 0,0183” * cos (2L1 - Ω) + 0,0113” * cos (2L1 + M1) -
0,0093” * cos (2L - M) - 0,0066” * cos (2L - Ω) = 0,0026583762126
U = T / 100 = -0,154974937205
ε = 23°26’21,448” - 0°00’4680,93” * (T / 100) - 0°00’01,55” * (T / 100)² + 0°00’1999,25” * (T / 100)³ - 0°00’51,38” * (T / 100)4 -
0°00’249,67” * (T / 100)5 - 0°00’39,05” * (T / 100)6 + 0°00’07,12” * (T / 100)7 +
0°00’27,87” * (T / 100)8 + 0°00’05,79” * (T / 100)9 + 0°00’02,45” * (T / 100)10 =
23,6517218873
α2 = (cos ε + sen ε * sen α1 * tan δ1) * Δψ - (cos α1 * tan δ1) * Δε = 0,0000838428716
δ2 = (sen ε * cos α1) * Δψ + (sen α1) * Δε = -0,001161078902
L = 279,69668° +
36000,76892° * T + 0,0003025° * T² = -557641,920938
M = 358,47583° +
35999,04975° * T - 0,00015 * T² - 0,0000033° * T³ = -557536,595356
C = + (1,91946° -
0,004789° * T - 0,000014° * T²) * sen M + (0,020094° - 0,0001° * T) * sen (2M)
+ 0,000293° * sen (3M) = 1,92610755743
Lv = L +
C = -557639,99483
α3= -0°00’20,49” * [(cos α1 * cos L1 * cos ε + sen α1 * sen L1)
/ cos δ1] = 0,0049687641091
δ3 = -0°00’20,49” * [cos L1
* cos ε * (tan ε * cos δ1 - sen α1 * sen δ1) + cos α1 * sen δ1* sen L1] =
-0,0018947049052
α4 = α1 + α2 + α3 = 180,0523092° =
12h 00m 12,55s
δ4 = δ1 + δ2 + δ3 = -2,14358437043°
= -2°08’36,9”.
Risulta quindi che Spica in data
20/03/350 d. C., alle ore 13h 00m 17s UT,
aveva la posizione apparente α4 12h 00m 12,55s e δ4
-2°08’36,9”.
Altri esempi numerici
1) calcolare le coordinate equatoriali
di β Tauri al mezzogiorno del 01/01/4061 a. C. JD 238143,0
α 1950 = 5h 23m
07,71s
δ
1950 = 26°34’01,74”
mpα
= 0,0019s
mpδ
= -0,175”
α4 01/01/4061
a. C. = 23h 50m 53,37s
δ4
01/01/4061 a. C. = 2°10’48,9”
2) calcolare le coordinate equatoriali
di ζ Tauri al mezzogiorno del 01/01/4061 a. C. JD 238143,0
α
1950 = 5h 34m 39,263s
δ 1950
= 21°06’50”
mpα
= 0,0001s
mpδ =
-0,022”
α4
01/01/4061 a. C. = 0h 08m 58,93s
δ4
01/01/4061 a. C. = -2°12’36,4”
3) calcolare le coordinate equatoriali
di α Ophiuchi (Ras al Hague) al mezzogiorno del 01/01/3000
a. C. JD 625674,0
α
1950 = 17h 32m 36,696s
δ
1950 = 12°35’41,92”
mpα
= 0,008s
mpδ
= -0,227”
α4
01/01/3000 a. C. = 13h 50m 25,08s
δ4
01/01/3000 a. C. = 28°16’08,29”.
Algoritmo sintetico
Di seguito viene dato l’intero algoritmo con le due varianti
di calcolo degli effetti del moto proprio delle stelle.
1) Calcolo dei Moti Propri
