Astronomia Teorica e Meccanica Celeste
Since summer 2010 this website, which involves Astrodynamics
and Celestial Mechanics, has a limited
section for English Speaking Users (ESU), where new enhancements
of SOFA routines are shown. As usually done up to now, much worked
examples will be presented with their complete codes (in Fortran90
and C) alongwith the related outputs in text format.
To download SOFA routines/functions please take the last version
directly from its Official Site.
And now, get a glance to: ESU Main Page
Elementi orbitali di un corpo celeste con 3 osservazioni angolari in Python
E' il classico problema di Gauss, risolto calcolando la seconda distanza eliocentrica (r2) tramite la nota equazione di Lagrange, un polinomio di ottavo grado che ammette o una oppure tre radici di r2. Sotto sono riportati sorgente del programma e database di un corpo, che per maggiore praticità sono stati compressi qui.
# --- pgaus.py ------ Ott-2020 ----------
# lettura parametri osservazioni da un file esterno e
# caricamento nelle rispettive variabili floating
#
# UNA radice positiva dell'equaz.di Lagrange di ottavo grado,
# per le altre due _eventuali_ c'e' l'opzione manuale
# nell'intervallo [0.2 to 2.0]
#
from math import pi, sqrt, sin, tan, cos, acos, asin, atan
rad= pi/180
mu= 1.0 # Parametro gravitaz. SOLE --> UA^3/UT0^2
kgauss= 0.01720209895 # costante di Gauss
uv0 = 29.78469169 # veloc. media Terra in km/s
tol= 1e-10
#
nome= input('\n Nome file senza ext. (es. baade, ceres, arend, kowal): ')
ch= input(' Ricerca altre due radici Eq.Lagrange (S/ ): ')
if ch=='S':
z=input(' Valore prima approx. [try 0.2 to 2.0]: ')
z=float(z)
else:
z=10
#
with open (nome+'.txt', "r") as myfile:
data = myfile.read().splitlines()
for k in range(0,len(data),10):
tle= data[k]+'\n'+data[k+1]+'\n'+data[k+2]+'\n'+ \
data[k+3]+'\n'+data[k+4]+'\n'+data[k+5]+'\n'+data[k+6]+'\n'+ \
data[k+7]+'\n'+data[k+8]+'\n'+data[k+9]+'\n'
#print(tle)
#
# ---------- Prospetto dei dati --------------------
print('\n'+data[31]+'\n'+data[32]+'\n'+data[33]+'\n'+data[34]+'\n'+data[35]+'\n' \
+data[36]+'\n'+data[37]+'\n'+data[38])
#
print ("\n CALCOLO DI UN'ORBITA PRELIMINARE ELLITTICA CON IL METODO DI GAUSS")
t1= float(data[0]); #print ('\n t1=',t1)
t2= float(data[10]); #print (' t2=',t2)
t3= float(data[20]); #print (' t3=',t3)
#
X1= float(data[1]); #print ('\n X1=',X1)
Y1= float(data[2]); #print (' Y1=',Y1)
Z1= float(data[3]); #print (' Z1=',Z1); R1= sqrt(X1**2+Y1**2+Z1**2); print (' R1=',R1)
#
X2= float(data[11]); #print ('\n X2=',X2)
Y2= float(data[12]); #print (' Y2=',Y2)
Z2= float(data[13]); #print (' Z2=',Z2); R2= sqrt(X2**2+Y2**2+Z2**2); print (' R2=',R2)
#
X3= float(data[21]); #print ('\n X3=',X3)
Y3= float(data[22]); #print (' Y3=',Y3)
Z3= float(data[23]); #print (' Z3=',Z3); R3= sqrt(X3**2+Y3**2+Z3**2); print (' R3=',R3)
#
# ======================================================================
Ar1ore= float(data[4]); Ar1minuti= float(data[5]); Ar1secondi= float(data[6]);
A1= Ar1ore+Ar1minuti/60+Ar1secondi/3600; #print ('\n Asc_retta_1=',A1,'gradi = ',A1*rad,'rad')
#
Ar2ore= float(data[14]); Ar2minuti= float(data[15]); Ar2secondi= float(data[16]);
A2= Ar2ore+Ar2minuti/60+Ar2secondi/3600; #print (' Asc_retta_2=',A2,'gradi = ',A2*rad,'rad')
#
Ar3ore= float(data[24]); Ar3minuti= float(data[25]); Ar3secondi= float(data[26]);
A3= Ar3ore+Ar3minuti/60+Ar3secondi/3600; #print (' Asc_retta_3=',A3,'gradi = ',A3*rad,'rad')
# ======================================================================
#
De1gra= float(data[7]); De1primi= float(data[8]); De1arcsec= float(data[9]);
D1=De1gra+De1primi/60+De1arcsec/3600; #print ('\n Declinaz_1=',D1,'gradi = ',D1*rad,'rad')
#
De2gra= float(data[17]); De2primi= float(data[18]); De2arcsec= float(data[19]);
D2=De2gra+De2primi/60+De2arcsec/3600; #print (' Declinaz_2=',D2,'gradi = ',D2*rad,'rad')
#
De3gra= float(data[27]); De3primi= float(data[28]); De3arcsec= float(data[29]);
D3=De3gra+De3primi/60+De3arcsec/3600; #print (' Declinaz_3=',D3,'gradi = ',D3*rad,'rad')
# ======================================================================
eps= float(data[30]); #print ('\n eps_eclittica=',eps,'gradi = ',eps*rad,'rad')
#
def FunzioneIpergeometricaHANSEN(h):
gs4 = 0.0
gs3= 11/9*h/(1+11/9*h/(1+11/9*h/(1+gs4)))
gs2= 11/9*h/(1+11/9*h/(1+11/9*h/(1+gs3)))
gs1= 11/9*h/(1+11/9*h/(1+11/9*h/(1+gs2)))
gs0= 11/9*h/(1+11/9*h/(1+11/9*h/(1+gs1)))
yyy= 1+10/9*h/(1+gs0)
return yyy
def FX(x):
fx= x**8 + KT*x**6 + LT*x**3 + MT
return fx
def FPriX(x):
fprix= 8*x**7 + 6*KT*x**5 + 3*LT*x**2
return fprix
def FSecX(x):
fsecx= 56*x**6 + 30*KT*x**4 + 6*LT*x
return fsecx
# --------- INIZIO Calcoli ----------
# Conversione (AR,De)_iesime in rad
A1=A1*15*rad; A2=A2*15*rad; A3=A3*15*rad
D1=D1*rad; D2=D2*rad; D3=D3*rad
#
# print(A1,A2,A3); print(D1,D2,D3) # eliminare
#
# ---- Tempi canonici ---------
tau1=(t3-t2)*kgauss; tau2=(t3-t1)*kgauss; tau3=(t2-t1)*kgauss
#
# ---- Coseni direttori l,m,n -------
l1=cos(D1)*cos(A1); m1=cos(D1)*sin(A1); n1=sin(D1)
l2=cos(D2)*cos(A2); m2=cos(D2)*sin(A2); n2=sin(D2)
l3=cos(D3)*cos(A3); m3=cos(D3)*sin(A3); n3=sin(D3)
#
# ----- Altri parametri -----
# Procedure Iterazioner2GAUSS;
# ============================
cc=2*sqrt(2)/3
R2cosPsi2=l2*X2+m2*Y2+n2*Z2
R2=sqrt(X2**2+Y2**2+Z2**2)
#
a1=m1/l1*l2-m2; a2=m3-m1/l1*l3; a3=n1/l1*l2-n2; a4=n3-n1/l1*l3
b1=Y1-m1/l1*X1; b2=m1/l1*X2-Y2; b3=Y3-m1/l1*X3
b4=Z1-n1/l1*X1; b5=n1/l1*X2-Z2; b6=Z3-n1/l1*X3
#
Eg=a3-a4/a2*a1
F1=b4-a4/a2*b1; F2=b5-a4/a2*b2; F3=b6-a4/a2*b3
c1zero=tau1/tau2; c3zero=tau3/tau2
nu1= tau1*tau3/6*(1+c1zero)
nu3= tau1*tau3/6*(1+c3zero)
Ag=(F1*c1zero+F2+F3*c3zero)/Eg
Bg=(F1*nu1+F3*nu3)/Eg
KT=2*Ag*R2cosPsi2-Ag*Ag-R2**2; LT=2*Bg*R2cosPsi2-2*Ag*Bg; MT=-Bg**2
#
# Salta1:
n= 8 #{Costante iterativa di LAGUERRE; Se n=1 Metodo Laguerre coincide con Newton-Raphson}
#z= 10 # Valore iniziale dell'Unica RADICE
while True:
#
#print ('\n KT=',KT)
#print (' LT=',LT)
#print (' MT=',MT)
#
# { Calcolo dell'UNICA Radice }
F0= FX(z); F1= FPriX(z); F2= FSecX(z) # { f, f' e f" }
RAX= (n - 1)**2 * F1**2 - n * (n - 1) * F0 * F2
if RAX < 0:
RAX=0
RAD1= F1 - sqrt(RAX); RAD2= F1 + sqrt(RAX)
if abs(RAD1) > abs(RAD2):
RADI= RAD1
else:
RADI= RAD2
zit= z - n * F0 / RADI # { Iterazione di LAGUERRE }
zax= abs(zit - z); z= zit
if zax < 1e-10:
break
#
print('\n Distanza Eliocentrica r2')
print(' r2 = {:.10f}'.format(zit),' Residuo = {:.2e}'.format(FX(zit)))
#
#
r2=zit # Radice r2
kount=0
# ----------------------------------------
while True:
kount+=1
c1= c1zero+nu1/r2**3
c3= c3zero+nu3/r2**3
RO2= Ag+Bg/r2**3
c3ro3=(b1*c1+b2+b3*c3-a1*RO2)/a2; RO3=c3ro3/c3
c1ro1=(X1*c1-X2+X3*c3+l2*RO2-l3*c3ro3)/l1; RO1=c1ro1/c1
#
x1=RO1*l1-X1; y1=RO1*m1-Y1; z1=RO1*n1-Z1;
x2=RO2*l2-X2; y2=RO2*m2-Y2; z2=RO2*n2-Z2;
x3=RO3*l3-X3; y3=RO3*m3-Y3; z3=RO3*n3-Z3;
r1=sqrt(x1**2+y1**2+z1**2); r3=sqrt(x3**2+y3**2+z3**2);
#
K1=sqrt(r2*r3+x2*x3+y2*y3+z2*z3)
K2=sqrt(r1*r3+x1*x3+y1*y3+z1*z3)
K3=sqrt(r1*r2+x1*x2+y1*y2+z1*z2)
#
h1=tau1*tau1/(K1*K1*(cc*K1+r2+r3))
h2=tau2*tau2/(K2*K2*(cc*K2+r1+r3))
h3=tau3*tau3/(K3*K3*(cc*K3+r1+r2))
#
y1t=FunzioneIpergeometricaHANSEN(h1)
y2t=FunzioneIpergeometricaHANSEN(h2)
y3t=FunzioneIpergeometricaHANSEN(h3)
#
c1it=c1zero*y2t/y1t; c3it=c3zero*y2t/y3t
c1ab=abs(c1it-c1); c3ab=abs(c3it-c3)
#
print('\n iter = {:2}'.format(kount),'eps(c1)={:.2e}'.format(c1ab),' eps(c3)={:.2e}'.format(c3ab))
print(' nu1 = {:.8f}'.format(nu1),' nu3 = {:.8f}'.format(nu3))
print(' c1 = {:.8f}'.format(c1it),' c3 = {:.8f}'.format(c3it))
print(' r1 = {:.10f}'.format(r1),' r3 = {:.10f}'.format(r3),' (distanze eliocentriche)')
#
nu1=(c1it-c1zero)*r2**3
nu3=(c3it-c3zero)*r2**3
#
if c1ab < tol and c3ab < tol:
break
#
# --------- Calcolo degli elementi orbitali del corpo celeste
# tramite il teorema di GIBBS, noti i 3 vettori (r1,r2,r3)
# che individuano il piano della traiettoria -----------------
#
si=sin(eps*rad); co=cos(eps*rad) # sin(eps) & cos(eps) in rad
ye1= y1*co+z1*si; ze1=-y1*si+z1*co
ye2= y2*co+z2*si; ze2=-y2*si+z2*co
ye3= y3*co+z3*si; ze3=-y3*si+z3*co
y1=ye1; z1=ze1; y2=ye2; z2=ze2; y3=ye3; z3=ze3
# Calcolo del Vettore D
d1=(y1*z2-z1*y2)+(y2*z3-z2*y3)+(y3*z1-z3*y1)
d2=(z1*x2-x1*z2)+(z2*x3-x2*z3)+(z3*x1-x3*z1)
d3=(x1*y2-y1*x2)+(x2*y3-y2*x3)+(x3*y1-y3*x1)
d= sqrt(d1**2+d2**2+d3**2)
# Calcolo del Vettore N
nn1= r3*(y1*z2-z1*y2)+r1*(y2*z3-z2*y3)+r2*(y3*z1-z3*y1)
nn2= r3*(z1*x2-x1*z2)+r1*(z2*x3-x2*z3)+r2*(z3*x1-x3*z1)
nn3= r3*(x1*y2-y1*x2)+r1*(x2*y3-y2*x3)+r2*(x3*y1-y3*x1)
nn= sqrt(nn1**2+nn2**2+nn3**2)
# Parametro p
p=nn/d
# Calcolo del Vettore S
s1=(r2-r3)*x1+(r3-r1)*x2+(r1-r2)*x3
s2=(r2-r3)*y1+(r3-r1)*y2+(r1-r2)*y3
s3=(r2-r3)*z1+(r3-r1)*z2+(r1-r2)*z3
s= sqrt(s1**2+s2**2+s3**2)
# Eccentricita', Semiasse magg., Dist.Perielio -> e,a,q
e=s/d; a=p/(1-e**2); q=a*(1-e)
# Calcolo del Vettore B2
bb1=d2*z2-d3*y2; bb2=d3*x2-d1*z2; bb3=d1*y2-d2*x2
lb=sqrt(1/(d*nn))
# Componenti della VELOCITA' V2
V2x= lb*(bb1/r2+s1); V2y= lb*(bb2/r2+s2); V2z= lb*(bb3/r2+s3)
V2= sqrt(V2x**2+V2y**2+V2z**2)
# Prodotto Scalare r2*V2
sigma2= x2*V2x+y2*V2y+z2*V2z;
# Componenti Vettori h,n,e
hi=y2*V2z-z2*V2y; hj=z2*V2x-x2*V2z; hk=x2*V2y-y2*V2x
h=sqrt(hi**2+hj**2+hk**2)
ni=-hj; nj=hi; n=sqrt(ni**2+nj**2)
ei=V2y*hk-V2z*hj-x2/r2; ej=V2z*hi-V2x*hk-y2/r2; ek=V2x*hj-V2y*hi-z2/r2
e=sqrt(ei**2+ej**2+ek**2)
#
# { Parametri Orbitali i,ê,ì }
Incl=acos(hk/h)
Nodo=acos(ni/n)
if hi<0:
Nodo=2*pi-Nodo
#
Om=acos((ni*ei+nj*ej)/(n*e))
if ek<0:
Om=2*pi-Om
print('\n\t\t RISULTATI FINALI - Parametri Orbitali')
print('\t\t -------------------------------------')
print('\t\t Semiasse magg. a = {:.8f}'.format(a),' UA')
print("\t\t Eccentricita' e = {:.8f}".format(e))
print('\t\t Dist.Perielio q = {:.8f}'.format(q),' UA');
print('\t\t Inclinazione i = {:.5f}'.format(Incl/rad),' gradi')
print('\t\t Long-Nodo-Asc Om = {:.5f}'.format(Nodo/rad),' gradi')
print('\t\t Argom-Perielio w = {:.5f}'.format(Om/rad),' gradi')
print('\t\t ------------------------
# EOF: pgaus.py
Il database dei dati osservati è composto da 40 righe, come il seguente preso ad esempio (cometa Arend-Rigaux), filename: arend.txt
5
+0.9518223
-0.2559233
-0.1109740
04
37
59.14
+18
23
49.6
13
+0.9834527
-0.1321738
-0.0573197
05
02
03.43
+21
15
38.5
23
+0.9962662
+0.0257442
+0.0111556
05
33
17.04
+24
11
45.6
23.4457889
B1950.0 Cometa Arend-Rigaux 1950 VII
------------------------------------------------------------------------------------------
Riepilogo:
Data Xsole Ysole Zsole AR De
1) 05.00000/03/1978 +0.9518223 -0.2559233 -0.1109740 04 37 59.14 +18 23 49.6
2) 13.00000/03/1978 +0.9834527 -0.1321738 -0.0573197 05 02 03.43 +21 15 38.5
3) 23.00000/03/1978 +0.9962662 +0.0257442 +0.0111556 05 33 17.04 +24 11 45.6
------------------------------------------------------------------------------------------
Arend B1950.0 (ultima riga di chiusura del database, NON cancellarla)
Esempio file arend
Nome file senza ext. (es. baade, ceres, arend, kowal): arend
Ricerca altre due radici Eq.Lagrange (S/ ):
B1950.0 Cometa Arend-Rigaux 1950 VII
------------------------------------------------------------------------------------------
Riepilogo:
Data Xsole Ysole Zsole AR De
1) 05.00000/03/1978 +0.9518223 -0.2559233 -0.1109740 04 37 59.14 +18 23 49.6
2) 13.00000/03/1978 +0.9834527 -0.1321738 -0.0573197 05 02 03.43 +21 15 38.5
3) 23.00000/03/1978 +0.9962662 +0.0257442 +0.0111556 05 33 17.04 +24 11 45.6
------------------------------------------------------------------------------------------
CALCOLO DI UN'ORBITA PRELIMINARE ELLITTICA CON IL METODO DI GAUSS
Distanza Eliocentrica r2
r2 = 1.5007501213 Residuo = -1.78e-15
iter = 1 eps(c1)=1.31e-05 eps(c3)=1.91e-05
nu1 = 0.00613744 nu3 = 0.00569905
c1 = 0.55738439 c3 = 0.44611141
r1 = 1.4785381135 r3 = 1.5346881432 (distanze eliocentriche)
iter = 2 eps(c1)=1.93e-08 eps(c3)=9.06e-08
nu1 = 0.00618158 nu3 = 0.00563444
c1 = 0.55738437 c3 = 0.44611132
r1 = 1.4784821143 r3 = 1.5347692270 (distanze eliocentriche)
iter = 3 eps(c1)=4.67e-10 eps(c3)=1.27e-09
nu1 = 0.00618152 nu3 = 0.00563414
c1 = 0.55738437 c3 = 0.44611131
r1 = 1.4784815082 r3 = 1.5347702388 (distanze eliocentriche)
iter = 4 eps(c1)=7.46e-12 eps(c3)=1.99e-11
nu1 = 0.00618151 nu3 = 0.00563413
c1 = 0.55738437 c3 = 0.44611131
r1 = 1.4784814988 r3 = 1.5347702545 (distanze eliocentriche)
RISULTATI FINALI - Parametri Orbitali
-------------------------------------
Semiasse magg. a = 3.52929084 UA
Eccentricita' e = 0.59216705
Dist.Perielio q = 1.43936062 UA
Inclinazione i = 17.83534 gradi
Long-Nodo-Asc Om = 121.60581 gradi
Argom-Perielio w = 328.87542 gradi
-------------------------------------
Ricerca altre orbite facendo variare manualmente il valore iniziale della radice di Lagrange. Il corpo scelto ne ammette altre due, eccole:
1..
Ricerca altre due radici Eq.Lagrange (S/ ): S
Valore prima approx. [try 0.2 to 2.0]: 0.4
Distanza Eliocentrica r2
r2 = 0.9928255051 Residuo = 4.44e-16
RISULTATI FINALI - Parametri Orbitali
-------------------------------------
Semiasse magg. a = 0.99780904 UA
Eccentricita' e = 0.01302043
Dist.Perielio q = 0.98481274 UA
Inclinazione i = 0.20725 gradi
Long-Nodo-Asc Om = 176.62297 gradi
Argom-Perielio w = 286.07555 gradi
-------------------------------------
2..
Ricerca altre due radici Eq.Lagrange (S/ ): S
Valore prima approx. [try 0.2 to 2.0]: 1.3
Distanza Eliocentrica r2
r2 = 1.2352599413 Residuo = 0.00e+00
RISULTATI FINALI - Parametri Orbitali
-------------------------------------
Semiasse magg. a = 1.18597596 UA
Eccentricita' e = 0.06433398
Dist.Perielio q = 1.10967862 UA
Inclinazione i = 13.67509 gradi
Long-Nodo-Asc Om = 133.95317 gradi
Argom-Perielio w = 222.71975 gradi
-------------------------------------
(17-Ott-2020)
Dati usati da F.Gauss per il calcolo degli elementi orbitali di Cerere
Sono quelli visibili in questa figura, che ho elaborato con il vecchio
ma sempreverde programma in Pascal P_GAUSS_PAS.
Sopra, su sfondo blu, sono riportati alcuni screenshot del procedimento di calcolo dell'orbita effettuato con il suddetto programma, che sarebbe opportuno convertire in un linguaggio più moderno (Python).
Confronto con i risultati ottenuti dalla fonte citata, in perfetta armonia.
Traiettorie Interplanetarie in Python
(Set-2020). Riprendo volentieri questo argomento di astrodinamica, che ha costituito da sempre (circa un quarto di secolo!) un settore di ricerca che mi ha appassionato. Il relativo software, sviluppato in QBasic e Pascal, non è più leggibile dai moderni computer a 64-bit di clock, per cui ho pensato bene di convertirlo, piano piano, nel potentissimo linguaggio del pitone.
Come primo argomento ho scelto il Teorema di Lambert, che risolve qualsiasi traiettoria del sistema gravitazionale dato, noti i raggi vettore (r1,r2) iniziali e finali, l'angolo tra essi compreso e il tempo di volo per percorrere il tragitto dalla posizione P1 a P2. E' di largo uso in astronautica quando si studia la traiettaria per i voli interplanetari di una sonda che si muove nel Sistema Solare.
Primo metodo: Battin
La teoria è trattata in questo capitolo
del libro Richard H. Battin - An Introduction to the Mathematics and Methods of Astrodynamics(1999)
da cui ho tratto l'esempio di figura con il sistema delle due equazioni non-lineari
da applicare.
# ------ battin_xyz.py --------- 12.09.2020 ------------
# Fonte dati: Problem 7-12 del libro di Battin, pag. 341
# Richard H. Battin - An Introduction to the Mathematics and Methods of Astrodynamics
# Revised Edition (Aiaa Education Series) (1999, Amer Inst of Aeronautics &)
#
# Dati del problema di Lambert:
# vettori (r1,r2) definiti con le rispettive componenti (x1,y1,z1), (x2,y2,z2)
# tempo di volo: t
# Calcolare: i parametri della traiettoria (a,p,e) e le due velocita' (V1,V2)
#
from math import pi, sqrt, sin, tan, cos, acos, asin, atan
# ------------
rad= pi/180
mu= 1.0 # Parametro gravitaz. SOLE --> UA^3/UT0^2
kgauss= 0.01720209895 # costante di Gauss
uv0 = 29.78469169 # veloc. media Terra in km/s
# ------------ Dati del problema (tratti da p.341)----------
# Coordinate dei vettori in AU, tempo in anni
x1 = 0.159321004; y1 = 0.579266185; z1 = 0.052359607 # componenti del vettore dato (r1)
x2 = 0.057594337; y2 = 0.605750797; z2 = 0.068345246 # componenti del vettore dato (r2)
t_anni = 0.010794065
# Risultati di Battin, solo Velocita' (v1) espressa in UA/year
v1x= -9.303603251; v1y= 3.018641330; v1z= 1.536362143
V1battin= sqrt(v1x**2 + v1y**2 + v1z**2)
coef_vel= 149597870.0/(365.25*86400) # coeff. di conversione da UA/year a km/s
tgiorni= t_anni*365.25 # tempo in giorni
t = tgiorni * kgauss
print ('\n V1_Battin= ',V1battin,'AU/y = ', V1battin*coef_vel, 'km/s')
print (' tempo_volo= ',t, 'giorni')
#
r1= sqrt(x1**2 + y1**2 + z1**2)
r2= sqrt(x2**2 + y2**2 + z2**2)
# Angolo tra i vettori (r1,r2)
teta= acos((x1*x2 + y1*y2 + z1*z2) /(r1*r2))
# --- Calcolo PARAMETRI Intermedi ---
c = sqrt(r1**2 + r2**2 - 2 * r1 * r2 * cos(teta)) # corda
s = (r1 + r2 + c) / 2 # semi-perimetro
lam = sqrt(r1 * r2) / s * cos(teta / 2) # lambda
Tadim = sqrt(8 * mu / s) * t / s
#
print ('\n r1=',r1,' r2= ',r2)
print(' Angolo tra i due vettori= ',teta, 'rad =', teta/rad, 'gradi')
print (' Lambda=',lam,' T=', Tadim)
#
cosf = lam * s * 2 / (r1 + r2); sinf = sqrt(1 - cosf**2)
tgf2 = (1 - cosf) / sinf
# Parametri Elle(l), rop, m
L = tgf2**2
rop = (r1 + r2 + 2 * lam * s) / 4
m = mu * t**2 / (8 * rop**3)
def f_csi(x):
eta = x / (1 + sqrt(1 + x))**2
gs4 = 0.0
gs3 = 169 / 25 / 27 * eta / (1 + 196 / 27 / 29 * eta / (1 + 225 / 29 / 31 * eta / (1 + gs4)))
gs2 = 100 / 19 / 21 * eta / (1 + 121 / 21 / 23 * eta / (1 + 144 / 23 / 25 * eta / (1 + gs3)))
gs1 = 49 / 13 / 15 * eta / (1 + 64 / 15 / 17 * eta / (1 + 81 / 17 / 19 * eta / (1 + gs2)))
gs0 = 16 / 7 / 9 * eta / (1 + 25 / 9 / 11 * eta / (1 + 36 / 11 / 13 * eta / (1 + gs1)))
csi = 8 * (1 + sqrt(1 + x)) / (3 + 1 / (5 + eta + 9 / 7 * eta / (1 + gs0))) # Formula (53)
return csi
def f_k(u):
fs4 = 0.0
fs3 = 928 / 3591 * u / (1 - 1054 / 4347 * u / (1 - 266 / 1035 * u / (1 - fs4)))
fs2 = 418 / 1755 * u / (1 - 598 / 2295 * u / (1 - 700 / 2907 * u / (1 - fs3)))
fs1 = 22 / 81 * u / (1 - 208 / 891 * u / (1 - 340 / 1287 * u / (1 - fs2)))
fs0 = 4 / 27 * u / (1 - 8 / 27 * u / (1 - 2 / 9 * u / (1 - fs1)))
ku = 1 / 3 / (1 - fs0) # Formula (58)
return ku
#
# -------- Costanti pre-ciclo ----------
kount=0 # contatore
eps= 1e-10 # tolleranza
# --------------------------------------
# ------ Successive approx. tramite il costrutto While True ------
x = L # valore iniziale di inizio iterazione
y = 1.0 # inizializzazione della variabile
x = sqrt(((1 - L) / 2)**2 + m / y**2) - (1 + L) / 2 # 1^ equaz. di GAUSS
print('')
while True:
kount += 1
# Param. ausiliario {eta} -> -1 < eta < +1 formula (52)
csi = f_csi(x)
# Parametri h1,h2
hden = (1 + 2 * x + L) * (4 * x + csi * (3 + x))
h1 = (L + x)**2 * (1 + 3 * x + csi) / hden; h2 = m * (x - L + csi) / hden
# Parametri ausiliari
B = 27 * h2 / (4 * (1 + h1)**3)
u = -B / 2 / (1 + sqrt(1 + B))
# Calcolo della Funzione Ipergeometrica K(u)
ku = f_k(u)
# 2^ Equaz. di Gauss -> Radice positiva y
yit = (1 + h1) / 3 * (2 + sqrt(1 + B) / (1 - 2 * u * ku**2))
#print()
print(' Iteraz= {:.0f}'.format(kount), ' --> y-iter= {:.8f}'.format(yit),\
' Residuo= {:.5e}'.format(abs(yit-y)))
# -------- parte finale -------
x = sqrt(((1 - L) / 2)**2 + m / y**2) - (1 + L) / 2 # 1^ equaz. di GAUSS
if abs(yit-y) > eps:
y= yit
else:
break # chiude il loop infinito di (while true) appena raggiunta la convergenza desiderata
#
#
# ------- Fine ciclo iterativo -----------
#
print('\n y-fin = {:.8f}'.format(yit), ' x-fin = {:.8f}'.format(x))
#
# Calcolo dei PARAMETRI ORBITALI {a,p,e}
a = m * s * (1 + lam)**2 / (8 * x * yit**2)
p = 2 * r1 * r2 * (y * (1 + x) * sin(teta / 2))**2 / (1 + lam)**2 / m / s
ecc = sqrt(1 - p / a)
# Calcolo delle VELOCITA' {V1, V2} - Metodo CLASSICO
# V1 = sqrt(2 * mu / r1 - mu / a); V2 = sqrt(2 * mu / r2 - mu / a)
# Calcolo delle VELOCITA' {V1, V2} tramite i parametri lagrangiani {f,g,f',g'}
# per verificare se i parametri lagrangiani di sotto sono CORRETTI
f = 1 - r2 / p * 2 * sin(teta / 2)**2
g = r1 * r2 * sin(teta) / sqrt(mu * p)
fpunto = (2 / p * sin(teta / 2)**2 - 1 / r1 - 1 / r2) * tan(teta / 2) * sqrt(mu / p)
gpunto = 1 - r1 / p * 2 * sin(teta / 2)**2
V1 = sqrt((r2**2 - 2 * f * r1 * r2 * cos(teta) + f**2 * r1**2) / g**2)
# Vettore V2 = fpunto * r1 + gpunto * V1
sigma1 = (r1 * r2 * cos(teta) - f * r1**2) / g
V2 = sqrt((fpunto * r1)**2 + (gpunto * V1)**2 + 2 * fpunto * gpunto * sigma1)
#
print ('\n ------ Risultati Traiettoria Interplanetaria -------')
print(" Semiasse maggiore: ",'a = {:.8f}'.format(a),'AU')
print(" Parametro focale: ",'p = {:.8f}'.format(p),'AU')
print(" Eccentricita': ",'e = {:.8f}'.format(ecc))
print('\n Lam = {:.8f}'.format(lam), ' T = {:.8f}'.format(Tadim))
print(" \n Velocita': ", ' V1 = {:.6f}'.format(V1),'AU/y = {:.3f}'.format(V1*uv0),'km/s')
print(' V2 = {:.6f}'.format(V2),'AU/y = {:.3f}'.format(V2*uv0),'km/s')
print(" \n Vel_Battin:", 'V1b = {:.6f}'.format(V1battin),'AU/y = {:.3f}'.format(V1battin*coef_vel),'km/s (in armonia con V1 calcolata sopra)')
print('\n ----------------- Fine Calcoli ----------------------')
#
# EOF: battin_xyz.py ------
'''
V1_Battin= 9.900991746569911 AU/y = 46.935358714681676 km/s
tempo_volo= 0.06781982972754777 giorni
r1= 0.6030539145426524 r2= 0.6123089158026407
Angolo tra i due vettori= 0.1745329249358501 rad = 9.999999984897816 gradi
Lambda= 0.9160269541241969 T= 0.35706944287290354
Iteraz= 1 --> y-iter= 1.00085682 Residuo= 8.56824e-04
Iteraz= 2 --> y-iter= 1.00085682 Residuo= 0.00000e+00
y-fin = 1.00085682 x-fin = 0.00065003
------ Risultati Traiettoria Interplanetaria -------
Semiasse maggiore: a = 1.20013342 AU
Parametro focale: p = 0.90003381 AU
Eccentricita': e = 0.50005520
Lam = 0.91602695 T = 0.35706944
Velocita': V1 = 1.575821 AU/y = 46.935 km/s
V2 = 1.559835 AU/y = 46.459 km/s
Vel_Battin: V1b = 9.900992 AU/y = 46.935 km/s (in armonia con V1 calcolata sopra)
----------------- Fine Calcoli ----------------------
'''
Secondo metodo: Bate
Teoria prelevata da questo paragrafo
del libro R. Bate, Donald D. Mueller, Jerry E. White - Fundamentals of Astrodynamics (1971, Dover)
insieme alla figura di riferimento.
# ------ lamb_bate_xyz.py --------- 29.09.2020 ------------
# stessi dati di battin_xyz.py
from math import pi, sqrt, sin, tan, cos, acos, asin, atan
# ------------
rad= pi/180
mu= 1.0 # Parametro gravitaz. SOLE --> UA^3/UT0^2
kgauss= 0.01720209895 # costante di Gauss
uv0 = 29.78469169 # veloc. media Terra in km/s
# ------------ Dati del problema (tratti da p.341 del Cap.7 Battin)----------
# Coordinate dei vettori in AU, tempo in anni
x1 = 0.159321004; y1 = 0.579266185; z1 = 0.052359607 # componenti del vettore dato (r1)
x2 = 0.057594337; y2 = 0.605750797; z2 = 0.068345246 # componenti del vettore dato (r2)
t_anni = 0.010794065
# Risultati di Battin, solo Velocita' (v1) espressa in UA/year
v1x= -9.303603251; v1y= 3.018641330; v1z= 1.536362143
V1battin= sqrt(v1x**2 + v1y**2 + v1z**2)
coef_vel= 149597870.0/(365.25*86400) # coeff. di conversione da UA/year a km/s
tgiorni= t_anni*365.25 # tempo in giorni
t = tgiorni * kgauss
print ('\n V1_Battin= ',V1battin,'AU/y = ', V1battin*coef_vel, 'km/s')
print (' tempo_volo= ',tgiorni, 'giorni')
#
r1= sqrt(x1**2 + y1**2 + z1**2)
r2= sqrt(x2**2 + y2**2 + z2**2)
# Angolo tra i vettori (r1,r2)
teta= acos((x1*x2 + y1*y2 + z1*z2) /(r1*r2)); tetagrad= teta/rad
#
print("\n PROBLEMA di LAMBERT risolto con le VARIABILI UNIVERSALI (Metodo BATE)")
print(" Determinazione degli elementi ORBITALI di una CONICA, NOTI:")
print(" 1 Due RAGGI Vettori r1,r2")
print(" 2 - L'angolo di trasferimento theta")
print(" 3 - Il tempo di Percorrenza tau=t2-t1")
#
print('\n Dati del problema')
print(' =================')
print(" Vettore r1: ",'r1 = {:.6f}'.format(r1),'AU')
print(" Vettore r2: ",'r2 = {:.6f}'.format(r2),'AU')
print(" Angolo tra vettori (r1,r2): ",'teta = {:.5f}'.format(tetagrad),'gradi')
print(" Tempo di volo: ",'t_volo = {:.3f}'.format(tgiorni),'giorni')
#
teta = tetagrad * rad; t = tgiorni * kgauss
if tetagrad<=180:
segno=1
else:
segno=-1
# --- Calcolo PARAMETRI Intermedi ---
amin=(r1+r2+sqrt(r1*r1+r2*r2-2*r1*r2*cos(teta)))/4
Ag=segno*sqrt(r1*r2*(1+cos(teta)))
#
# -------- Costanti pre-ciclo ----------
eps= 1e-10 # tolleranza
z=1.0 # Valore iniziale di z=x/a = Parametro adimensionale
#
# ----- Inizio ciclo: successive approx. tramite il costrutto While True ------
print('')
#
while True:
S=1/6*((((((((z/342-1)*z/272+1)*z/210-1)*z/156+1)*z/110-1)*z/72+1)*z/42-1)*z/20+1)
C=1/2*((((((((z/306-1)*z/240+1)*z/182-1)*z/132+1)*z/90-1)*z/56+1)*z/30-1)*z/12+1)
Spri=(C-3*S)/(2*z); Cpri=(1-z*S-2*C)/(2*z)
# Anomalia eccentrica x e parametro y
y=r1+r2-Ag*(1-z*S)/sqrt(C)
if y<0:
z=z/10
continue
x=sqrt(abs(y/C))
F=x*x*x*S+Ag*sqrt(abs(y))-t
F1=x*x*x*(Spri-3*S*Cpri/(2*C))+Ag/8 *(3*S*sqrt(abs(y))/C+Ag/x)
zit= z - F / F1 # Iterazione di NEWTON
print(' z-iter = {:.8f}'.format(zit), ' x-iter = {:.8f}'.format(x))
#
if abs(zit-z) > eps:
z= zit
else:
break # chiude il loop infinito di (while true) appena raggiunta la convergenza desiderata
# ------ Parametri FINALI della conica ------ (a,p,e) ---------
a=x*x/z; p=r1*r2/y*(1-cos(teta)); ecc=sqrt(1-p/a)
# Velocita' agli estremi dell'arco
V1=sqrt(2/r1-1/a); V2=sqrt(2/r2-1/a)
#
print ('\n ------ Risultati Traiettoria Interplanetaria -------')
print(" Semiasse maggiore: ",'a = {:.8f}'.format(a),'AU')
print(" Parametro focale: ",'p = {:.8f}'.format(p),'AU')
print(" Eccentricita': ",'e = {:.8f}'.format(ecc))
#
print(" \n Velocita': ", ' V1 = {:.3f}'.format(V1*uv0),'km/s;', ' V2= {:.3f}'.format(V2*uv0),'km/s')
# EOF: lamb_bate_xyz.py
# ------- OUTPUT ---------
'''
C:\Training>python lamb_bate_xyz.py
V1_Battin= 9.900991746569911 AU/y = 46.935358714681676 km/s
tempo_volo= 3.94253224125 giorni
PROBLEMA di LAMBERT risolto con le VARIABILI UNIVERSALI (Metodo BATE)
Determinazione degli elementi ORBITALI di una CONICA, NOTI:
1 Due RAGGI Vettori r1,r2
2 - L'angolo di trasferimento theta
3 - Il tempo di Percorrenza tau=t2-t1
Dati del problema
=================
Vettore r1: r1 = 0.603054 AU
Vettore r2: r2 = 0.612309 AU
Angolo tra vettori (r1,r2): teta = 10.00000 gradi
Tempo di volo: t_volo = 3.943 giorni
z-iter = -0.45385757 x-iter = 0.57666886
z-iter = 0.00874141 x-iter = 0.08912353
z-iter = 0.01037934 x-iter = 0.10942535
z-iter = 0.01039601 x-iter = 0.11167601
z-iter = 0.01039601 x-iter = 0.11169869
z-iter = 0.01039601 x-iter = 0.11169869
------ Risultati Traiettoria Interplanetaria -------
Semiasse maggiore: a = 1.20013342 AU
Parametro focale: p = 0.90003381 AU
Eccentricita': e = 0.50005520
Velocita': V1 = 46.935 km/s; V2= 46.459 km/s
'''
In conclusione, può essere utile avere gli stessi programmi python con i raggi vettori in modulo e l'aggiunta dell'angolo di tresferimento theta, che sopra veniva calcolato in automatico tramite le loro coordinate. Ecco il link.
... 29-Set-2020Astronomia da Almanacco
(Mag/Giu-2020). E' il classico problema inverso, tipico dei calendari astronomici,
in cui, assegnato un intervallo temporale, occorre stabilire quando avvengono alcuni
fenomeni tra cui i) Sorgere/Tramonto Pianeti,
ii) Equinozi/Solstizi, iii) Fasi Lunari, iv) Crepuscoli.
Alcuni linguaggi di programmazione, come il potentissimo Python, permettono di
raggiungere elevati standard di precisione nelle successive approssimazioni
richieste per risolvere tali problemi.
Nel primo caso, oltre ai due istanti in cui il pianeta sorge e tramonta, vengono
calcolati i loro azimut nord-est, l'altezza e l'orario al momento del transito.
Nel secondo caso, quando per gli equinozi
la longitudine del Sole Apparente sia nulla o pari a 180 gradi e per i solstizi
90/270 gradi; nel terzoo l'approssimazione rigurderà l'età
della Luna, ovvero (0, 0.25, 0.5, 0.75) del ciclo mensile, per le sue fasi
(nuova, primo quarto, piena, ultimo quarto); infine nel quarto caso,
i tre tipi di crepuscolo (civile, nautico, astronomico) si raggiungono quando
l'altezza del Sole Apparente si trova (-6, -12 -18) gradi dall'orizzonte.
Viene utilizzata la libreria (ephem) di questo pacchetto
che si installa con la nota routine pip install pyephem
Raccolta di alcuni script che utilizzano la libreria (ephem) in ambiente Python-2. Nello specifico, ecco quelli planetari
.Primo Esempio: Sorgere/Tramonto Pianeti
# UT_mars_ra_dec_az.py -------- May.2020 ---------------
# ------- Mars Ephemerides Table ---------
# -------- Python 3 environment ---------
# Windows: python UT_mars_ra_dec_az.py > out.txt
# Linux: python3 UT_mars_ra_dec_az.py > out.txt
import ephem
from datetime import datetime, timedelta
curtime = datetime(2002, 1, 5, 0, 0, 0) # start time
endtime = datetime(2002, 1, 30, 0, 0, 0) # end time
mars = ephem.Mars()
observer = ephem.Observer() # ficticious place nearby Rome
observer.lon = '12.0'
observer.lat = '42.0'
observer.elevation = 0
print('\n Observer Coords (Lon,Lat,H):', observer.lon, observer.lat, observer.elevation)
print(' --------------------------------------------------------')
print(' Celestial Body: Mars Instants in Universal Time\n')
print(' -----------------------------------------------------------------------------------------------------------------------------')
print(' Rising(UT) Transit(UT) Setting(UT) RA Dec Altit. Azimuth Azimuth')
print(' Transit Transit Transit Rising Setting')
print(' -----------------------------------------------------------------------------------------------------------------------------')
while curtime <= endtime:
observer.date = curtime
# ---------------------------------------
t2= observer.next_transit(mars); # Transit time, as reference for previous.rising and next.setting
observer.date = t2
t1= observer.previous_rising(mars);
azris= mars.az
observer.date = t2
t3= observer.next_setting(mars);
azset= mars.az
#----------------------------------------
observer.date = t2 # to compute RA, Dec, Altitude at Transit
mars.compute(observer)
raa = mars.ra
ddec = mars.dec
str1= str(ddec).zfill(11)
if str1[0] == '0':
str1=str1.replace('0','+',1)
altit= mars.alt
#-------------
str2= str(azris).zfill(11) # to compute Azimuth at Rising
str3= str(azset).zfill(11) # to compute Azimuth at Setting
#-----------
print ('',t1.datetime().strftime('%Y/%m/%d %H:%M:%S'),''\
,t2.datetime().strftime('%Y/%m/%d %H:%M:%S'),''\
,t3.datetime().strftime('%Y/%m/%d %H:%M:%S'),''\
,str(raa).zfill(11),'',str1,'',altit,'',str2,'',str3)
curtime += timedelta(days = 5) # timedelta parameters: days, hours, minutes, seconds
#
# EOF: UT_mars_ra_dec_az.py -----------
'''
Observer Coords (Lon,Lat,H): 12:00:00.0 42:00:00.0 0.0
--------------------------------------------------------
Celestial Body: Mars Instants in Universal Time
-----------------------------------------------------------------------------------------------------------------------------
Rising(UT) Transit(UT) Setting(UT) RA Dec Altit. Azimuth Azimuth
Transit Transit Transit Rising Setting
-----------------------------------------------------------------------------------------------------------------------------
2002/01/05 09:50:30 2002/01/05 15:37:01 2002/01/05 21:24:04 23:25:15.05 -04:24:41.0 43:36:18.9 095:31:59.8 264:39:53.6
2002/01/10 09:38:35 2002/01/10 15:30:37 2002/01/10 21:23:12 23:38:33.01 -02:52:58.3 45:07:58.5 093:28:34.7 266:43:32.2
2002/01/15 09:26:38 2002/01/15 15:24:11 2002/01/15 21:22:17 23:51:48.54 -01:21:03.1 46:39:50.8 091:24:57.5 268:47:20.7
2002/01/20 09:14:40 2002/01/20 15:17:43 2002/01/20 21:21:19 00:05:02.20 +00:10:43.9 48:11:35.0 089:21:31.6 270:50:55.6
2002/01/25 09:02:42 2002/01/25 15:11:14 2002/01/25 21:20:19 00:18:14.60 +01:42:02.5 49:42:50.9 087:18:39.8 272:53:53.8
2002/01/30 08:50:45 2002/01/30 15:04:44 2002/01/30 21:19:16 00:31:26.39 +03:12:32.5 51:13:18.4 085:16:45.2 274:55:52.2
'''
Secondo Esempio: Equinozi/Solstizi
# --------- equisol.py --------- Equinoxes and Solstices - Nov.2019 ----
# ----- Istanti in Temi Universali (UT) ------
# -------- Python 3 environment ---------
# Windows: python equisol.py
# Linux: python3 equisol.py
import ephem
d1 = ephem.next_equinox('2020')
print ('\n Equinozio=', d1)
d2 = ephem.next_solstice(d1)
print (' Solstizio=', d2)
d3 = ephem.next_equinox(d2)
print (' Equinozio=', d3)
d4 = ephem.next_solstice(d3)
print (' Solstizio=', d4)
#
d5= ephem.next_vernal_equinox('2021')
print ('\n Equinozio successivo= ', d5)
# EOF: -------- Equisol.py -------
'''
---------- OUTPUT --------------
Equinozio= 2020/3/20 03:49:33
Solstizio= 2020/6/20 21:43:47
Equinozio= 2020/9/22 13:30:38
Solstizio= 2020/12/21 10:02:09
Equinozio successivo= 2021/3/20 09:37:26
'''
Terzo esempio: Fasi Lunari
# fasi_luna.py -------- Jan.2020 ----------
# ----- Istanti in Temi Universali (UT) ------
# -------- Python 3 environment ---------
# Windows: python fasi_luna.py
# Linux: python3 fasi_luna.py
import ephem
from datetime import datetime, timedelta
curtime = datetime(2020, 1, 15, 0, 0, 0) # start time
endtime = datetime(2021, 1, 15, 0, 0, 0) # end time
moon = ephem.Moon()
print (' ------------------------------------------------------------------------------------------')
print (' Nuova Primo Quarto Piena Ultimo Quarto ')
print (' ------------------------------------------------------------------------------------------')
while curtime <= endtime:
d1= ephem.previous_new_moon(curtime)
d2= ephem.next_first_quarter_moon(d1)
d3= ephem.next_full_moon(d2)
d4= ephem.next_last_quarter_moon(d3)
print ('Fasi Luna -> ',d1.datetime().strftime('%Y/%m/%d %H:%M:%S'),' ',d2.datetime().strftime('%Y/%m/%d %H:%M:%S'),\
' ',d3.datetime().strftime('%Y/%m/%d %H:%M:%S'),' ',d4.datetime().strftime('%Y/%m/%d %H:%M:%S'))
curtime += timedelta(days = 30) # timedelta parameters: days, hours, minutes, seconds
'''
------------------------------------------------------------------------------------------
Nuova Primo Quarto Piena Ultimo Quarto
------------------------------------------------------------------------------------------
Fasi Luna -> 2019/12/26 05:13:06 2020/01/03 04:45:24 2020/01/10 19:21:17 2020/01/17 12:58:23
Fasi Luna -> 2020/01/24 21:41:58 2020/02/02 01:41:39 2020/02/09 07:33:15 2020/02/15 22:17:11
Fasi Luna -> 2020/02/23 15:31:59 2020/03/02 19:57:22 2020/03/09 17:47:42 2020/03/16 09:34:10
Fasi Luna -> 2020/03/24 09:28:11 2020/04/01 10:21:13 2020/04/08 02:35:03 2020/04/14 22:56:06
Fasi Luna -> 2020/04/23 02:25:49 2020/04/30 20:38:19 2020/05/07 10:45:11 2020/05/14 14:02:40
Fasi Luna -> 2020/05/22 17:38:49 2020/05/30 03:29:52 2020/06/05 19:12:20 2020/06/13 06:23:38
Fasi Luna -> 2020/06/21 06:41:24 2020/06/28 08:15:38 2020/07/05 04:44:22 2020/07/12 23:28:57
Fasi Luna -> 2020/07/20 17:32:54 2020/07/27 12:32:31 2020/08/03 15:58:43 2020/08/11 16:44:44
Fasi Luna -> 2020/08/19 02:41:37 2020/08/25 17:57:35 2020/09/02 05:22:01 2020/09/10 09:25:41
Fasi Luna -> 2020/09/17 11:00:10 2020/09/24 01:54:49 2020/10/01 21:05:13 2020/10/10 00:39:29
Fasi Luna -> 2020/10/16 19:31:00 2020/10/23 13:22:53 2020/10/31 14:49:07 2020/11/08 13:46:03
Fasi Luna -> 2020/11/15 05:07:09 2020/11/22 04:44:59 2020/11/30 09:29:39 2020/12/08 00:36:34
Fasi Luna -> 2020/12/14 16:16:33 2020/12/21 23:41:11 2020/12/30 03:28:11 2021/01/06 09:37:11
'''
Quarto esempio: Crepuscoli
# --------- crepus3.py ---------- Apr/Jul.2020 ----------
# Calcolo dei 3 tipi di crepuscolo: civile, nautico e astronomico
# oltre a sorgere, transito, tramonto Sole
# -------- Python 3 environment ---------
# Windows: python crepus3.py
# Linux: python3 crepus3.py
import ephem
#Make an observer
observer = ephem.Observer()
#PyEphem takes and returns only UTC times. Paris (France)
observer.date = "1961-08-22"
#Location of paris
observer.lon = str(2.330833) #Note that lon should be in string format
observer.lat = str(48.868889) #Note that lat should be in string format
print('Longitudine ', observer.lon ,'/ Latitudine ', observer.lat, '/ data ',observer.date ,'UT')
print( '')
#Elevation of paris in metres
observer.elev = 32.0
#To get U.S. Naval Astronomical Almanac values, use these settings
observer.pressure= 0
observer.horizon = '-0:34'
sunrise=observer.next_rising(ephem.Sun()) #Sunrise
noon =observer.next_transit (ephem.Sun()) #Solar noon , start=sunrise
sunset =observer.next_setting (ephem.Sun()) #Sunset
#We relocate the horizon to get twilight times
observer.horizon = '-18' #-6=civil twilight, -12=nautical, -18=astronomical
beg_twilight=observer.next_rising(ephem.Sun(), use_center=True) #Begin civil twilight
end_twilight=observer.next_setting (ephem.Sun(), use_center=True) #End civil twilight
end1=end_twilight
print(' Crepuscolo astronomico mattino = ',beg_twilight, 'UT')
#
observer.horizon = '-12' #-6=civil twilight, -12=nautical, -18=astronomical
beg_twilight=observer.next_rising(ephem.Sun(), use_center=True) #Begin civil twilight
end_twilight=observer.next_setting (ephem.Sun(), use_center=True) #End civil twilight
print(' Crepuscolo nautico mattino = ',beg_twilight, 'UT')
end2=end_twilight
#
observer.horizon = '-6' #-6=civil twilight, -12=nautical, -18=astronomical
beg_twilight=observer.next_rising(ephem.Sun(), use_center=True) #Begin civil twilight
end_twilight=observer.next_setting (ephem.Sun(), use_center=True) #End civil twilight
end3=end_twilight
#
print(' Crepuscolo civile mattino = ',beg_twilight, 'UT')
#
print( '')
print( ' Sole sorge = ', sunrise, 'UT')
print( ' Sole transita = ', noon, 'UT')
print( ' Sole tramonta = ', sunset, 'UT')
print( '')
print( ' Crepuscolo civile sera = ',end3, 'UT')
print( ' Crepuscolo nautico sera = ',end2, 'UT')
print( ' Crepuscolo astronomico sera = ',end1, 'UT')
#
# EOF: crepus3.py
'''
#----------------output ---------------------------------
Longitudine 2:19:51.0 / Latitudine 48:52:08.0 / data 1961/8/22 00:00:00 UT
Crepuscolo astronomico mattino = 1961/8/22 02:48:58 UT
Crepuscolo nautico mattino = 1961/8/22 03:36:18 UT
Crepuscolo civile mattino = 1961/8/22 04:18:31 UT
Sole sorge = 1961/8/22 04:52:35 UT
Sole transita = 1961/8/22 11:53:33 UT
Sole tramonta = 1961/8/22 18:53:32 UT
Crepuscolo civile sera = 1961/8/22 19:27:27 UT
Crepuscolo nautico sera = 1961/8/22 20:09:24 UT
Crepuscolo astronomico sera = 1961/8/22 20:56:14 UT
'''
La Teoria Planetaria VSOP87 con i Termini Periodici Completi
(Febbraio 2020). In passato erano disponibili solo quelli pubblicati da J.Meeus
in forma cartacea, questi consentono di avere maggiori approfondimenti come si può
notare scaricandoli insieme al programma principale (pmenu.exe) che richiama gli eseguibili
compilati con Fortran-90; il link è
il seguente.
Il codice ausiliario di avvio (12.exe) serve a trasformare la data dell'effemeride
in quella giuliana e a trasferirla su un file esterno (dummy.txt) che viene
letto dal menu' di scelta del pianeta. Un file ausiliario esterno (deg2dms.exe)
permette di trasformare gli angoli da decimali in gradi, primi e secondi.
C'è pure la versione Linux scaricabile da qui.
Lista dei Pianeti
------- MENU -------
1 Mercurio
2 Venere
3 Sole Apparente
4 Marte
5 Giove
6 Saturno
7 Urano
8 Nettuno
12 Scelta Data
--------------------
Pianeta nr. -> {20 per uscire} =
Data dell'Effemeride
Pianeta nr. -> {20 per uscire} = 12
Pianeta nr. 12
--------------------------------------------
Data del calendario e Ora (Tempi Dinamici)
--------------------------------------------
Immettere giorno (g): 20
mese (m): 2
anno (y): 2020
ore (ha): 17
minuti (mn): 25
secondi (sc): 0
--------------------------------------------
JDE= 2458900.2256944445
JDE trasferito su file esterno
cliccare su MENU-Pianeti.bat
per avviare il PROGRAMMA
----- Premere un tasto (+invio) per uscire ----
Scelta del Pianeta
Pianeta nr. -> {20 per uscire} = 2
Pianeta nr. 2
========= Prima parte ===========
Data Giuliana dell'Effemeride: 2458900.22569
Data Calendario: 2020- 2-20.72569
Pianeta VENERE - Solo POSIZIONE geometrica
L= Longitudine Eclittica Eliocentrica (Equin.Data)= 208.84129997 rad = 85.72507654 deg
B= Latitudine Eclittica Eliocentrica (Equin.Data)= 0.00913968 rad = 0.52366529 deg
R= Raggio Vettore Eliocentrico (Equin.Data)= 0.71992190 UA
--------------------------------------------------
Xecl= 0.05366238 Yecl= 0.71788899 Zecl= 0.00657977 | R= 0.71992190 [UA]
Pianeta TERRA - Solo POSIZIONE geometrica
L= Longitudine Eclittica Eliocentrica (Equin.Data)= 128.30861557 rad = 151.54214728 deg
B= Latitudine Eclittica Eliocentrica (Equin.Data)= -0.00000020 rad = -0.00001171 deg
R= Raggio Vettore Eliocentrico (Equin.Data)= 0.98872546 UA
--------------------------------------------------
Xecl=-0.86925566 Yecl= 0.47113971 Zecl=-0.00000020 | R= 0.98872546 [UA]
------------------------------------------------------------
Coordinate ECLITTICHE GEOcentriche del Pianeta (Equin.Data)
------------------------------------------------------------
Xgeo= 0.92291804 Ygeo= 0.24674928 Zgeo= 0.00657997 | dgeo= 0.95535659 [UA]
Long GEOC.= 374.96838 Latit GEOG.= 0.39462 Elong.(E/W)= 43.4 E
Ascensione retta= 0.90861 ore Declinaz.= 6.26034 deg
========= Seconda parte ===========
Data Giuliana dell'Effemeride: 2458900.22018 (considerato il tempo-luce)
Data Calendario: 2020- 2-20.72018
Pianeta VENERE - Con ABERRAZIONE e NUTAZIONE
L= Longitudine Eclittica Eliocentrica (Equin.Data)= 208.84114448 rad = 85.71616749 deg
B= Latitudine Eclittica Eliocentrica (Equin.Data)= 0.00913057 rad = 0.52314314 deg
R= Raggio Vettore Eliocentrico (Equin.Data)= 0.71992245 UA
--------------------------------------------------
Xecl= 0.05377405 Yecl= 0.71788124 Zecl= 0.00657321 | R= 0.71992245 [UA]
Pianeta TERRA - con ABERRAZIONE e NUTAZIONE
L= Longitudine Eclittica Eliocentrica (Equin.Data)= 128.30851845 rad = 151.53658283 deg
B= Latitudine Eclittica Eliocentrica (Equin.Data)= -0.00000021 rad = -0.00001189 deg
R= Raggio Vettore Eliocentrico (Equin.Data)= 0.98872425 UA
--------------------------------------------------
Xecl=-0.86920884 Yecl= 0.47122355 Zecl=-0.00000021 | R= 0.98872425 [UA]
------------------------------------------------------------
Coordinate ECLITTICHE GEOcentriche del Pianeta (Equin.Data)
------------------------------------------------------------
Xgeo= 0.92298289 Ygeo= 0.24665769 Zgeo= 0.00657342 | dgeo= 0.95539554 [UA]
Long GEOC.= 374.96207 Latit GEOG.= 0.39422 Elong.(E/W)= 43.4 E
Ascensione retta= 0.90823 ore Declinaz.= 6.25752 deg
----- Premere un tasto (+invio) per uscire -----
Alcuni file ausiliari (Python)
Riguardano la traformazione di gradi decimali in sessagesimali e la differenza (TT-UT), calcolata tramite la libreria pymeeus.
----------- (1) ---------------
# ------- deg2dms.py -------------
# Avvio: python deg2dms.py
# --------------------------------
def DecimaltoDMS(Decimal):
d = int(Decimal)
m = int((Decimal - d) * 60)
s = (Decimal - d - m/60) * 3600.00
z= round(s, 2)
if d >= 0:
print (" +", abs(d), "º", abs(m), "'", abs(z), '"')
else:
print (" -", abs(d), "º", abs(m), "'", abs(z), '"')
while True:
x = float(input("\n Angolo Decimale {0 per uscire}: "))
if x == 0:
break
DecimaltoDMS(x) # non serve il print() perche' dentro la funzione
# -------- OUTPUT -------
'''
Angolo Decimale {0 per uscire}: -36.89567
- 36 ° 53 ' 44.41 "
Angolo Decimale {0 per uscire}:
'''
----------- (2) ---------------
# --------- Delta-T.py -----------
# Avvio: python Delta-T.py
# --------------------------------
from pymeeus.Epoch import Epoch
print('\n Differenza tra Tempo Dinamico e Tempo Universale (TT-UT)')
print(' --------------------------------------------------------')
while True:
y = int(input("\n Anno {99 per uscire} (y): "))
if y == 99:
break
m = int(input(" Mese (m): "))
print("\n DeltaT (TT - UT)= ", round(Epoch.tt2ut(y, m), 1), 'secondi')
# ------ OUTPUT ----------
# Differenza tra Tempo Dinamico e Tempo Universale (TT-UT)
# --------------------------------------------------------
#
# Anno {99 per uscire} (y): 2045
# Mese (m): 3
#
# DeltaT (TT - UT)= 88.9 secondi
#
# Anno {99 per uscire} (y): 99
... (25-Feb-2020)
Eclissi di Sole senza Besseliani in Basic e Fortran
(Ottobre 2019). In questi giorni ho messo a punto il magnifico software che David Eagle
ha scritto tanto tempo fa in Basic e
reperibile qui.
Con gli attuali computer a 64bit di clock è stato necessario ricorrere al pacchetto
qb64, che, rispetto al vecchio QBASIC interpretato, permette di ricavare la versione compilata
in .EXE di ogni programma in Basic.
A dire il vero, lo stesso autore americano ha recentemente rilasciato
qui
una versione dello stesso proramma in MatLab/Octave, ma, come si evince dalle note inserite,
esso è affetto da errori.
La caratteristica peculiare del programa è la determinazione molto accurata dei
tempi dei 3 eventi nelle eclissi parziali di sole riferite ad una qualsiasi
località terrestre di note coordinate geografiche e dei tempi dei 5 eventi
per le eclissi totali.
Il metodo analitico di David Eagle
E' basato essenzialmente sul calcolo dei punti di tangenza dei cerchi
di Sole e Luna, che sono di tipo esterno per gli istanti di ingresso/uscita
del disco lunare (penumbra) e di tipo interno nella fase di totalità
(umbra). Vengono sfruttate due funzioni di minimizzazione delle distanze
relative sole-luna che l'astronomo americano ha sapientemente utilizzato
a partire da un ampio intervallo di giorni (Duration days= 30) intorno al
giorno selezionato dell'eclisse.
Ció consente di evitare l'utilizzo dei noti termini besseliani,
riducendo ad un problema geometrico la determinazione delle circostanze
principali delle eclissi di sole. Un abile artificio che consente grandi
accuratezze sui loro tempi (ingresso/uscita, fase massima), paragonati a
quelli rigorosi dei besseliani.
Elenco Eclissi di Sole nel periodo 1951-2050
Questa lista puó essere utile per le applicazioni del programma che vedremo più in basso
| July 10, 1972 | January 4, 1973 | June 30, 1973 | December 24, 1973 | June 20, 1974 | December 13, 1974 | | May 11, 1975 | November 3, 1975 | April 29, 1976 | October 23, 1976 | April 18, 1977 | October 12, 1977 | | April 7, 1978 | October 2, 1978 | February 26, 1979 | August 22, 1979 | February 16, 1980 | August 10, 1980 | | February 4, 1981 | July 31, 1981 | January 25, 1982 | June 21, 1982 | July 20, 1982 | December 15, 1982 | | June 11, 1983 | December 4, 1983 | May 30, 1984 | November 22, 1984 | May 19, 1985 | November 12, 1985 | | April 9, 1986 | October 3, 1986 | March 29, 1987 | September 23, 1987 | March 18, 1988 | September 11, 1988 | | March 7, 1989 | August 31, 1989 | January 26, 1990 | July 22, 1990 | January 15, 1991 | July 11, 1991 | | January 4, 1992 | June 30, 1992 | December 24, 1992 | May 21, 1993 | November 13, 1993 | May 10, 1994 | | November 3, 1994 | April 29, 1995 | October 24, 1995 | April 17, 1996 | October 12, 1996 | March 9, 1997 | | September 1, 1997 | February 26, 1998 | August 22, 1998 | February 16, 1999 | August 11, 1999 | February 5, 2000 | | July 1, 2000 | July 31, 2000 | December 25, 2000 | June 21, 2001 | December 14, 2001 | June 10, 2002 | | December 4, 2002 | May 31, 2003 | November 23, 2003 | April 19, 2004 | October 14, 2004 | April 8, 2005 | | October 3, 2005 | March 29, 2006 | September 22, 2006 | March 19, 2007 | September 11, 2007 | February 7, 2008 | | August 1, 2008 | January 26, 2009 | July 22, 2009 | January 15, 2010 | July 11, 2010 | January 4, 2011 | | June 1, 2011 | July 1, 2011 | November 25, 2011 | May 20, 2012 | November 13, 2012 | May 10, 2013 | | November 3, 2013 | April 29, 2014 | October 23, 2014 | March 20, 2015 | September 13, 2015 | March 9, 2016 | | September 1, 2016 | February 26, 2017 | August 21, 2017 | February 15, 2018 | July 13, 2018 | August 11, 2018 | | January 6, 2019 | July 2, 2019 | December 26, 2019 | June 21, 2020 | December 14, 2020 | June 10, 2021 | | December 4, 2021 | April 30, 2022 | October 25, 2022 | April 20, 2023 | October 14, 2023 | April 8, 2024 | | October 2, 2024 | March 29, 2025 | September 21, 2025 | February 17, 2026 | August 12, 2026 | February 6, 2027 | | August 2, 2027 | January 26, 2028 | July 22, 2028 | January 14, 2029 | June 12, 2029 | July 11, 2029 | | December 5, 2029 | June 1, 2030 | November 25, 2030 | May 21, 2031 | November 14, 2031 | May 9, 2032 | | November 3, 2032 | March 30, 2033 | September 23, 2033 | March 20, 2034 | September 12, 2034 | March 9, 2035 | | September 2, 2035 | February 27, 2036 | July 23, 2036 | August 21, 2036 | January 16, 2037 | July 13, 2037 | | January 5, 2038 | July 2, 2038 | December 26, 2038 | June 21, 2039 | December 15, 2039 | May 11, 2040 | | November 4, 2040 | April 30, 2041 | October 25, 2041 | April 20, 2042 | October 14, 2042 | April 9, 2043 | | October 3, 2043 | February 28, 2044 | August 23, 2044 | February 16, 2045 | August 12, 2045 | February 5, 2046 | | August 2, 2046 | January 26, 2047 | June 23, 2047 | July 22, 2047 | December 16, 2047 | June 11, 2048 | | December 5, 2048 | May 31, 2049 | November 25, 2049 | May 20, 2050 | November 14, 2050 | April 11, 2051 |
Sorgente del programma (seclful3.bas) ed eseguibile (seclful3.exe)
I due file sono zippati in questa cartella e in ambiente Windows l'eseguibile seclful3.exe, puó essere avviato direttamente da console. Due note importanti:
1.. Le coordinate geografiche possono essere date o in (gradi sessadecimali,0,0)
oppure tradizionalmente in (gradi, primi, secondi);
qb64 ammette l'input con soli (gradi sessadecimali)
2.. I risultati che compaiono sulle mascherine blu sono STAMPATI
anche sul file esterno ecl-out.txt
Eclisse parziale del 2027, località italiana
I dati completi sono riportati all'interno del file out-2027.txt della cartella dei risultati Qui sotto sono evidenziati input e tempi di parzialità e fase massima
Modalita' di inserimento dei dati
---------------------------------
Latitude
37° 01' 58" North
Longitude
15° 03' 54" East
Altitude (m)
370
Time zone ( 0 = Greenwich, West long negative )
1
Daylight Savings Time ( 1 = yes, 0 = no )
1
Initial calendar date
2027 8 2
Circostanze dell'eclisse parziale: inizio, fase massima, fine
---------------------------------
Enter Penumbra Shadow
Local civil time 10 hours 2.23 minutes
Maximum Eclipse Conditions
Local civil time 11 hours 17.76 minutes
Magnitude 0.927
Obscuration 92.36 %
Exit Penumbra Shadow
Local civil time 12 hours 36.71 minutes
Eclisse totale del 1983, località indonesiana
Come sopra, riportiamo i tempi dell'evento, che al completo sono conservati nel file out-1983.txt
Modalita' di inserimento dei dati
---------------------------------
Latitude
6° 14' 42" South
Longitude
114° 09' 54" East
Altitude (m)
0
Time zone ( 0 = Greenwich, West long negative )
7
Daylight Savings Time ( 1 = yes, 0 = no )
1
Initial calendar date
1983 6 11
Circostanze dell'eclisse totale: inizio parzialità, inizio totalità, fase massima, fine totalità, fine parzialità
-------------------------------
Enter Penumbra Shadow
Local civil time 11 hours 5.12 minutes
Enter Umbra Shadow
Local civil time 12 hours 41.11 minutes
Maximum Eclipse Conditions
Local civil time 12 hours 43.67 minutes
Magnitude 1.023
Obscuration 100.00 %
Exit Umbra Shadow
Local civil time 12 hours 46.22 minutes
Exit Penumbra Shadow
Local civil time 14 hours 22.94 minutes
Eclisse anulare del 2019, località Singapore
Qui l'anularità è completa, come si vede dai 5 tempi dell'evento riportati sotto e in dettaglio conservati nel file out-2019.txt
Modalita' di inserimento dei dati
---------------------------------
Latitude
1° 17' 25" North
Longitude
103° 51' 07" East
Altitude (m)
15
Time zone ( 0 = Greenwich, West long negative )
8
Daylight Savings Time ( 1 = yes, 0 = no )
0
Initial calendar date
2019 12 26
Circostanze dell'eclisse anulare:
--------------------------------
Enter Penumbra Shadow
Local civil time 11 hours 26.98 minutes
Enter Umbra Shadow
Local civil time 13 hours 22.52 minutes
Magnitude 0.969
Obscuration 94.15 %
Maximum Eclipse Conditions
Local civil time 13 hours 23.58 minutes
Magnitude 0.971
Obscuration 94.22 %
Exit Umbra Shadow
Local civil time 13 hours 24.63 minutes
Magnitude 0.969
Obscuration 94.16 %
Exit Penumbra Shadow
Local civil time 15 hours 18.50 minutes
E' interessante riportare, nel grafico seguente, il confronto tra
quanto dedotto da un software professionale e i risultati ottenuti
con questo programma.
La versione dell'algoritmo di David Eagle in Fortran
E' dovuta alla grande abilità e applicazione dell'astrofilo e amico Aldo Nicola, che ha certosinamente convertito il codice Basic in un linguaggio più moderno, il Fortran90, il cui compilatore (gfortran) è reperibile su tutti i repository delle distribuzioni Linux e sul web per gli utenti Windows. E con molto piacere ne ospito qui il suo contributo.
Sorgente del programma (seclipse.f90) e subroutine (111.f90)
I due file sono zippati in questa cartella
insieme all'eseguibile Windows, che si chiama a.exe, e agli output di alcuni esempi applicativI
in formato testo.
La compilazione avviene tramite il comando da terminale gfortran -w seclipse.f90 111.f90
I risultati delle 3 applicazioni viste sopra
1° esempio
----------
Latitude(degs)= 37.03280 Longitude(degs)= 15.06500 Altitude= 370.0 meters
Time zone= 1 [0 = Greenwich, West long negative]
DST= 1 [1= yes Daylight Savings Time; 0= no]
Initial Calendar Date: Year=2027 Month= 8 Day= 2
Enter Penumbra Shadow
Local civil time HR 10: MIN 2: SEC 13.9
Maximum Eclipse Conditions
Local civil time HR 11: MIN 17: SEC 45.7
Magnitude 0.927
Obscuration 92.36 %
Exit Penumbra Shadow
Local civil time HR 12: MIN 36: SEC 42.6
2° esempio
----------
Latitude(degs)= -6.24500 Longitude(degs)= 114.16500 Altitude= 0.0 meters
Time zone= 7 [0 = Greenwich, West long negative]
DST= 1 [1= yes Daylight Savings Time; 0= no]
Initial Calendar Date: Year=1983 Month= 6 Day=11
Enter Penumbra Shadow
Local civil time HR 11: MIN 5: SEC 7.4
Enter Umbra Shadow
Local civil time HR 12: MIN 41: SEC 6.8
Magnitude 1.001
Obscuration 100.00 %
Maximum Eclipse Conditions
Local civil time HR 12: MIN 43: SEC 40.2
Magnitude 1.023
Obscuration 100.00 %
Exit Umbra Shadow
Local civil time HR 12: MIN 46: SEC 13.5
Magnitude 1.000
Obscuration 100.00 %
Exit Penumbra Shadow
Local civil time HR 14: MIN 22: SEC 56.2
3° esempio
----------
Latitude(degs)= 1.29027 Longitude(degs)= 103.85196 Altitude= 15.0 meters
Time zone= 8 [0 = Greenwich, West long negative]
DST= 0 [1= yes Daylight Savings Time; 0= no]
Initial Calendar Date: Year=2019 Month=12 Day=26
Enter Penumbra Shadow
Local civil time HR 11: MIN 26: SEC 59.1
Enter Umbra Shadow
Local civil time HR 13: MIN 22: SEC 31.4
Magnitude 0.969
Obscuration 94.15 %
Maximum Eclipse Conditions
Local civil time HR 13: MIN 23: SEC 34.7
Magnitude 0.971
Obscuration 94.22%
Exit Umbra Shadow
Magnitude 0.969
Obscuration 94.16 %
Exit Penumbra Shadow
Local civil time HR 12: MIN 36: SEC 42.6
E' tutto!
... fine (18-10-2019), agg. fortran (22-10-2019)Le effemeridi planetarie di Aldo Nicola, astrofilo e amico
(Agosto 2019). La collaborazione intellettuale con Aldo
è iniziata nel 2012 quando egli fu incuriosito da questa sezione del sito
Standards of Fundamental Astronomy (SOFA)
e mi pose dei quesiti su come si potessero applicare le potentissime librerie
SOFA (scritte in Fortran) sulle posizioni e velocità di Pianeti, Luna e Sole Apparente.
Il problema era molto arduo, perché bisognava tener conto di svariati
fattori, ovvero l'introduzione di un nuovo modello di precessione-nutazione,
il complicatissimo moto istantaneo del Polo Terrestre fino alle correzioni relativistiche
e alle oscillazioni dovute alle maree e soprattutto il tipo di approccio mentale
nell'applicare il calcolo di pos/vel planetari completamente diverso dal precedente,
imposto dal nuovo modello.
Va da sé che i primi tentativi dell'attuale versione di Planet430B furono
fatti nel periodo 2014-15 ed io accettai volentieri la mansione di beta-tester
che è continuata senza interruzione. C'era inoltre un altro problema pratico
da superare, l'utilizzo di un enorme database binario (JPLEPH.430), da inglobare nella stessa
cartella, ma di difficile fruibilità nello scambio con altri astrofili.
La faccenda è stata risolta da poco grazie a questo
big repository,
che permette di effettuare facilmente i necessari e periodici aggiornamenti
di file di supporto e di scaricare in pochi secondi la cartella Planet430B-New
come indicato nella seguente figura:
Gli utenti Windows possono direttamente avviare l'eseguibile a-1919.exe,
mentre i linuxiani hanno necessità di scaricare il compilatore
gfortran presente in ogni distribuzione, effettuando la compilazione da terminale
con il comando indicato nel file 111_gfortran_command.txt e avviando quindi
il programma dalla stessa console.
Primo Esempio
L'effemeride scelta riguarda una data particolare [2027-08-02, 09:18:00 UTC], quella dell'eclisse solare del 2027, che nella località di riferimento [φ=37°.0327, λ=15.065°, H=370m] risulterà al limite della totalità [m=0.966, secondo Fred Espenak].Noi calcoleremo (Azimut, Altezza) di Sole e Luna, oltre ai rispettivi valori di (AR,Dec). Come descritto prima, si avvia l'eseguibile a-1919.exe e si inseriscono i dati nella sequenza richiesta, per semplicità fare riferimento al relativo file zippato qui. Esso contiene anche i tantissimi risultati che il programma mostra a video, quelli che a noi interessano si possono estrapolare e conservare su file. Eccoli:
Sun: righe 199-202
==================================================================================================
Apparent Topocentric position no refraction (IAU 2000 EQUINOX BASED)
R.A Dec. Azimut (*)Elevation . (*)Hour Angle
h , m , s o , ' , '' degs degs degs
8 49 19.25918 + 17 46 9.4191 121.717890 59.436393 332.981714
==================================================================================================
Moon: righe 405-408
==================================================================================================
Apparent Topocentric position no refraction (IAU 2000 EQUINOX BASED)
R.A Dec. Azimut (*)Elevation . (*)Hour Angle
h , m , s o , ' , '' degs degs degs
8 49 12.06536 + 17 43 5.3506 121.828904 59.420861 333.011689
==================================================================================================
Nota importante
I risultati della RA della Luna (il cui moto e' molto complesso per le miriadi di correzioni necessarie) risultano accurati fino a 0.0002 arcsec e quelli della Declinazione fino a 0.006 arcsec, rispetto a software professionali di altissimo livello. Strepitoso successo per un astrofilo
Secondo esempio
Anche qui facciamo riferimento a una data particolare, la grande opposizione di Marte del 2018, [2018-07-27, 05:07:11 UTC]. Input e output sono raccolti nel relativo file dello stesso zip.
Mars: righe 169-180
---------------------
Heliocentric rectangular position and velocity mean equator and equinox J2000.
X Y Z
POS 0.7802379594649730 -1.0483574299257052 -0.5019126451465140 AU
VEL 0.0121444498477254 0.0083032798082133 0.0034807019908563 AU/day
Heliocentric ecliptic coordinate mean equator and equinox J2000 ( Lam, B, RV )
o , ' , '' o , ' , '' RV AU Rdot (Km/s)
303 53 28.59772 - 1 46 47.9578 1.3999074532 -1.207511
Heliocentric ecliptic coordinate mean equator and equinox of the date ( Lam, B, RV )
o , ' , '' o , ' , '' RV AU Rdot (Km/s)
304 9 2.57541 - 1 46 54.7335 1.3999074532 -1.207511
------------------------
... fine (19-08-2019)
Python: linguaggio ideale per effemeridi satellitari
grazie alle librerie disponibili
(Maggio/Giugno 2019). Negli esempi in basso [#6] è stata utilizzata questa.
Primo esempio
# ---- eci_pv_ele.py ---------- Jun.30-Jul.01, 2019 ------
# Main source: https://bwinkel.github.io/pycraf/satellite/index.html
# Source for elements calcs: https://github.com/aerospaceresearch/orbitdeterminator [file state_kep.py]
from numpy import sqrt, arccos, pi, arctan2, dot, cross, clip, array
from math import acos, degrees
from astropy.coordinates import EarthLocation
from astropy import time
from pycraf import satellite
#
# ECI stands for Earth Center Inertial (frame)
mu = 398600.4405 # Earth gravitational parameter [km^3/s^2]
# --------- Step 0 -------------
# Given prediction datetime
year=2019; month=6; day=29; utch=11; utcm=3; utcs=57
print ('\n Prediction date: ', year, month, day )
print (' Prediction time: ', utch, utcm, utcs, ' UTC' )
print (' ================= ')
# --------- Step 1 -------------
with open ("3obj.txt", "r") as myfile:
data = myfile.read().splitlines()
for k in range(0,len(data),3):
tle= data[k]+'\n'+data[k+1]+'\n'+data[k+2]+'\n'
satname, sat = satellite.get_sat(tle)
print ('\n Name = ', satname,'\n Epoch= ', sat.epoch, '\n Numb = ', sat.satnum)
# --------- Step 2 ------------
# using sgp4 directly, to get position and velocity in ECI coordinates
r, v = sat.propagate(year, month, day, utch, utcm, utcs)
# --------- Step 3 ------------
#
# Computing orbital elements OTTIMO !!!!
#
mag_r = sqrt(dot(r,r))
mag_v = sqrt(dot(v,v))
#
#
print (' r_magnit= ', mag_r, ' km')
print (r, ' km')
print (' v_magnit= ', mag_v, ' km/s')
print (v, ' km/s')
a = 1 / ((2 / mag_r) - (mag_v**2 / mu))
print (' a = ', a, ' km')
#
h = cross(r, v)
mag_h = sqrt(dot(h,h))
#
e = cross(v, h) / mu - r / mag_r
mag_e = sqrt(dot(e,e))
print (' Eccentr.= ', mag_e, ' ')
#
i = acos(clip(h[2] / mag_h, -1, 1))
i = degrees(i)
if i >= 360.0:
i = i - 360
print (' Inclin.= ', i, ' degs')
#
true_anom = acos(clip(dot(e,r)/(mag_r * mag_e), -1, 1))
if dot(r, v) < 0:
true_anom = 2 * pi - true_anom
true_anom = degrees(true_anom)
print (' TrueAnom= ', true_anom, ' degs')
#
n = array([-h[1], h[0], 0])
mag_n = sqrt(dot(n,n))
raan = acos(clip(n[0] / mag_n, -1, 1))
if n[1] < 0:
raan = 2 * pi - raan
raan = degrees(raan) # right ascension
if raan >= 360.0:
raan = raan - 360
print (' RAAN = ', raan, ' degs')
#
per = acos(clip(dot(n, e) / (mag_n * mag_e), -1, 1))
if e[2] < 0:
per = 2 * pi - per
per = degrees(per)
if per >= 360.0:
per = per - 360
print (' ArgPerig= ', per, ' degs')
# EOF: eci_pv_ele.py ----------
# -------- OUTPUT ------------
'''
Prediction date: 2019 6 29
Prediction time: 11 3 57 UTC
=================
Name = STARLINK BH
Epoch= 2019-06-29 13:07:58.276991
Numb = 44290
r_magnit= 6933.665921356643 km
(-1355.213422491646, -4262.1234760085545, -5298.435916734078) km
v_magnit= 7.57644238427609 km/s
(7.367928372488619, -0.1610445730426039, -1.75788945773786) km/s
a = 6923.403686250501 km
Eccentr.= 0.0015107464813216766
Inclin.= 52.99124330893892 degs
TrueAnom= 168.87119780307017 degs
RAAN = 9.106093762412517 degs
ArgPerig= 84.25437668126654 degs
Name = STARLINK BJ
Epoch= 2019-06-28 22:40:57.025344
Numb = 44291
r_magnit= 6934.320213672662 km
(1630.9316727225473, -3921.3155040793654, -5481.618686312147) km
v_magnit= 7.5752836986374925 km/s
(7.298726436778112, 1.8394431010552108, 0.8543796544274769) km/s
a = 6922.597171559926 km
Eccentr.= 0.0016991442780485491
Inclin.= 52.98960924345745 degs
TrueAnom= 175.31411363712274 degs
RAAN = 9.236579780294168 degs
ArgPerig= 102.81355374111698 degs
Name = STARLINK BK
Epoch= 2019-06-29 17:45:11.129760
Numb = 44292
r_magnit= 6935.210651584592 km
(2899.7760915879658, -3518.157271489358, -5225.994145956705) km
v_magnit= 7.575246446991064 km/s
(6.811936832940996, 2.6405012692003664, 2.002405644224909) km/s
a = 6924.304370351646 km
Eccentr.= 0.0015752292006502545
Inclin.= 52.99260403068732 degs
TrueAnom= 180.80696215081082 degs
RAAN = 9.2649954945731 degs
ArgPerig= 108.52221619274525 degs
'''
Find below 3obj.txt database
STARLINK BH 1 44290U 19029BH 19180.54720228 .00000858 00000-0 78230-4 0 9994 2 44290 53.0070 8.7063 0002316 135.1721 224.9455 15.05479827 4858 STARLINK BJ 1 44291U 19029BJ 19179.94510446 .00000461 00000-0 50934-4 0 9997 2 44291 53.0078 11.5594 0002670 138.6184 221.5007 15.05488288 4770 STARLINK BK 1 44292U 19029BK 19180.73971215 .00000883 00000-0 80344-4 0 9998 2 44292 53.0074 8.0294 0002634 127.9210 232.2016 15.05287249 6737
Secondo esempio
# ---- azel_TLE_F.py ---------- Jun.23, 2019 ------ F stands for Final ----------
# Azimuth, elevation, distance froma TLE's elements
# Source: https://bwinkel.github.io/pycraf/satellite/index.html
import datetime
import numpy as np
from astropy.coordinates import EarthLocation
from astropy import time
from pycraf import satellite
with open ("short.txt", "r") as myfile:
data = myfile.read().splitlines()
for k in range(0,len(data),3):
tle= data[k]+'\n'+data[k+1]+'\n'+data[k+2]+'\n'
satname, sat = satellite.get_sat(tle)
print ('\n Sat Name = ', satname,'\n Sat Epoch= ', sat.epoch, '\n Sat Numb = ', sat.satnum)
# Datetime input
year=2019; month=6; day=16; utch=13; utcm=7; utcs=4
print ('\n Prediction date: ', year, month, day )
print (' Prediction time: ', utch, utcm, utcs, ' UTC' )
dt = datetime.datetime(year, month, day, utch, utcm, utcs)
obstime = time.Time(dt)
# define observer location
Lambda=15.0650; Phi=37.0328; H=370.0
location = EarthLocation(Lambda, Phi, H)
# create a SatelliteObserver instance
sat_obs = satellite.SatelliteObserver(location)
#
az, el, dist = sat_obs.azel_from_sat(tle, obstime)
print('\n Observer location: Lambda=', Lambda, 'deg Phi=', Phi, 'deg H=', H,' m')
print(' Azimuth : {:.5f}'.format(az))
print(' Elevation: {:.5f}'.format(el))
print(' Distance : {:.3f}'.format(dist))
print(' -------------------------')
# EOF azel_TLE_F.py ----------
# -------- OUTPUT ------------
'''
Sat Name = FALCON 9 DEB
Sat Epoch= 2019-06-16 20:44:25.611647
Sat Numb = 44298
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Observer location: Lambda= 15.065 deg Phi= 37.0328 deg H= 370.0 m
Azimuth : -102.93305 deg
Elevation: -6.12281 deg
Distance : 3167.504 km
-------------------------
Sat Name = STARLINK AB
Sat Epoch= 2019-06-06 23:49:07.999679
Sat Numb = 44260
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Observer location: Lambda= 15.065 deg Phi= 37.0328 deg H= 370.0 m
Azimuth : 72.53671 deg
Elevation: -72.28934 deg
Distance : 12507.391 km
-------------------------
Sat Name = STARLINK E
Sat Epoch= 2019-06-16 03:12:20.408255
Sat Numb = 44239
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Observer location: Lambda= 15.065 deg Phi= 37.0328 deg H= 370.0 m
Azimuth : 41.06985 deg
Elevation: 0.45210 deg
Distance : 2666.146 km
-------------------------
Sat Name = STARLINK P
Sat Epoch= 2019-06-16 20:40:12.943488
Sat Numb = 44248
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Observer location: Lambda= 15.065 deg Phi= 37.0328 deg H= 370.0 m
Azimuth : 47.41357 deg
Elevation: -8.69939 deg
Distance : 3830.295 km
-------------------------
Sat Name = STARLINK D
Sat Epoch= 2019-06-16 20:36:33.732864
Sat Numb = 44238
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Observer location: Lambda= 15.065 deg Phi= 37.0328 deg H= 370.0 m
Azimuth : 51.74955 deg
Elevation: -17.69062 deg
Distance : 5262.991 km
-------------------------
Sat Name = STARLINK H
Sat Epoch= 2019-06-16 17:21:59.422751
Sat Numb = 44242
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Observer location: Lambda= 15.065 deg Phi= 37.0328 deg H= 370.0 m
Azimuth : 54.88012 deg
Elevation: -27.18010 deg
Distance : 6902.716 km
-------------------------
'''
Find below short.txt database
FALCON 9 DEB 1 44298U 19029BR 19167.86418532 .00001313 00000-0 38063-4 0 9991 2 44298 53.0043 62.2948 0004912 13.4119 346.7015 15.43512957 4901 STARLINK AB 1 44260U 19029AB 19157.99245370 .01269439 00000-0 29305-1 0 9993 2 44260 53.0001 109.7412 0006821 81.1213 255.2409 15.39387046 2245 STARLINK E 1 44239U 19029E 19167.13356954 -.00013522 00000-0 -91175-3 0 9994 2 44239 53.0083 69.0416 0002637 113.5217 246.6049 15.05529924 3542 STARLINK P 1 44248U 19029P 19167.86126092 -.00617922 00000-0 -42907-1 0 9993 2 44248 52.9974 64.8224 0004816 42.9633 317.1731 15.07410453 3682 STARLINK D 1 44238U 19029D 19167.85872376 -.00444761 00000-0 -30673-1 0 9990 2 44238 52.9982 64.8538 0002160 82.0309 278.0925 15.07669053 2123 STARLINK H 1 44242U 19029H 19167.72360443 -.01130229 00000-0 -86918-1 0 9993 2 44242 52.9937 65.5606 0008561 56.4128 303.7676 15.04350216 3717
Terzo esempio
# ---- eci_angle.py ---------- May.24, 2019 ------
# Source: https://bwinkel.github.io/pycraf/satellite/index.html
from numpy import sqrt, arccos, pi
from astropy.coordinates import EarthLocation
from astropy import time
from pycraf import satellite
#
# ECI stands for Earth Center Inertial (frame)
# --------- Step 1 -------------
tle_string = '''ISS (ZARYA)
1 25544U 98067A 19141.19851322 .00000663 00000-0 18020-4 0 9997
2 25544 51.6411 135.2638 0001720 14.4836 102.9344 15.52691834171066'''
satname, sat = satellite.get_sat(tle_string)
print ('\n Sat Name = ', satname,'\n Sat Epoch= ', sat.epoch, '\n Sat Numb = ', sat.satnum)
# --------- Step 2 ------------
# Given prediction datetime
year=2019; month=5; day=24; utch=19; utcm=5; utcs=28
print ('\n Prediction date: ', year, month, day )
print (' Prediction time: ', utch, utcm, utcs, ' UTC' )
# using sgp4 directly, to get position and velocity in ECI coordinates
position, velocity = sat.propagate(year, month, day, utch, utcm, utcs)
print ('\n ECI-sat-pos= ', position, ' km')
posvector= sqrt(position[0]**2+position[1]**2+position[2]**2)
print(' ECI-Position vector magnitude: {:.3f}'.format(posvector), ' km')
print ('\n ECI-sat-vel= ', velocity, ' km/s')
velvector= sqrt(velocity[0]**2+velocity[1]**2+velocity[2]**2)
print(' ECI-Velocity vector magnitude: {:.3f}'.format(velvector), ' km/s')
# --------- Step 3 ------------
# Given elapsed time to compute ECI angle separation
elapt= 120 # seconds
print('\n Given elapsed time= ', elapt, 's')
pos2,vel2 = sat.propagate(year, month, day, utch, utcm, utcs+elapt)
print (' New position= ',pos2, ' km')
pos2vector= sqrt(pos2[0]**2+pos2[1]**2+pos2[2]**2)
print(' ECI-Pos2 vector magnitude: {:.3f}'.format(pos2vector), ' km')
# Compute angle between two vectors (ECI frame)
dotp= position[0]*pos2[0]+position[1]*pos2[1]+position[2]*pos2[2]
anglerad= arccos(dotp/(posvector*pos2vector))
print('\n ECI-angle btw. 2 pos-vectors= {:.4f}'.format (anglerad*180/pi), ' degs')
# EOF: eci_angle.py ----------
Sat Name = ISS (ZARYA) Sat Epoch= 2019-05-21 04:45:51.542207 Sat Numb = 25544 Prediction date: 2019 5 24 Prediction time: 19 5 28 UTC ECI-sat-pos= (-4357.6879564476385, -488.8784918567767, 5170.24402837909) km ECI-Position vector magnitude: 6779.371 km ECI-sat-vel= (2.452799491572253, -7.132278954255017, 1.3908781477772678) km/s ECI-Velocity vector magnitude: 7.669 km/s Given elapsed time= 120 s New position= (-4024.2879354447355, -1337.6492649089546, 5289.0839039880875) km ECI-Pos2 vector magnitude: 6779.278 km ECI-angle btw. 2 pos-vectors= 7.7781 degs
Quarto esempio
# ---- eci_long_F.py ---------- Jun.22, 2019 ------ F stands for Final ----
# Source: https://bwinkel.github.io/pycraf/satellite/index.html
from numpy import sqrt, arccos, pi, arctan2
from astropy.coordinates import EarthLocation
from astropy import time
from pycraf import satellite
#
# ECI stands for Earth Center Inertial (frame)
# --------- Step 1 -------------
with open ("short.txt", "r") as myfile:
data = myfile.read().splitlines()
for k in range(0,len(data),3):
tle= data[k]+'\n'+data[k+1]+'\n'+data[k+2]+'\n'
satname, sat = satellite.get_sat(tle)
print ('\n Sat Name = ', satname,'\n Sat Epoch= ', sat.epoch, '\n Sat Numb = ', sat.satnum)
# --------- Step 2 ------------
# Given prediction datetime
year=2019; month=6; day=16; utch=13; utcm=7; utcs=4
print ('\n Prediction date: ', year, month, day )
print (' Prediction time: ', utch, utcm, utcs, ' UTC' )
# using sgp4 directly, to get position and velocity in ECI coordinates
position, velocity = sat.propagate(year, month, day, utch, utcm, utcs)
#
Long_Sat= arctan2(position[1],position[0])
if Long_Sat < 0.0:
Long_Sat= Long_Sat + 2*pi
print ('\n Sat Position= ', position, ' km')
print (' Sat Longitude= {:.4f}'.format (Long_Sat*180/pi), ' degs')
# EOF: eci_long_F.py ----------
# -------- OUTPUT ----------
'''
Sat Name = FALCON 9 DEB
Sat Epoch= 2019-06-16 20:44:25.611647
Sat Numb = 44298
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Sat Position= (419.33398633374367, 6069.554240285104, 3054.3438753021123) km
Sat Longitude= 86.0478 degs
Sat Name = STARLINK AB
Sat Epoch= 2019-06-06 23:49:07.999679
Sat Numb = 44260
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Sat Position= (-988.5971434343431, -6181.851363990614, -2443.805353782833) km
Sat Longitude= 260.9142 degs
Sat Name = STARLINK E
Sat Epoch= 2019-06-16 03:12:20.408255
Sat Numb = 44239
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Sat Position= (-3302.402781352367, 2721.3833338073237, 5437.730905861728) km
Sat Longitude= 140.5094 degs
Sat Name = STARLINK P
Sat Epoch= 2019-06-16 20:40:12.943488
Sat Numb = 44248
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Sat Position= (-3870.559624372596, 1534.0515298425196, 5516.94297946272) km
Sat Longitude= 158.3797 degs
Sat Name = STARLINK D
Sat Epoch= 2019-06-16 20:36:33.732864
Sat Numb = 44238
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Sat Position= (-4394.828100428529, 1.4455234644697357, 5335.388637821106) km
Sat Longitude= 179.9812 degs
Sat Name = STARLINK H
Sat Epoch= 2019-06-16 17:21:59.422751
Sat Numb = 44242
Prediction date: 2019 6 16
Prediction time: 13 7 4 UTC
Sat Position= (-4702.101942994376, -1822.7655877255017, 4741.511105182952) km
Sat Longitude= 201.1888 degs
'''
Find below short.txt database
FALCON 9 DEB 1 44298U 19029BR 19167.86418532 .00001313 00000-0 38063-4 0 9991 2 44298 53.0043 62.2948 0004912 13.4119 346.7015 15.43512957 4901 STARLINK AB 1 44260U 19029AB 19157.99245370 .01269439 00000-0 29305-1 0 9993 2 44260 53.0001 109.7412 0006821 81.1213 255.2409 15.39387046 2245 STARLINK E 1 44239U 19029E 19167.13356954 -.00013522 00000-0 -91175-3 0 9994 2 44239 53.0083 69.0416 0002637 113.5217 246.6049 15.05529924 3542 STARLINK P 1 44248U 19029P 19167.86126092 -.00617922 00000-0 -42907-1 0 9993 2 44248 52.9974 64.8224 0004816 42.9633 317.1731 15.07410453 3682 STARLINK D 1 44238U 19029D 19167.85872376 -.00444761 00000-0 -30673-1 0 9990 2 44238 52.9982 64.8538 0002160 82.0309 278.0925 15.07669053 2123 STARLINK H 1 44242U 19029H 19167.72360443 -.01130229 00000-0 -86918-1 0 9993 2 44242 52.9937 65.5606 0008561 56.4128 303.7676 15.04350216 3717
Quinto esempio
''' kep2pv.py [Jul.08,2019]
Pos/Vel satellite vectors known 6 keplerian elements: (a, e, i, omega, Omega, M)
Source: https://github.com/aerospaceresearch/orbitdeterminator [file new_tle_kep_state.py]
'''
from numpy import sqrt, sin, cos, arccos, pi, arctan2, dot, cross, clip, array, zeros
from math import acos, degrees, atan2, radians
from scipy.optimize import fsolve
mu = 398600.4418
def TtoE(T,e):
E = atan2((1-e**2)**0.5*sin(T),e+cos(T))
E = E%(2*pi)
return E
def kep_to_state(kep):
""" This function converts from keplerian elements to the position and velocity vector
Args:
kep(1x6 numpy array): kep contains the following variables
kep[0] = semi-major axis (kms)
kep[1] = eccentricity (number)
kep[2] = inclination (degrees)
kep[3] = argument of perigee (degrees)
kep[4] = right ascension of ascending node (degrees)
kep[5] = true anomaly (degrees)
Returns:
r: 1x6 numpy array which contains the position and velocity vector
r[0],r[1],r[2] = position vector [rx,ry,rz] km
r[3],r[4],r[5] = velocity vector [vx,vy,vz] km/s
"""
r = zeros((6,1))
# unload orbital elements array
sma = kep[0] # semi major axis
ecc = kep[1] # eccentricity
inc = radians(kep[2]) # inclination
argp = radians(kep[3]) # argument of perigee
raan = radians(kep[4]) # right ascension of ascending node
eanom = TtoE(radians(kep[5]), ecc) # calling function TtoE
smb = sma * sqrt(1-ecc**2)
x = sma * (cos(eanom) - ecc)
y = smb * sin(eanom)
# calculate position and velocity in orbital frame
m_dot = (mu/sma**3)**0.5
e_dot = m_dot/(1 - ecc*cos(eanom))
x_dot = -sma * sin(eanom) * e_dot
y_dot = smb * cos(eanom) * e_dot
# rotate them by argp degrees
x_rot = x * cos(argp) - y * sin(argp)
y_rot = x * sin(argp) + y * cos(argp)
x_dot_rot = x_dot * cos(argp) - y_dot * sin(argp)
y_dot_rot = x_dot * sin(argp) + y_dot * cos(argp)
# convert them into 3D coordinates
r[0] = x_rot * cos(raan) - y_rot * sin(raan) * cos(inc)
r[1] = x_rot * sin(raan) + y_rot * cos(raan) * cos(inc)
r[2] = y_rot * sin(inc)
r[3] = x_dot_rot * cos(raan) - y_dot_rot * sin(raan) * cos(inc)
r[4] = x_dot_rot * sin(raan) + y_dot_rot * cos(raan) * cos(inc)
r[5] = y_dot_rot * sin(inc)
return r
# ------ Input array (manually): a, e, i, omega, Omega, M ----------
ele = array([6773.731405415944,0.0004993439967060411, 52.99724192977879, 13.163778174190766, 322.91113283192016, 39.37047972093648])
r = kep_to_state(ele)
#
print ('\n Given satellite orbital elements')
print(' a = 6773.731405415944 km')
print(' Eccentr.= 0.0004993439967060411')
print(' Inclin.= 52.99724192977879 degs')
print(' ArgPerig= 13.163778174190766 degs')
print(' RAAN = 322.91113283192016 degs')
print(' TrueAnom= 39.37047972093648 degs')
#
print ('\n Computed pos/vel vectors and related magnitudes')
print(' Pos-vector: (rx,ry,rz)= ', r[0],r[1],r[2], ' km')
print(' Vel-vector: (vx,vy,z)= ', r[3],r[4],r[5], ' km/s')
#
pos= sqrt(r[0]**2+r[1]**2+r[2]**2)
vel= sqrt(r[3]**2+r[4]**2+r[5]**2)
print('\n Pos.Magnitude= ',pos, ' km')
print(' Vel.Magnitude= ',vel, ' km/s')
print('\n')
'''
Given satellite orbital elements
a = 6773.731405415944 km
Eccentr.= 0.0004993439967060411
Inclin.= 52.99724192977879 degs
ArgPerig= 13.163778174190766 degs
RAAN = 322.91113283192016 degs
TrueAnom= 39.37047972093648 degs
Computed pos/vel vectors and related magnitudes
Pos-vector: (rx,ry,rz)= [5236.17488962] [96.3779873] [4291.99189931] km
Vel-vector: (vx,vy,z)= [-3.16266643] [5.91433255] [3.7294375] km/s
Pos.Magnitude= [6771.11590913] km
Vel.Magnitude= [7.67401411] km/s
'''
Sesto esempio
# ------ rapid.py ---------
from math import pi
nsat = float(input('\n Satellite rounds per day (source TLE), try 15.56790542: '))
#
period= 1440.0/nsat
No= nsat*pi/720
a= (No*13.44685206)**(-2.0/3.0)
akm= a*6378.137
#
print ('\n Period T= {:.3f}'.format(period),' minutes \n Semi-major axis a= {:.3f}'.format(akm), ' km')
# EOF: rapid.py
'''
Satellite rounds per day (source TLE), try 15.56790542: 15.56790542
Period T= 92.498 minutes
Semi-major axis a= 6775.090 km
'''
Cosa fa questo script Python in più
rispetto a un normale planetario?
(Maggio 2019). Semplicemente mette in grafico la posizione
di una stella, tra le migliaia disponibili sul database, confrontandola
con azimut e altezza di Sole e Luna e valutando la zona del crepuscolo
per individuare l'orario di migliore visibilità.
L'ascissa del grafico mostra le 12 ore antecedenti la mezzanotte e
quelle successive, mentre il relativo codice Python (.py) è riportato
di sotto.
Per settare una nuova stella (nome oppure #SAO), basta modificare le seguenti righe:
11 body= ... inserire il nuovo nome 13 cambio coordinate geografiche, quando serve 14 cambio ora solare/legale [1 oppure 2] 15 data di riferimento, sempre a mezzanotte 34,45 trasferire nome e data nelle etichette
# ------ final.py -------- May.10-11, 2019 ------
# Necessita apposito settaggio in ambiente Linux
# Warning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
# % get_backend())
import numpy as np
import matplotlib.pyplot as plt
from astropy.visualization import astropy_mpl_style
import astropy.units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy.coordinates import get_sun
from astropy.coordinates import get_moon
#
body = SkyCoord.from_name('deneb')
#
loc = EarthLocation(lat=37.0328*u.deg, lon=15.0650*u.deg, height=370*u.m)
utcoffset = 2*u.hour # TMEC Daylight Time
midnight = Time('2019-5-11 00:00:00') - utcoffset
#
# ---------- Get sun, moon and body ---------
# sun
delta_midnight = np.linspace(-12, 12, 1000)*u.hour
times = midnight + delta_midnight
frame = AltAz(obstime=times, location=loc)
sunaltazs = get_sun(times).transform_to(frame)
# moon
moon = get_moon(times) # calling moon (+10Mb), for more precision
moonaltazs = moon.transform_to(frame)
# body
bodyaltazs = body.transform_to(frame)
#
# ----------------- Plotting Altaz Graphs [Sun,Moon,body] --------------
plt.style.use(astropy_mpl_style)
plt.plot(delta_midnight, sunaltazs.alt, color='r', label='Sun')
plt.plot(delta_midnight, moonaltazs.alt, color=[0.75]*3, ls='--', label='Moon')
plt.scatter(delta_midnight, bodyaltazs.alt,
c=bodyaltazs.az, label='body=deneb', lw=0, s=8,
cmap='viridis_r')
plt.fill_between(delta_midnight.to('hr').value, 0, 90,
sunaltazs.alt < -0*u.deg, color='0.5', zorder=0)
plt.fill_between(delta_midnight.to('hr').value, 0, 90,
sunaltazs.alt < -18*u.deg, color='k', zorder=0)
plt.colorbar().set_label('Azimuth [deg]')
plt.legend(loc='upper left')
plt.xlim(-12, 12)
plt.xticks(np.arange(13)*2 -12)
plt.ylim(0, 90)
plt.xlabel('Hours from TMEC Midnight [2019-5-11]')
plt.ylabel('Altitude [deg]')
plt.show()
# EOF: final.py ----------
May.11, 2019
La potente libreria AstroPy di Python per il calcolo di
effemeridi
di precisione tramite il dbase DE430
(Maggio 2019). La versatilità del linguaggio Python, che è multi-piattaforma,
interpretato e direttamente avviabile da console, permette all'utente di ridurre drasticamente la
lunghezza dei programmi informatici. Su due fronti, input delle date delle effemeridi, e,
grazie alla libreria astropy,
utilizzo diretto dell'enorme database DE430 (di circa 120Mb) senza bisogno di scaricarlo nel computer, ma
semplicemente 'agganciandolo' da comando come nell'esempio seguente.
Nel suo output ho anche aggiunto un confronto con il risultato ottenuto con
l'eccellente codice Planet dell'astrofilo e amico Aldo Nicola,
compilato da Fortran90 tramite le insuperabili librerie SOFA.
# -------- siteupload.py --------- May.04, 2019 -----------
# Installazione libreria: pip install astropy
# Avvio da sistemi Linux: sudo python3 siteupload.py
# Avvio da sistemi Windows: python siteupload.py
from astropy.time import Time
import astropy.units as u
from astropy.coordinates import solar_system_ephemeris, EarthLocation
from astropy.coordinates import get_body_barycentric, get_body, get_body_barycentric_posvel
#
t = Time('2017-08-21 17:28:11.284', format='iso', scale='tt', precision= 6)
#
with solar_system_ephemeris.set('de430'):
jup = get_body('jupiter', t)
print('\n Jupiter Geocentric [RA,Dec,Dist]= \n', jup)
#
casamia = EarthLocation(lat=37.0328*u.deg,lon=15.0650*u.deg,height=370*u.m)
jup2 = get_body('jupiter', t, casamia)
print('\n Jupiter Topocentric [RA,Dec,Dist]= \n', jup2)
#
jup3= get_body_barycentric_posvel('jupiter', t)
print('\n Jupiter Geocentric Position [x,y,z] and Velocity [Vx,Vy,Vz]= ')
print(jup3)
#
print ('\n',solar_system_ephemeris.bodies) # nomenclatura pianeti
# EOF: siteupload.py -------
Output
Jupiter Geocentric [RA,Dec,Dist]=
<SkyCoord (GCRS: obstime=2017-08-21 17:28:11.284000, obsgeoloc=(0., 0., 0.) m,
obsgeovel=(0., 0., 0.) m / s):
(ra, dec, distance) in (deg, deg, km)
(198.83064461, -6.75384426, 9.00946238e+08)>
Jupiter Topocentric [RA,Dec,Dist]=
<SkyCoord (GCRS: obstime=2017-08-21 17:28:11.284000,
obsgeoloc=(-1999901.58608911, -4686828.13847154, 3823744.31696533) m,
obsgeovel=(341.7577692, -146.30786813, -0.58520951) m / s):
(ra, dec, distance) in (deg, deg, AU)
(198.82993429, -6.75479066, 6.0222937)>
Jupiter Geocentric Position [x,y,z] and Velocity [Vx,Vy,Vz]=
(<CartesianRepresentation (x, y, z) in AU
(-4.79425873, -2.40930992, -0.91613978)>,
<CartesianRepresentation (x, y, z) in AU / d
(0.00348154, -0.00575187, -0.00255017)>)
('earth', 'sun', 'moon', 'mercury', 'venus', 'earth-moon-barycenter',
'mars', 'jupiter', 'saturn', 'uranus', 'neptune')
Dal software Planet, by Aldo Nicola, astrofilo torinese (Chieri)
Barycentric equatorial position and velocity
mean equator and equinox of J2000.
X Y Z
POS -4.7943811026733227 -2.4094085431640209 -0.9160860398906552 AU
VEL 0.0034822044780244 -0.0057500941888184 -0.0025493991894452 AU/day
-------------
Breve nota: i lievissimi scarti (pos/vel) sono dovuti ai due differenti equinozi:
della data il primo e J2000 il secondo
-------------
(05 Mag.2019)
La Grande Opposizione di Marte:
27-Luglio-2018, 05h 07m 11.5s TDT
(Luglio 2018). Come nel 2003, questo evento può essere
monitorato utilizzando lo stesso software di allora scaricabile da
questo file gzippato. Si ottiene la
tabella seguente, con l'istante esatto, calcolato fino a far coincidere 6 decimali di grado
in longitudine eliocentrica di Marte e Terra, e la distanza reciproca pari a 0.386158 UA = 57,768,424 km.
(Aggiornamento 29-07-18). Inserito in forma zippata il codice fortran del
programma: eccolo.
Il compilatore usato nell'esempio è gfortran della distribuzione Ubuntu-16.04.
Istante della Grande Opposizione di Marte
------------------------------------------
27 / 07 / 2018
05h 07m 11.5s TDT
che si deduce dagli scarti delle longitudini eliocentriche in basso
------------------------------------------------------------------------
1..
Data : 27 / 07 / 2018 Giorno Giuliano : 2458326.713322
Ora : 05h 07m 11s TDT
D - M A R T E R
Distanza Marte-Terra calcolata con VSOP87 del BdL
*** Usati in totale 12497 termini delle serie ***
Coordinate Eclittiche Eliocentriche (Equinozio della Data)
Xmarte : 0.7854888857 ua Xterra : 0.5700921835 ua
Ymarte : -1.1579521698 ua Yterra : -0.8404186816 ua
Zmarte : -0.0435294520 ua Zterra : -0.0000001471 ua
rm : 1.3999074363 ua rt : 1.0155336814 ua
Long-ma: 304.150737 gr Long-te: 304.150735 gr [Scarto long.elioc. +2E-6 gradi]
------------------------------------------------------------------------
Distanza: 0.3861580710 ua = 57768425. km
Distanza angolare: 24.330 arcsec
------------------------------------------------------------------------
Scelta del tempo : avanti -> + indietro -> - stop -> 0
Scelta ? +
2..
Data : 27 / 07 / 2018 Giorno Giuliano : 2458326.713333
Ora : 05h 07m 12s TDT
D - M A R T E R
Distanza Marte-Terra calcolata con VSOP87 del BdL
*** Usati in totale 12497 termini delle serie ***
Coordinate Eclittiche Eliocentriche (Equinozio della Data)
Xmarte : 0.7854890258 ua Xterra : 0.5700923450 ua
Ymarte : -1.1579520649 ua Yterra : -0.8404185706 ua
Zmarte : -0.0435294532 ua Zterra : -0.0000001471 ua
rm : 1.3999074282 ua rt : 1.0155336802 ua
Long-ma: 304.150744 gr Long-te: 304.150746 gr [Scarto long.elioc. -2E-6 gradi]
------------------------------------------------------------------------
Distanza: 0.3861580643 ua = 57768424. km
Distanza angolare: 24.330 arcsec
------------------------------------------------------------------------
La distanza minima, invece, si avrà qualche giorno dopo, come mostrato in questa tabella:
Data : 31 / 07 / 2018 Giorno Giuliano : 2458330.826389
Ora : 07h 50m 00s TDT
Distanza: 0.3849629045 ua = 57589631. km
Distanza angolare: 24.405 arcsec
MatLab/Octave: effemeridi planetarie di grande precisione
(28 Gennaio 2018). Le potenti funzioni che il prof. Mahooti ha pubblicato
nel repository MathWorks di figura sono state processate opportunamente per
tabulare i parametri più immediati, ovvero le coordinate rettangolari e polari,
sia eliocentriche che geocentriche dei nove pianeti del sistema solare, oltre
a Luna e Sole apparente. Il main file che le richiama (000_main.m)
è inserito nel pacchetto contenente tutte le funzioni e può
essere prelevato da qui, mentre i
risultati di un esempio applicativo sono tabulati sotto.
Package including all functions and the main listing (000_main.m)
are available here.
Ephemeris Date: 2017,8,21
Time: 17,28,11.284 TT
------- Heliocentric equatorial Coordinates --------
in Astronomical Units
------------------------------
1) |r|= 0.41736838 AU (Mercury)
x= 0.29057210
y= -0.25070626
z= -0.16404447
------------------------------
2) |r|= 0.72152571 AU (Venus)
x= 0.33233464
y= 0.59164277
z= 0.24517723
------------------------------
4) |r|= 1.65508612 AU (Mars)
x= -1.17862515
y= 1.04377710
z= 0.51057044
------------------------------
5) |r|= 5.44819162 AU (Jupiter)
x= -4.79690347
y= -2.41444970
z= -0.91811855
------------------------------
6) |r|= 10.06113755 AU (Saturn)
x= -0.65142839
y= -9.28942712
z= -3.80902529
------------------------------
7) |r|= 19.91679327 AU (Uranus)
x= 17.95595033
y= 7.98410750
z= 3.24291987
------------------------------
8) |r|= 29.94682316 AU (Neptune)
x= 28.56036259
y= -8.06374648
z= -4.01172018
------------------------------
9) |r|= 33.39128916 AU (Pluto)
x= 10.37234786
y= -29.27504096
z= -12.26232315
------------------------------
------- GEOCENTRIC Equatorial and Polar Coordinates --------
in Astronomical Units and Equinox J2000
1G) |r|= 0.62003111 AU (Mercury_GEOC)
x= -0.57285543
y= 0.23280702
z= 0.04556457
RA= 157.88330 deg | long_ecl= 157.97660 deg
DEC= 4.21433 deg | lat_ecl= -4.69965 deg
-----------------------------------------------------
2G) |r|= 1.28251734 AU (Venus_GEOC)
x= -0.53109289
y= 1.07515605
z= 0.45478627
RA= 116.28791 deg | long_ecl= 114.46363 deg
DEC= 20.76920 deg | lat_ecl= -0.46526 deg
-----------------------------------------------------
4G) |r|= 2.64976481 AU (Mars_GEOC)
x= -2.04205268
y= 1.52729038
z= 0.72017948
RA= 143.20648 deg | long_ecl= 140.42667 deg
DEC= 15.77081 deg | lat_ecl= 1.15108 deg
-----------------------------------------------------
5G) |r|= 6.02244537 AU (Jupiter_GEOC)
x= -5.66033100
y= -1.93093642
z= -0.70850950
RA= 198.83630 deg | long_ecl= 199.93952 deg
DEC= -6.75620 deg | lat_ecl= 1.12305 deg
-----------------------------------------------------
6G) |r|= 9.63300080 AU (Saturn_GEOC)
x= -1.51485592
y= -8.80591384
z= -3.59941624
RA= 260.23911 deg | long_ecl= 260.95031 deg
DEC= -21.94119 deg | lat_ecl= 1.19199 deg
-----------------------------------------------------
7G) |r|= 19.38491406 AU (Uranus_GEOC)
x= 17.09252280
y= 8.46762078
z= 3.45252891
RA= 26.35377 deg | long_ecl= 28.14086 deg
DEC= 10.25934 deg | lat_ecl= -0.59290 deg
-----------------------------------------------------
8G) |r|= 28.96612152 AU (Neptune_GEOC)
x= 27.69693505
y= -7.58023321
z= -3.80211114
RA= 344.69382 deg | long_ecl= 343.00129 deg
DEC= -7.54244 deg | lat_ecl= -0.93590 deg
-----------------------------------------------------
9G) |r|= 32.62881471 AU (Pluto_GEOC)
x= 9.50892033
y= -28.79152768
z= -12.05271410
RA= 288.27674 deg | long_ecl= 286.94469 deg
DEC= -21.67792 deg | lat_ecl= 0.69269 deg
-----------------------------------------------------
------- Geocentric equatorial Coordinates --------
of Moon in kilometers and Sun apparent in Astronomical Units
------------------------------
a) |d|= 371950.219 km (Moon)
x= -315540.767
y= 179700.981
z= 80551.521
long_ecl= 150.33842 deg
lat_ecl= 12.50738 deg
---------------------------------------
b) |d|= 1.01154740 AU (Sun Apparent)
x= -0.86342753
y= 0.48351327
z= 0.20960904
long_ecl= 150.75148 deg
lat_ecl= 11.95927 deg
---------------------------------------
(13-16 Dicembre 2016). Operazione di restyling del codice asteroidale aster16.c pubblicato sotto. Per la lettura delle stringhe del database non viene usato il laborioso comando fseek ma fscanf che si avvale di un file di deposito esterno per copiarvi (in una riga) i dati del record intercettato. Ciò permette di modificare rapidamente il codice C in caso di variazione del template dei campi (larghezza e numero). L'algoritmo è stato provato sul file AsterSelect.c ed applicato successivamente allo script finale AsterDIC16.c, entrambi reperibili qui, insieme allo stesso dbase di 20,000 pianetini.
(28 Settembre 2016). Una interessante variante del software asteroidale è questa. I due Database degli elementi orbitali, richiamati dai codici AstGeocP.c e AstNEAp.c, contengono rispettivamente 20,000 asteroidi del catalogo MPC e 14,690 NEA, acronimo di Near Earth Asteroids; essi vengono scandagliati scegliendo come parametro di screening la distanza geocentrica; per una migliore consultazione, i risultati ottenuti sono copiati su file testo. Nel caso dei NEA, si percepisce subito quali siano i più vicini alla Terra alla data dell'effemeride, come nell'esempio seguente:
Distanza Geocentrica limite assunta= 0.02500 Data Effemeride= 2016 9 24 2016 RM20 Distanza GEOCENTRICA d= 0.01991001 (UA) Distanza Eliocentr. r= 1.01175656 (UA) Long.Eclitt.Eliocen. l= 1.07918 (gradi) Lati.Eclitt.Eliocen. b= -1.01969 (gradi) Ascensione Retta AR= 36.12525 (gradi) = 2H 24M 30S Declinazione DE= -58.36825 (gradi) psi= -114.98831 (gradi) Elongaz. Ovest Angolo di Fase beta= 63.98968 (gradi) Magnitudine= 20.0 2016 RD34 Distanza GEOCENTRICA d= 0.01327780 (UA) Distanza Eliocentr. r= 1.00928176 (UA) Long.Eclitt.Eliocen. l= 1.82229 (gradi) Lati.Eclitt.Eliocen. b= -0.01367 (gradi) Ascensione Retta AR= 62.36386 (gradi) = 4H 9M 27S Declinazione DE= 19.96857 (gradi) psi= -116.99746 (gradi) Elongaz. Ovest Angolo di Fase beta= 62.33090 (gradi) Magnitudine= 20.1 ============================= N.ro corpi selezionati= 2 =============================
Altro esempio, per gentile concessione dell'autore di AstGeoc.f90 [Fortran 90],
l'astrofilo piemontese Aldo Nicola.
Date and Time UTC
==============================
Input : Year , Month, Day (yyyy,mm,dd)
2016,10,1
Input : Maximum Geocentric distance 1.00 [UA]
(2016, revised by G. Matarazzo)
Time of Ephemeris ID 2457662.5007891473 Time of epoch JD 2457600.5000000000
--------------------------------------------------------------------------------------------------
R.A.2000 Decl.2000 Delta r Phase Mag Elong.
num body hh mm ss.ss dd pp ss.s A.U. A.U. deg V deg
--------------------------------------------------------------------------------------------------
00433 Eros 21 0 12.87 - 2 21 50.50 0.846252 1.657753 28.7 12.75 127.4 E
01863 Antinous 0 4 18.17 + 16 11 57.34 0.919574 1.905601 7.5 17.30 165.6 E
01865 Cerberus 5 55 50.21 + 19 54 23.11 0.685827 1.298868 49.6 18.44 99.0 W
01943 Anteros 1 37 11.49 + 25 41 24.22 0.783340 1.733142 15.6 17.27 152.2 W
02100 Ra-Shalom 3 7 42.15 - 25 51 57.79 0.162216 1.116143 41.8 14.03 131.9 W
03103 Eger 8 12 41.29 - 19 25 44.57 0.520085 0.909234 84.2 16.70 64.6 W
03200 Phaethon 13 58 8.90 + 82 51 22.89 0.403303 1.056688 71.0 15.24 86.6 E
03352 McAuliffe 3 1 54.05 + 6 45 50.91 0.818840 1.721235 21.0 17.58 141.9 W
03753 Cruithne 6 50 56.68 - 13 46 43.52 0.644850 1.133080 61.5 17.11 84.0 W
04179 Toutatis 18 52 11.08 - 23 21 37.69 0.633991 1.221312 54.9 16.86 93.9 E
04341 Poseidon 23 50 16.60 - 31 14 20.35 0.674437 1.598263 21.4 17.21 144.3 E
04581 Asclepius 17 8 41.34 - 19 0 22.81 0.727503 1.015213 67.8 22.43 69.9 E
04953 1990 MU 6 5 16.61 - 3 36 31.94 0.773186 1.325591 48.7 15.98 95.8 W
05381 Sekhmet 9 35 54.04 - 51 52 52.45 0.649299 0.885928 79.8 18.19 60.6 W
05653 Camarillo 4 18 22.87 + 27 1 12.91 0.934055 1.679978 30.9 18.40 120.4 W
05786 Talos 23 16 17.22 - 18 4 25.08 0.956216 1.899822 14.3 19.20 152.1 E
05836 1993 MF 3 10 7.85 + 30 24 38.95 0.545095 1.431187 30.6 15.58 133.3 W
05863 Tara 18 7 28.72 - 12 24 16.33 0.785987 1.206008 55.6 17.71 84.0 E
06239 Minos 16 7 49.78 - 23 12 17.01 0.946415 0.921460 64.8 20.50 56.4 E
06569 Ondaatje 3 18 54.45 - 34 37 49.91 0.794182 1.596334 30.9 18.24 125.1 W
07088 Ishtar 0 51 33.27 - 16 12 50.45 0.708463 1.684057 11.8 17.80 159.9 W
07341 1991 VK 0 18 39.83 + 19 14 29.00 0.497297 1.485093 10.9 16.72 163.7 E
07888 1993 UC 6 29 48.04 - 12 57 4.86 0.814011 1.281518 51.4 17.09 89.2 W
10150 1994 PN 16 20 16.84 - 56 9 18.88 0.904571 1.095014 59.2 17.50 69.9 E
10636 1998 QK56 23 10 39.21 - 10 15 51.21 0.539845 1.510947 15.5 18.01 156.2 E
11500 Tomaiyowit 5 21 31.17 + 33 9 0.18 0.516005 1.245886 50.6 19.22 105.9 W
15745 Yuliya 3 6 15.24 - 11 39 45.48 0.878535 1.757120 22.3 19.31 138.3 W
16834 1997 WU22 19 57 10.94 + 14 32 16.72 0.682906 1.407616 41.3 17.13 112.0 E
END OF PROGRAM
(Settembre 2016).Con il software asteroidale pubblicato qui appresso,
voglio rendere omaggio all'amico, collega ed astrofilo Sergio Foglia
ben conosciuto, in Italia e all'estero, per la sua poliedrica attività
di scopritore di pianetini e cultore di programmi astronomici di rara
eleganza e precisione.
Nel lontano 1996 ho avuto il privilegio di ricevere una personal release
del codice MPEPH.c per il calcolo delle effemeridi degli asteroidi tramite
lettura del Dbase degli elementi orbitali pubblicato dal Minor Planet Center (MPC).
Per avere un'idea più dettagliata basta leggere questo suo esaustivo
paper sull'argomento.
Va da sé, però, che nel corso degli anni il MPC ha cambiato lo standard del
database, aumentando considerevolmente il numero dei campi (fissi) da cui
estrarre i dati e le modalità di lettura di alcuni: stringhe invece di numeri.
Ho quindi rivisto la parte del codice da cambiare e con il permesso dell'autore
procedo alla pubblicazione della sua opera intellettuale, scaricabile da
qui; il file zippato comprende
il più recente dbase limitato ai primi 20mila asteroidi del catalogo MPC
e il nuovo listato aster16.c, per la cui compilazione è sufficiente
usare (gcc) di licenza GNU, gratuito e di facile reperibilità sia per sistemi
Windows che Linux.
Grazie di cuore, Sergio.
+-------------------------------------------+
| MINOR PLANET EPHEMERIS |
| M P E P H v. 2.1 |
| 1996, S.Foglia, Serafino Zani Observatory |
| (2016, revisited by G. Matarazzo) |
+-------------------------------------------+
Minor Planet catalogue number [1 to 20,000]: 1620
Anomalia Media M= 322.32014 (gradi)
Argomento Peri w= 276.86945 (gradi)
Longitud. Nodo O= 337.21625 (gradi)
Inclinazione i= 13.33747 (gradi)
Eccentricita' e= 0.3353926
Moto Medio n= 0.70916069 (gradi/g)
Semiasse Magg. a= 1.2453926 (UA)
Numero Catalogo : (1620)
Nome Pianetino : Geographos
Parametro H= 15.60
Parametro G= 0.15
Sigla Secolo : K [K=20; J=19]
Anni-Secolo = 2000
Sigla Anni : 16 [01-99]
Anni = 16
Sigla Mese : 7 [1-9; A=10; B=11; C=12]
Mese = 7
Sigla Giorni : V [1-9; A=10; B=11; ... U=30; V=31]
Giorno = 31
------------------------------------
Epoca Anomalia Media: anno= 2016
mese= 7
gior= 31
------------------------------------
Giorno Giuliano Epoca Anomalia Media= 2457600.50000
Ephemeris begin
year month day : 2016 09 14
Ephemeris end
year month day : 2016 09 30
Ephemeris step (days) : 2
Numero Catalogo: (1620)
Nome Pianetino: Geographos
MPEPH: Minor Planet Ephemeris
1996, S.Foglia, Serafino Zani Observatory
(2016, revisited by G. Matarazzo)
------------------------------------------------------------------------------
Date R.A.2000 Decl.2000 Delta r Phase Mag Elong.
year mo day hh mm ss.ss dd pp ss.s A.U. A.U. deg V deg
------------------------------------------------------------------------------
2016 9 14 13 10 48.65 -15 25 8.9 1.50147 0.83246 39.3 17.6 31.6E
2016 9 16 13 19 47.21 -16 22 43.0 1.48924 0.83041 39.9 17.6 32.0E
2016 9 18 13 28 56.08 -17 19 8.1 1.47726 0.82893 40.6 17.6 32.5E
2016 9 20 13 38 15.45 -18 14 13.7 1.46556 0.82803 41.2 17.6 32.9E
2016 9 22 13 47 45.46 -19 7 48.8 1.45417 0.82770 41.9 17.6 33.4E
2016 9 24 13 57 26.19 -19 59 42.4 1.44313 0.82795 42.5 17.6 33.9E
2016 9 26 14 7 17.66 -20 49 43.3 1.43247 0.82878 43.1 17.6 34.4E
2016 9 28 14 17 19.79 -21 37 39.9 1.42223 0.83018 43.7 17.6 34.9E
2016 9 30 14 27 32.44 -22 23 21.0 1.41244 0.83216 44.2 17.7 35.4E
------------------------------------------------------------------------------
(Agosto 2016). Rivisitazione del Teorema di Lambert con l'elegante algoritmo di Battin. Codice in linguaggio C e risultati sono riportati sotto.
/* Battin.c :: Teorema di Lambert (metodo di Battin) in C :: 13/15-Ago-2016
* convertito da BestBatt.bas */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.141592653589793238462643 /* in double conserva 15 cifre significative */
#define RAD 0.01745329251994329576923691 /* in double conserva 15 cifre significative */
#define KGAUSS 0.01720209895 /* in double mantiene le 11 cifre designate */
#define MU 1 /* parametro gravitazionale unitario espresso in UA^3/UT0^2 */
#define uv0 29.78469169 /* costante di conversione velocità da UC in km/s */
int main(void)
{
float r1,r2,tetagrad,tgiorni;
// double r1,r2,tetagrad,tgiorni;
double teta,t,c,s,lam,Tadim, cosf,sinf,tgf2,L,rop,m;
double x,y, yit, eta,gs0,gs1,gs2,gs3,gs4, csi;
double hden,h1,h2,U,B, fs0,fs1,fs2,fs3,fs4, ku;
double a,p,ecc, V1,V2, f,g,fpunto,gpunto,sigma1,V1L,V2L;
/*
+-------------------------------------------------------+
. Da BATTIN: An elegant Lambert algorithm - 1984 .
. (+Vaughan) (3^ parte dell'articolo- p.669-670) .
. Rimozione singolarità θ=180° e RAPIDA CONVERGENZA .
. con uso della frazione continua K(u) nonchè dell'ar- .
. ticio del calcolo di csi(x) tramite 'eta' .
+-------------------------------------------------------+
*/
printf("\tProblema di LAMBERT: calcolo della traiettoria di trasferimento\n");
printf("\tdi una sonda noti a) la posizione geometrica -> r1, r2, teta\n");
printf("\t b) il tempo di percorrenza -> t\n");
printf("\t\t\t [Moto Centrale Eliocentrico]\n");
printf("\n Raggio vettore r1 [UA], p.es. 1.0= "); scanf("%f", &r1);
printf(" Raggio vettore r2 [UA], p.es. 1.8= "); scanf("%f", &r2);
printf(" Angolo tra i due vettori teta [gradi], p.es. 109.884= "); scanf("%f", &tetagrad);
printf(" Tempo di percorrenza t [giorni], p.es. 253.86= "); scanf("%f", &tgiorni);
//
printf("\n\t r1= %.3f UA\t r2= %.3f UA\t teta= %.3f gradi\t delta-t= %.3f giorni\n\n", r1,r2,tetagrad,tgiorni);
/* ----------- Inizio calcoli ---------- */
teta= tetagrad*RAD; // angolo tra i due vettori in radianti
t= tgiorni*KGAUSS; // tempo in UC (Unità canoniche)
c= sqrt(r1*r1 + r2*r2 -2*r1*r2*cos(teta)); // c= corda del triangolo dei vettori dati
s= (r1 + r2 + c)/2 ; // semiperimetro del triangolo
//
lam= sqrt(r1*r2)/s * cos(teta/2) ; // 1^ parametro adimensionale
Tadim= sqrt(8*MU/s) * t/s; // 2^ parametro adimensionale
//
cosf= sqrt(r1*r2)*cos(teta/2)/((r1+r2)/2); sinf= sqrt(1-cosf*cosf);
tgf2= (1 - cosf)/sinf; L= tgf2*tgf2;
rop= (r1 + r2 + 2*lam*s)/4;
m= MU*t*t/(8*rop*rop*rop);
printf("\tParametri intermedi\n");
printf("\tc= %.14f UA s= %.14f UA\n\tlam=%.14f adim. Tadim=%1.14f adim. \n\n", c,s,lam,Tadim);
/* Da qui inizia il ciclo iterativo*/
x=L; // valore iniziale di x
int count=0;
goto label11;
for (;;) // ciclo infinito che è interrotto dal break
{
x= sqrt(((1 - L)/2)*((1 - L)/2) + m /(y*y)) - (1 + L)/2; /* 1^ equaz. di GAUSS */
label11:
{count++; // contatore iterazioni
eta= x/(1+sqrt(1+x))/(1+sqrt(1+x)); /* Param. ausiliario {eta} -> -1 < eta < +1 formula (52)*/
/* Calcolo della Funzione Ipergeometrica Csi = f(x,eta) */
gs3= 169.0/25.0/27.0*eta/(1 + 196.0/27.0/29.0*eta/(1 + 225.0/29.0/31.0*eta/(1 + gs4)));
gs2= 100.0/19.0/21.0*eta/(1 + 121.0/21.0/23.0*eta/(1 + 144.0/23.0/25.0*eta/(1 + gs3)));
gs1= 49.0/13.0/15.0*eta/(1 + 64.0/15.0/17.0* eta/(1 + 81.0/17.0/19.0*eta/(1 + gs2)));
gs0= 16.0/7.0/9.0*eta/(1 + 25.0/9.0/11.0* eta/(1 + 36.0/11.0/13.0*eta/(1 + gs1)));
csi= 8*(1 + sqrt(1 + x))/(3 + 1/(5 + eta + 9.0/7.0*eta/(1 + gs0))); // Formula (53)
/* Parametri ausiliari h1,h2 */
hden= (1 + 2*x + L)*(4*x + csi*(3 + x));
h1 = (L + x)*(L + x)*(1 + 3*x + csi)/hden;
h2 = m * (x - L + csi)/hden;
/* Parametri ausiliari B,U */
B= 27*h2/(4*(1+h1)*(1+h1)*(1+h1));
U= -B/2/(1 + sqrt(1 + B));
/*Calcolo della Funzione Ipergeometrica K(u) */
fs3= 928.0/3591.0*U/(1 - 1054.0/4347.0*U/(1 - 266.0/1035.0*U/(1 - fs4)));
fs2= 418.0/1755.0*U/(1 - 598.0/2295.0*U/(1 - 700.0/2907.0*U/(1 - fs3)));
fs1= 22.0/81.0*U/(1 - 208.0/891.0*U/(1 - 340.0/1287.0*U/(1 - fs2)));
fs0= 4.0/27.0*U/(1 - 8.0/27.0*U/(1 - 2.0/9.0*U/ (1 - fs1)));
ku= 1.0/3.0/(1 - fs0); //Formula (58)
/* 2^ Equaz. di Gauss -> Radice positiva y */
yit= (1 + h1)/3*(2 + sqrt(1 + B)/(1 - 2*U*ku*ku));
printf("\t %d) \t yit = %.8f\t Residuo= %.2e \n", count,yit,yit-y);
} // chiude da count
if (fabs(yit-y) < 1E-6)
break;
else
y=yit;
} // fine ciclo infinito di for
printf("\n\t\t y-fin = %.8f\t x-fin= %.8f \n", yit,x);
/* Calcolo dei PARAMETRI ORBITALI {a,p,e} */
a= m*s*(1 + lam)*(1 + lam)/(8 * x * yit*yit);
p= 2*r1*r2*(y*(1 + x)*sin(teta/2))*(y*(1 + x)*sin(teta/2))/((1 + lam)*(1 + lam))/m/s;
ecc= sqrt(1 - p/a);
/* Calcolo delle VELOCITA' {V1, V2} - Metodo CLASSICO */
V1= sqrt(2* MU/r1 - MU/a); V2= sqrt(2* MU/r2 - MU/a);
/* Calcolo delle VELOCITA' {V1, V2} tramite i parametri lagrangiani {f,g,f',g'}
per verificare se i parametri lagrangiani di sotto sono CORRETTI */
f= 1 - r2/p*2*sin(teta/2)*sin(teta/2);
g= r1*r2*sin(teta)/sqrt(MU*p);
fpunto= (2/p*sin(teta/2)*sin(teta/2) - 1/r1 - 1/r2)*tan(teta/2)*sqrt(MU/p);
gpunto= 1 - r1/p*2*sin(teta/2)*sin(teta/2);
//
V1L = sqrt((r2*r2 - 2*f*r1*r2*cos(teta) + f*f * r1*r1)/(g*g));
/* Vettore V2 = fpunto * r1 + gpunto * V1 */
sigma1= (r1*r2*cos(teta) - f*r1*r1)/g;
V2L= sqrt((fpunto*r1)*(fpunto*r1) + (gpunto*V1L)*(gpunto*V1L) + 2*fpunto*gpunto*sigma1);
/* OUTPUT finali */
printf("\n\n\t RISULTATI della CONICA (parametri di forma e velocita')\n");
printf("\t -------------------------------------------------------\n");
printf("\t\t a= semiasse magg.= %.10f UA\n\t\t p= param.focale= %.10f UA\n\t\t e= eccentricita'= %.10f\n", a,p,ecc);
printf("\n\t\tCalcolo delle velocita' con il metodo classico\n");
printf("\t\t V1= velocita' in P1= %.5f km/s\n\t\t V2= velocita' in P2= %.5f km/s\n\n", V1*uv0,V2*uv0);
printf("\t\tCalcolo delle velocita' con i parametri lagrangiani\n");
printf("\t\t V1= velocita' in P1= %.5f km/s\n\t\t V2= velocita' in P2= %.5f km/s\n", V1L*uv0,V2L*uv0);
printf("\t ------------------- Fine Calcoli ----------------------\n\n");
return (EXIT_SUCCESS);
}
Problema di LAMBERT: calcolo della traiettoria di trasferimento
di una sonda noti a) la posizione geometrica -> r1, r2, teta
b) il tempo di percorrenza -> t
[Moto Centrale Eliocentrico]
Raggio vettore r1 [UA], p.es. 1.0= 1
Raggio vettore r2 [UA], p.es. 1.8= 1.8
Angolo tra i due vettori teta [gradi], p.es. 109.884= 109.884
Tempo di percorrenza t [giorni], p.es. 253.86= 253.86
r1= 1.000 UA r2= 1.800 UA teta= 109.884 gradi delta-t= 253.860 giorni
Parametri intermedi
c= 2.33761010582490 UA s= 2.56880502907059 UA
lam=0.30000160296500 adim. Tadim=3.00001894683454 adim.
1) yit = 1.25093289 Residuo= 1.25e+000
2) yit = 1.24381067 Residuo= -7.12e-003
3) yit = 1.24380628 Residuo= -4.39e-006
4) yit = 1.24380628 Residuo= -1.68e-012
y-fin = 1.24380628 x-fin= 0.50885139
RISULTATI della CONICA (parametri di forma e velocita')
-------------------------------------------------------
a= semiasse magg.= 1.2853374112 UA
p= param.focale= 1.0495667071 UA
e= eccentricita'= 0.4282884337
Calcolo delle velocita' con il metodo classico
V1= velocita' in P1= 32.92514 km/s
V2= velocita' in P2= 17.19032 km/s
Calcolo delle velocita' con i parametri lagrangiani
V1= velocita' in P1= 32.92514 km/s
V2= velocita' in P2= 17.19032 km/s
------------------- Fine Calcoli ----------------------
(Luglio 2016). Raccolta dei miei lavori (works) anno per anno.
1997
File PDF da aprire o scaricare
Aggiornamento del codice di MOID in MatLab, prelevabile da questo file zippato e contenente il relativo DataBase. Di sotto sono riportati listato e output di un esempio applicativo.
% ---------- moid_new.m --------- 21 Luglio 2016
% Nota: la creazione della funzione m_iter.m è opera dell'amico
% e astrofilo Aldo Nicola, che ringrazio.
% ---------------------------------------------------------------------------
%
clc; clear all;
rad=pi/180;
format long
%
disp ('----------------------------------------------');
disp (' Calcolo del MOID ');
disp (' DataBase di 20,000 pianetini ');
disp ('----------------------------------------------');
disp ('');
%1^ parte: viene selezionata la riga del dBase, scartando la colonna con %*f
i=input('Numero di catalogo dell''asteroide [da 1 a 20,000]: ');
[perielio nodo inclinazione eccentricita semiassemaggiore nome] = ...
textread('DB20k.txt', '%f %f %f %f %*f %f %s',i);
%
nome(i)
fprintf ('\n omega=%10.5f Omega=%10.5f i=%9.5f e=%10.7f a=%10.7f\n', ...
perielio(i), nodo(i), inclinazione(i), eccentricita(i), semiassemaggiore(i));
disp('');
% ---------------- Fine lettura dei parametri del Dbase -----------
a= semiassemaggiore(i); e= eccentricita(i);
% ------- angoli in radianti --------
lnodo=nodo(i)*rad; om=perielio(i)*rad; incl=inclinazione(i)*rad;
%
% Terra (Fonte: Astronomical Algorithms di J.Meeus)
et = .01670862; lp = 102.937348*rad; at = 1.000001018;
%'****************************************************************************
p = a * (1 - e ^ 2); pt = at * (1 - et ^ 2);
%'****************************************************************************
% Angoli di EULERO {i, Omega, omega} --> Per Calcolo MATRICE di Trasformazione
R11 = cos(lnodo) * cos(om) - sin(lnodo) * cos(incl) * sin(om);
R12 = -cos(lnodo) * sin(om) - sin(lnodo) * cos(incl) * cos(om);
R21 = sin(lnodo) * cos(om) + cos(lnodo) * cos(incl) * sin(om);
R22 = -sin(lnodo) * sin(om) + cos(lnodo) * cos(incl) * cos(om);
R31 = sin(incl) * sin(om);
R32 = sin(incl) * cos(om);
%
% a) scelta della coppia dei valori iniziali (te,el) vicino al nodo ascendente
te = -om; el = lnodo;
% ----------------------------------------------
[x,y,z,xt,yt] = m_iter (R11,R21,R31,R12,R22,R32,pt,e,te,el,et,lp,p);
% ----------------------------------------------
moid1=sqrt((x-xt)^2+(y-yt)^2+z^2);
fprintf('\n d1= minima distanza tra le orbite nello spazio [nodo Asc]= %8.6f UA\n\n',moid1);
%
% b) scelta della coppia dei valori iniziali (te,el) vicino al nodo discendente
te = -om+pi; el = lnodo+pi;
% ----------------------------------------------
[x,y,z,xt,yt] = m_iter (R11,R21,R31,R12,R22,R32,pt,e,te,el,et,lp,p);
% ----------------------------------------------
moid2=sqrt((x-xt)^2+(y-yt)^2+z^2);
fprintf('\n d2= minima distanza tra le orbite nello spazio [nodo Disc]= %8.6f UA\n\n',moid2);
fprintf(' ----------------------------------------------------- \n');
if moid1 > moid2
moid = moid2;
else
moid = moid1;
end
%
fprintf(' MOID= Distanza minima Terra-Pianetino= %8.6f UA\n',moid);
fprintf(' ----------------------------------------------------- \n');
% ---- EOF: moid_new.m
function [x,y,z,xt,yt] = m_iter (R11,R21,R31,R12,R22,R32,pt,e,te,el,et,lp,p)
%
% --- file: m_iter.m da copiare nella stassa cartella di moid_new.m
%
ddte=1; ddel=1; % inizializzazione variabili di controllo
while (abs(ddte) > 1e-8) && (abs(ddel) > 1e-8)
%
x = (R11 * cos(te) + R12 * sin(te)) * p/(1+e*cos(te));
y = (R21 * cos(te) + R22 * sin(te)) * p/(1+e*cos(te));
z = (R31 * cos(te) + R32 * sin(te)) * p/(1+e*cos(te));
xt = pt * cos(el)/(1+et*cos(el-lp));
yt = pt * sin(el)/(1+et*cos(el-lp));
%----------------------------------------------------------------------------
Atil = -R11 * sin(te) + R12 * (cos(te) + e);
Btil = -R21 * sin(te) + R22 * (cos(te) + e);
Ctil = -R31 * sin(te) + R32 * (cos(te) + e);
Dtil = et * (-sin(el) * cos(el - lp) + cos(el) * sin(el - lp)) - sin(el);
Etil = et * (cos(el) * cos(el - lp) + sin(el) * sin(el - lp)) + cos(el);
%
A = Atil*p/(1+e*cos(te))^2;
B = Btil*p/(1+e*cos(te))^2;
C = Ctil*p/(1+e*cos(te))^2;
D = Dtil*pt/(1+et*cos(el-lp))^2;
E = Etil*pt/(1+et*cos(el-lp))^2;
%
Aast = -R11 * cos(te) - R12 * sin(te);
Bast = -R21 * cos(te) - R22 * sin(te);
Cast = -R31 * cos(te) - R32 * sin(te);
Dast = -cos(el);
East = -sin(el);
%-----------------------------------------------------------------------
F1 = Atil*(x-xt)+Btil*(y-yt)+Ctil*z;
F2 = Dtil*(x-xt)+Etil*(y-yt);
%
F11 = Aast*(x-xt) + Atil*A + Bast*(y-yt) + Btil*B + Cast*z + Ctil*C;
F12 = -Atil*D - Btil*E;
F21 = Dtil*A + Etil*B;
F22 = Dast*(x-xt) - Dtil*D + East*(y-yt) - Etil*E;
%
W0 = F11*F22 - F21*F12;
W1 = F1*F22 - F2*F12;
W2 = F11*F2 - F21*F1;
%
ddte= W1/W0; ddel= W2/W0;
te1= te - ddte; el1 = el - ddel; %Iterazione di Newton
%
fprintf ('\n te_i = %12.10f el_i = %12.10f', te1,el1);
te = te1; el=el1;
end
--------------------------------------------------------------------------
PHA's - Potential Hazardous Asteroids nell'intervallo [1-20,000]
--------------------------------------------------------------------------
(1566) Icarus (4034) Vishnu (6491) 1991 OA
(1620) Geographos (4179) Toutatis (7335) 1989 JA
(1862) Apollo (4183) Cuno (7341) 1991 VK
(1981) Midas (4450) Pan (7482) 1994 PC1
(2101) Adonis (4486) Mithra (7753) 1988 XB
(2102) Tantalus (4581) Asclepius (7822) 1991 CS
(2135) Aristaeus (4660) Nereus (8014) 1990 MF
(2201) Oljato (4769) Castalia (8566) 1996 EN
(2340) Hathor (4953) 1990 MU (9856) 1991 EE
(3122) Florence (5011) Ptah (10115) 1992 SK
(3200) Phaethon (5189) 1990 UQ (11500) Tomaiyowit
(3361) Orpheus (5604) 1992 FE (12538) 1998 OH
(3362) Khufu (5693) 1993 EA (12923) Zephyr
(3671) Dionysus (6037) 1988 EG (13651) 1997 BR
(3757) Anagolay (6239) Minos (14827) Hypnos
(4015) Wilson-Harrington (6489) Golevka (16960) 1998 QS52
--------------------------------------------------------------------------
----------------------------------------------
Calcolo del MOID
DataBase di 20,000 pianetini
----------------------------------------------
Numero di catalogo dell'asteroide [da 1 a 20,000]: 4179
ans =
'(4179)-Toutatis'
omega= 278.73541 Omega= 124.37315 i= 0.44711 e= 0.6294515 a= 2.5342052
te_i = -5.3537474518 el_i = 1.6797376225
te_i = -5.6123064511 el_i = 1.4225598679
te_i = -5.7380636521 el_i = 1.2973889389
te_i = -5.7782101638 el_i = 1.2573642509
te_i = -5.7826308079 el_i = 1.2529544594
te_i = -5.7826834982 el_i = 1.2529018907
te_i = -5.7826835057 el_i = 1.2529018833
d1= minima distanza tra le orbite nello spazio [nodo Asc]= 0.006111 UA
te_i = -1.0867007754 el_i = 5.9450770007
te_i = -0.7960267800 el_i = 6.2383106446
te_i = -0.6426328717 el_i = 6.3926840315
te_i = -0.5825931891 el_i = 6.4529380460
te_i = -0.5722515375 el_i = 6.4633071040
te_i = -0.5719498614 el_i = 6.4636095082
te_i = -0.5719496075 el_i = 6.4636097626
te_i = -0.5719496075 el_i = 6.4636097626
d2= minima distanza tra le orbite nello spazio [nodo Disc]= 0.007131 UA
-----------------------------------------------------
MOID= Distanza minima Terra-Pianetino= 0.006111 UA
-----------------------------------------------------
1997
File PDF da aprire o scaricare
Aggiornamento del codice delle Variabili Universali in C. Di sotto sono riportati listato (variabuniv.c) e output di un'orbita parabolica.
/* variabuniv.c, Rif. 1997b --- 07.Ago.2016 ---
* ' Dati iniziali (Elementi UNIVERSALI) nell'ordine: al, q, i, Om, om, tau
* */
#include ⟨stdio.h⟩
#include ⟨stdlib.h⟩
#include ⟨math.h⟩
#define PI 3.141592653589793238462643 /* in double conserva 15 cifre significative */
#define RAD 0.01745329251994329576923691 /* in double conserva 15 cifre significative */
#define KGAUSS 0.01720209895 /* in double mantiene le 11 cifre designate */
#define MU 1
int main(void)
{
float AL, Q, EI, BOM, OM, TAU;
double TOL = 1E-10;
double EIrad,BOMrad,OMrad, PX,PY,PZ, QX,QY,QZ, VP,XP,YP,ZP, VXP,VYP,VZP;
double PSIA,PSI, B,S3PARZ,S3, S2PARZ,S2, S1,S0;
double R,F,G, F1,G1, X,Y,Z, VX,VY,VZ, V;
/* Esempio orbita ELLITTICA
AL = -1.321604534; Q = .455635165; EI = 6.87631;
BOM = 63.07572; OM = 205.6752; TAU = 1.457528167;
*/
// Esempio orbita PARABOLICA
AL = 0; Q = 0.22432; EI = 122.639;
BOM = 188.943; OM = 131.202; TAU = -0.632503976;
/* Esempio orbita IPERBOLICA
AL = 4.98736E-4; Q = 0.555404; EI = 72.5488;
BOM = 237.8971; OM = 276.7690; TAU = 1.157986814;
*/
printf("\t\t ---------- I 6 Dati Orbitali Universali ------------\n");
printf(" alfa=(-mu/a)= %.8f q=dist.peri.= %.8f UA tau= t-Tp= %.8f TU\n",AL,Q,TAU);
printf(" Om= Long.Nodo= %.8f gr. w= Arg.peri.= %.8f gr. incl= %.8f gr.\n",BOM,OM,EI);
printf("\t\t ----------------------------------------------------\n\n");
EIrad= EI*RAD; BOMrad= BOM*RAD; OMrad= OM*RAD;
PX= cos(OMrad)*cos(BOMrad) - sin(OMrad)*cos(EIrad)*sin(BOMrad);
PY= cos(OMrad)*sin(BOMrad) + sin(OMrad)*cos(EIrad)*cos(BOMrad);
PZ= sin(OMrad)*sin(EIrad);
QX= -sin(OMrad)*cos(BOMrad) - cos(OMrad)*cos(EIrad)*sin(BOMrad);
QY= -sin(OMrad)*sin(BOMrad) + cos(OMrad)*cos(EIrad)*cos(BOMrad);
QZ= cos(OMrad)*sin(EIrad);
VP= sqrt(2*MU/Q + AL); /*Velocità al Perielio*/
XP= Q*PX; YP= Q*PY; ZP= Q*PZ;
VXP= VP*QX; VYP= VP*QY; VZP= VP*QZ;
// Risoluzione equazione di Keplero espressa con le V.U.
PSIA= TAU/Q; /* PSIA= valore iniziale dell'iterazione*/
for (;;) // ciclo infinito che è interrotto dal break
{
B= AL*PSIA*PSIA;
S3PARZ= ((((B/342 + 1)*B/272+1)*B/210 + 1)*B/156 + 1);
S3= PSIA*PSIA*PSIA/6*((((S3PARZ*B/110 + 1)*B/72 + 1)*B/42 + 1)*B/20 + 1);
S2PARZ= ((((B/306 + 1)*B/240 + 1)*B/182 + 1)*B/132 + 1);
S2= PSIA*PSIA/2*((((S2PARZ*B/90 + 1)*B/56 + 1)*B/30 + 1)*B/12 + 1);
S1 = PSIA + AL*S3;
S0 = 1 + AL*S2;
PSI= PSIA - (Q*S1 + MU*S3 - TAU)/(Q*S0 + MU*S2);
printf("\t PSI_iterato= %.10f\n",PSI);
if (fabs(PSI - PSIA) ⟨ TOL)
break;
else
PSIA= PSI;
} // fine ciclo infinito di for
printf("\n\t PSI_finale (PSIA)= %.10f\n",PSIA);
printf(" \t PSI_finale ( PSI)= %.10f\n\n",PSI);
/* Calcolo dei Vettori Posizione e Velocita'*/
R= Q*S0 + MU*S2; F= 1 - MU*S2/Q; G= TAU - MU*S3;
F1= -MU*S1/(R*Q); G1= 1 - MU*S2/R;
X= F*XP + G*VXP;
Y= F*YP + G*VYP;
Z= F*ZP + G*VZP;
VX= F1*XP + G1*VXP;
VY= F1*YP + G1*VYP;
VZ= F1*ZP + G1*VZP;
V= sqrt(VX*VX + VY*VY + VZ*VZ);
printf("\t X= %.8f Y= %.8f Z= %.8f --> R= %.8f UA\n",X,Y,Z,R);
printf("\t VX= %.8f VY= %.8f VZ= %.8f --> V= %.8f UA/UT\n",VX,VY,VZ,V);
return (EXIT_SUCCESS);
}
---------- I 6 Dati Orbitali Universali ------------
alfa=(-mu/a)= 0.00000000 q=dist.peri.= 0.22431999 UA tau= t-Tp= -0.63250399 TU
Om= Long.Nodo= 188.94299316 gr. w= Arg.peri.= 131.20199585 gr. incl= 122.63899994 gr.
----------------------------------------------------
PSI_iterato= -1.9299712990
PSI_iterato= -1.4514440116
PSI_iterato= -1.2927892425
PSI_iterato= -1.2761833321
PSI_iterato= -1.2760124511
PSI_iterato= -1.2760124331
PSI_iterato= -1.2760124331
PSI_finale (PSIA)= -1.2760124331
PSI_finale ( PSI)= -1.2760124331
X= -1.02901223 Y= -0.09682577 Z= 0.10041273 --> R= 1.03842386 UA
VX= 1.23710597 VY= 0.46747707 VZ= 0.42074911 --> V= 1.38780251 UA/UT
2002
File PDF da aprire o scaricare
2002
File PDF da aprire o scaricare
2003
File PDF da aprire o scaricare
2003
File PDF da aprire o scaricare
2003
File PDF da aprire o scaricare
2004
File PDF da aprire o scaricare
2004
File PDF da aprire o scaricare
2004
File PDF da aprire o scaricare
File PDF da aprire o scaricare
2004
2004
File PDF da aprire o scaricare
2004
File PDF da aprire o scaricare
2005
File PDF da aprire o scaricare
2005
File PDF da aprire o scaricare
2005
File PDF da aprire o scaricare
Aggiornamento del codice degli elemnti orbitali noti 2 Vettori e l'angolo con la tangente in C. Di sotto sono riportati listato (orbit-r1-r2.c) e output.
/* orbit-r1-r2.c :: Riferimento paper 2005c.pdf ---- 02-Ago-2016 --- */
#include ⟨stdio.h⟩
#include ⟨stdlib.h⟩
#include ⟨math.h⟩
#define PI 3.141592653589793238462643 /* in double conserva 15 cifre significative */
#define RAD 0.01745329251994329576923691 /* in double conserva 15 cifre significative */
#define KGAUSS 0.01720209895 /* in double mantiene le 11 cifre designate */
#define MU 1
int main(void)
{
/* Dati i vettori (r1,r2) con coordinate in float */
float x1 = -0.106418; float y1 = 0.137154; float z1 = 1.637343;
float x2 = -2.60002887; float y2 = 1.62023766; float z2 = 2.21048897;
float be = 63.54333316; /* angolo beta in gradi */
/* Inizio Calcoli*/
double r1= sqrt(x1*x1 + y1*y1 + z1*z1);
double r2= sqrt(x2*x2 + y2*y2 + z2*z2);
double cosalfa = (x1*x2 + y1*y2 + z1*z2)/(r1*r2);
double alfa = acos(cosalfa);
/* Costanti */
double K2 = cosalfa; double K1 = sqrt(1-K2*K2);
double B = r1 * sin(be*RAD); double K0 = 1/(B*B);
double K4 = r1/r2 + K1*sqrt(K0*r1*r1 - 1);
/* Calcolo di x=e*COS(te1), y=e^2, z=p */
double x = (1 - K4)/(K4 - K2);
double y = K0*r1*r1 * (1+x)*(1+x) - 2*x - 1;
double z = r1*(1+x);
double ecc = sqrt(y);
double p = z;
double a = p/(1-ecc*ecc);
/* Calcolo Anomalie (te1,te2) */
double te1 = acos(x/ecc);
double te2 = te1 + alfa;
/* Calcolo dei Versori Perifocali P=(Px,Py,Pz) e Q=(Qx,Qy,Qz) */
double a1= r1*cos(te1); double b1= r1*sin(te1);
double a2= r2*cos(te2); double b2= r2*sin(te2);
double Det= a1*b2 - b1*a2;
double Px= (x1*b2 - x2*b1)/ Det;
double Py= (y1*b2 - y2*b1)/ Det;
double Pz= (z1*b2 - z2*b1)/ Det;
double Qx= (a1*x2 - a2*x1)/ Det;
double Qy= (a1*y2 - a2*y1)/ Det;
double Qz= (a1*z2 - a2*z1)/ Det;
/* Calcolo delle Velocità (V1,V2) */
double V1= sqrt(MU*(1 + 2*ecc*cos(te1) + ecc*ecc)/p);
double V2= sqrt(MU*(1 + 2*ecc*cos(te2) + ecc*ecc)/p);
/* Calcolo degli elementi orbitali (i,Omega,w) */
double om= atan2(Pz,Qz);
if (om ⟨ 0)
om = om + 2*PI;
double num1= Py*cos(om) - Qy*sin(om);
double den1= Px*cos(om) - Qx*sin(om);
double nodo= atan2(num1,den1);
if (nodo ⟨ 0)
nodo = nodo + 2*PI;
double num2= Pz/sin(om);
double den2= -(Px*sin(om) + Qx*cos(om))/sin(nodo);
double incl= atan2(num2,den2);
/* OUTPUTS */
printf("\t Modulo Vettore r1= %.8f UA\n",r1);
printf("\t Modulo Vettore r2= %.8f UA\n\n",r2);
printf("\t Angolo tra i due Vettori= %.8f gradi \n\n",alfa/RAD);
printf("\t a= %.8f UA ecc= %.8f p= %.8f UA\n\n",a,ecc,p);
if (a ⟨ 0)
printf("\t La Traiettoria e' IPERBOLICA\n\n");
else
printf("\t La Traiettoria e' ELLITTICA\n\n");
printf("\t Anomalia-1= %.8f gradi\n",te1/RAD);
printf("\t Anomalia-2= %.8f gradi\n\n",te2/RAD);
printf("\t Arg.Peri(w)= %.8f gradi\n",om/RAD);
printf("\t Long.Nodo(O)= %.8f gradi\n",nodo/RAD);
printf("\t Inclinaz.(i)= %.8f gradi\n\n",incl/RAD);
printf("\t\t Fine Calcoli \n");
return (EXIT_SUCCESS);
}
Modulo Vettore r1= 1.64652005 UA
Modulo Vettore r2= 3.77777468 UA
Angolo tra i due Vettori= 48.54151166 gradi
a= -1.88461051 UA ecc= 1.73559597 p= 3.79238923 UA
La Traiettoria e' IPERBOLICA
Anomalia-1= 41.33077865 gradi
Anomalia-2= 89.87229031 gradi
Arg.Peri(w)= 54.28322750 gradi
Long.Nodo(O)= 329.70534118 gradi
Inclinaz.(i)= 87.73564140 gradi
2005
File PDF da aprire o scaricare
2005
File PDF da aprire o scaricare
2008
File PDF da aprire o scaricare
2012
File PDF da aprire o scaricare
Aggiornamento del codice delle Variabili Universali in C. Di sotto sono riportati listato (variabuniv.c) e output di un'orbita parabolica.
/* variabuniv.c, Rif. 1997b --- 07.Ago.2016 ---
* ' Dati iniziali (Elementi UNIVERSALI) nell'ordine: al, q, i, Om, om, tau
* */
#include ⟨stdio.h⟩
#include ⟨stdlib.h⟩
#include ⟨math.h⟩
#define PI 3.141592653589793238462643 /* in double conserva 15 cifre significative */
#define RAD 0.01745329251994329576923691 /* in double conserva 15 cifre significative */
#define KGAUSS 0.01720209895 /* in double mantiene le 11 cifre designate */
#define MU 1
int main(void)
{
float AL, Q, EI, BOM, OM, TAU;
double TOL = 1E-10;
double EIrad,BOMrad,OMrad, PX,PY,PZ, QX,QY,QZ, VP,XP,YP,ZP, VXP,VYP,VZP;
double PSIA,PSI, B,S3PARZ,S3, S2PARZ,S2, S1,S0;
double R,F,G, F1,G1, X,Y,Z, VX,VY,VZ, V;
/* Esempio orbita ELLITTICA
AL = -1.321604534; Q = .455635165; EI = 6.87631;
BOM = 63.07572; OM = 205.6752; TAU = 1.457528167;
*/
// Esempio orbita PARABOLICA
AL = 0; Q = 0.22432; EI = 122.639;
BOM = 188.943; OM = 131.202; TAU = -0.632503976;
/* Esempio orbita IPERBOLICA
AL = 4.98736E-4; Q = 0.555404; EI = 72.5488;
BOM = 237.8971; OM = 276.7690; TAU = 1.157986814;
*/
printf("\t\t ---------- I 6 Dati Orbitali Universali ------------\n");
printf(" alfa=(-mu/a)= %.8f q=dist.peri.= %.8f UA tau= t-Tp= %.8f TU\n",AL,Q,TAU);
printf(" Om= Long.Nodo= %.8f gr. w= Arg.peri.= %.8f gr. incl= %.8f gr.\n",BOM,OM,EI);
printf("\t\t ----------------------------------------------------\n\n");
EIrad= EI*RAD; BOMrad= BOM*RAD; OMrad= OM*RAD;
PX= cos(OMrad)*cos(BOMrad) - sin(OMrad)*cos(EIrad)*sin(BOMrad);
PY= cos(OMrad)*sin(BOMrad) + sin(OMrad)*cos(EIrad)*cos(BOMrad);
PZ= sin(OMrad)*sin(EIrad);
QX= -sin(OMrad)*cos(BOMrad) - cos(OMrad)*cos(EIrad)*sin(BOMrad);
QY= -sin(OMrad)*sin(BOMrad) + cos(OMrad)*cos(EIrad)*cos(BOMrad);
QZ= cos(OMrad)*sin(EIrad);
VP= sqrt(2*MU/Q + AL); /*Velocità al Perielio*/
XP= Q*PX; YP= Q*PY; ZP= Q*PZ;
VXP= VP*QX; VYP= VP*QY; VZP= VP*QZ;
// Risoluzione equazione di Keplero espressa con le V.U.
PSIA= TAU/Q; /* PSIA= valore iniziale dell'iterazione*/
for (;;) // ciclo infinito che è interrotto dal break
{
B= AL*PSIA*PSIA;
S3PARZ= ((((B/342 + 1)*B/272+1)*B/210 + 1)*B/156 + 1);
S3= PSIA*PSIA*PSIA/6*((((S3PARZ*B/110 + 1)*B/72 + 1)*B/42 + 1)*B/20 + 1);
S2PARZ= ((((B/306 + 1)*B/240 + 1)*B/182 + 1)*B/132 + 1);
S2= PSIA*PSIA/2*((((S2PARZ*B/90 + 1)*B/56 + 1)*B/30 + 1)*B/12 + 1);
S1 = PSIA + AL*S3;
S0 = 1 + AL*S2;
PSI= PSIA - (Q*S1 + MU*S3 - TAU)/(Q*S0 + MU*S2);
printf("\t PSI_iterato= %.10f\n",PSI);
if (fabs(PSI - PSIA) ⟨ TOL)
break;
else
PSIA= PSI;
} // fine ciclo infinito di for
printf("\n\t PSI_finale (PSIA)= %.10f\n",PSIA);
printf(" \t PSI_finale ( PSI)= %.10f\n\n",PSI);
/* Calcolo dei Vettori Posizione e Velocita'*/
R= Q*S0 + MU*S2; F= 1 - MU*S2/Q; G= TAU - MU*S3;
F1= -MU*S1/(R*Q); G1= 1 - MU*S2/R;
X= F*XP + G*VXP;
Y= F*YP + G*VYP;
Z= F*ZP + G*VZP;
VX= F1*XP + G1*VXP;
VY= F1*YP + G1*VYP;
VZ= F1*ZP + G1*VZP;
V= sqrt(VX*VX + VY*VY + VZ*VZ);
printf("\t X= %.8f Y= %.8f Z= %.8f --> R= %.8f UA\n",X,Y,Z,R);
printf("\t VX= %.8f VY= %.8f VZ= %.8f --> V= %.8f UA/UT\n",VX,VY,VZ,V);
return (EXIT_SUCCESS);
}
---------- I 6 Dati Orbitali Universali ------------
alfa=(-mu/a)= 0.00000000 q=dist.peri.= 0.22431999 UA tau= t-Tp= -0.63250399 TU
Om= Long.Nodo= 188.94299316 gr. w= Arg.peri.= 131.20199585 gr. incl= 122.63899994 gr.
----------------------------------------------------
PSI_iterato= -1.9299712990
PSI_iterato= -1.4514440116
PSI_iterato= -1.2927892425
PSI_iterato= -1.2761833321
PSI_iterato= -1.2760124511
PSI_iterato= -1.2760124331
PSI_iterato= -1.2760124331
PSI_finale (PSIA)= -1.2760124331
PSI_finale ( PSI)= -1.2760124331
X= -1.02901223 Y= -0.09682577 Z= 0.10041273 --> R= 1.03842386 UA
VX= 1.23710597 VY= 0.46747707 VZ= 0.42074911 --> V= 1.38780251 UA/UT
(06-Lug-13). L'uploading dei miei lavori finisce qui
