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
Equinozi e Solstizi in Python con arrotondamento dei decimali
Un piccolo neo contraddistingue la potente funzione datetime di Python-3, la stringa dei tempi può essere tagliata fino ai microsecondi, come ho fatto in basso nella sezione Astronomia da Almanacco (secondo esempio), ma non arrotondata ai decimali voluti. Il problema é stato risolto grazie al contributo di Giosue' Ruscica, a cui mi legano stretti vincoli di parentela.
# ---- equi_solst.py ------- 03-Gen-2021 ------ # --- script reso _professionale_ da Giosue' Ruscica, che ringrazio --- # # Riguarda l'arrotondamento dei decimali della stringa temporale # di (datetime), che normalmente l'opzione '{:.2f}'.format() visualizza # con i decimali di calcolo tagliati invece che arrotondati. # import ephem import datetime as dt DECIMALS = 2 DATE_FMT = '%Y-%m-%d %H:%M:%S.%f' # lasciare sempre %f finale def round_dt_string(t, ndigits=DECIMALS, fmt=DATE_FMT): t_stamp = t.timestamp() r = dt.datetime.fromtimestamp(round(t_stamp, ndigits)) if ndigits < 1: ndigits = -1 right_end = -(6 - ndigits) if ndigits < 6 else None return r.strftime(fmt)[:right_end] if __name__ == "__main__": a = input("\n Scegliere l'anno: ") d1 = ephem.next_equinox(a) print('\n Equinozio=', round_dt_string(d1.datetime())) d2 = ephem.next_solstice(d1) print(' Solstizio=', round_dt_string(d2.datetime())) d3 = ephem.next_equinox(d2) print(' Equinozio=', round_dt_string(d3.datetime())) d4 = ephem.next_solstice(d3) print(' Solstizio=', round_dt_string(d4.datetime())) d5 = ephem.next_vernal_equinox(str(int(a)+1)) print('\n Equinozio successivo= ', round_dt_string(d5.datetime())) # EOF: equi_solst.py ''' C:\Training>python equi_solst.py Scegliere l'anno: 2021 Equinozio= 2021-03-20 09:37:26.12 Solstizio= 2021-06-21 03:32:14.90 Equinozio= 2021-09-22 19:21:08.96 Solstizio= 2021-12-21 15:59:05.44 Equinozio successivo= 2022-03-20 15:33:19.38 '''
Riporto qui appresso l'output dello script in basso, con i solstizi a decimali 'tagliati'
# EOF: -------- Equisol.py ------- ''' Equinozio= 2021-03-20 09:37:26.12 Solstizio= 2021-06-21 03:32:14.89 Equinozio= 2021-09-22 19:21:08.96 Solstizio= 2021-12-21 15:59:05.43 Equinozio successivo= 2022-03-20 15:33:19.37 '''
Istante quando una parete di un edificio rimane in ombra (script Python)
Tra le applicazioni dell'astronomia di posizione che in passato avrei voluto risolvere c'è quella presentata in questo paragrafo. La limitazione era dovuta all'impossibilità di conoscere con ottima approssimazione l'azimut della parete di riferimento (muro del terrazzo di casa) su cui prima batte il sole e dopo, nel corso della giornata, va in ombra. Ebbene, la conoscenza dell'istante in cui il centro del sole oltrepassa l'allineamento della parete adesso si può risolvere con uno strumento grafico efficiente invece che con un teodolite, grazie alle mappe telematiche disponibili nel web e processate con il formidabile (e gratuito) pacchetto Geogebra di geometria dinamica, come nella seguente figura.

Il coefficiente angolare della retta che congiunge due punti della proiezione orizzontale della parete definisce l'angolo che essa forma con l'asse x e il relativo azimut. Quindi il passo successivo è quello di disporre di un programma computerizzato (script) che calcoli l'istante di passaggio del sole al raggiungimento di tale azimut. Il linguaggio Python dispone dell'eccellente funzione fsolve() che calcola (fino alle frazioni di secondo) il tempo quando si annulla la differenza (az_sole-az_parete). Con un copia-e-incolla dello script riportato sotto il lettore può eseguire il calcolo cercato inserendovi i dati necessari commentati all'interno del listato.
# ------- azi_parete_jd.py ---------- Dec.18, 2020 --------- # --- calcolo dell'azimut del Sole quando uguaglia quello di un allineamento assegnato (parete) ---- from scipy.optimize import fsolve import ephem from datetime import datetime, timedelta from math import floor, pi # def calendario(jd): z= floor(jd+0.5); f=jd+0.5-z if z>= 2299161: alfa=floor((z-1867216.25)/36524.25) a= z+1+alfa-floor(alfa/4) else: a=z b= a+1524; c= floor((b-122.1)/365.25) d= floor(365.25*c); e= floor((b-d)/30.6001) giorno= b-d-floor(30.6001*e)+f # if e<13.5: mese= e-1 else: mese= e-13 # if mese>2.5: anno= c-4716 else: anno= c-4715 # return giorno, mese, anno def equation(p): x = p az0=205.15 # azimut della parete di riferimento sun = ephem.Sun() observer = ephem.Observer() # ------ Coordinate geografiche (longitudine, latitudine) in gradi, altezza slm in metri observer.lon = '15.0617' observer.lat = '37.0350' observer.elevation = 360 # observer.date = x sun.compute(observer) azi= float(sun.az) * 180.0/pi return (azi-az0) # ---- Inserimento di anno, mese, giorno (lasciare immutati i tempi) ---------- x0=ephem.Date('2021/3/30 11:0:0') x = fsolve(equation, x0) # jd_converge= float(x)+2415020.0 # JD=2415020.0 = inizio scala dei tempi in JD # --- Istante in cui il centro del Sole raggiunge l'azimut della parete --- print("\nIstante convergenza(jd)=",jd_converge) print(' Data civile=', calendario(jd_converge)) jd = float(x) + 2415020.0 day = calendario(jd)[0] # frazd= day - floor(day); orefraz= frazd*24.0; oreint= floor(orefraz) minfraz= (orefraz-oreint)*60.0; minint=floor(minfraz) sec=(minfraz-minint)*60 print(' T_Univ(ore, min,sec)=',oreint, minint,'{:.0f}'.format(sec)) print ('\nFormat unico (Data+Tempo_Univ):',ephem.Date(x)) # (x) utilizza la scala dei tempi di Dublino (0= 1899-12-31,12:00:00) # # EOF: azi_parete_jd.py # ''' Istante convergenza(jd)= 2459304.0010969983 Data civile= (30.50109699834138, 3, 2021) T_Univ(ore, min,sec)= 12 1 35 Format unico (Data+Tempo_Univ): 2021/3/30 12:01:35 '''
Istanti riportati in tabella
Ad intervalli desiderati, come nello script seguente, da inizio a fine data con la scansione indicata, per es. timedelta(days = 10). Il passaggio del Sole al meridiano si ottiene settando a 180 gradi l'azimut az0, come riportato sotto.
# -------- sun_transit_terrazzo.py ----------- 20-26.12.2020 ---------- # --- azimut parete posto a 180°, passaggio del Sole al meridiano ---- from scipy.optimize import fsolve import ephem from datetime import datetime, timedelta from math import pi def equation(p): x = p az0=180.0 # azimut del meridiano locale sun = ephem.Sun() observer = ephem.Observer() # ------ Coordinate geografiche (longitudine, latitudine) in gradi, altezza slm in metri observer.lon = '15.0650' observer.lat = '37.0328' observer.elevation = 370 # observer.date = x sun.compute(observer) azi= float(sun.az) * 180.0/pi return (azi-az0) curtime = datetime(2021, 1, 20, 11, 0, 0) # start time, Modificare solo la data, NON le ore che servono per la convergenza endtime = datetime(2021, 4, 20, 11, 0, 0) # end time, come sopra print('\nOre in tempo universale (TU)\n-----------------------') while curtime <= endtime: x0=ephem.Date(curtime) x = fsolve(equation, x0) stringa= ephem.Date(x) string2= stringa.datetime() # serve per poter eliminare le frazioni di secondo print (str(string2)[0:22]) # a 2 decimali di secondo per mostrare la validita' della manipolazione # settare a [0:19] per i secondi senza decimali #print(stringa.datetime()) # output a 6 'irrealistici' decimali di secondo #print(stringa.datetime().strftime("%Y/%m/%d %H:%M:%S")) # taglia le frazioni di secondo curtime += timedelta(days = 10) # timedelta parameters: days, hours, minutes, seconds # # EOF: sun_transit_terrazzo.py ''' Ore in tempo universale (TU) ----------------------- 2021-01-20 11:10:49.23 2021-01-30 11:13:02.53 2021-02-09 11:13:55.61 2021-02-19 11:13:32.09 2021-03-01 11:12:01.28 2021-03-11 11:09:41.51 2021-03-21 11:06:51.81 2021-03-31 11:03:50.74 2021-04-10 11:01:00.04 2021-04-20 10:58:37.42 '''
(16-26.Dic.2020)
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 alcuni corpi, che per maggiore praticità sono stati compressi qui. All'avvio appare un prompt che indica all'utente il nome del file dati da scegliere e l'opzione di ricerca di altre radici. Si preme invio per ottenere quella certa, per le altre basta fare riferimento al riepilogo in alto, dove sono suggeriti i valori da digitare
Esempio 1: Pianetino 1924-TD-Baade (1 radice, starting value generalizzato) Esempio 2: B1950.0 Cometa Arend-Rigaux 1950 VII (3 radici, starting values: 0.2, 1.3, 2.0)
# --- 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)
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 - Dic.2020 ---- # ----- Istanti in Tempi Universali (UT) ------ # -------- Python3 environment --------- # Windows: python equisol.py # Linux: python3 equisol.py import ephem d1 = ephem.next_equinox('2021') stringa= d1.datetime(); string2=str(stringa)[0:22] print ('\n Equinozio=', string2) d2 = ephem.next_solstice(d1) stringa= d2.datetime(); string2=str(stringa)[0:22] print (' Solstizio=', string2) d3 = ephem.next_equinox(d2) stringa= d3.datetime(); string2=str(stringa)[0:22] print (' Equinozio=', string2) d4 = ephem.next_solstice(d3) stringa= d4.datetime(); string2=str(stringa)[0:22] print (' Solstizio=', string2) # d5= ephem.next_vernal_equinox('2022') stringa= d5.datetime(); string2=str(stringa)[0:22] print ('\n Equinozio successivo= ', string2) # EOF: -------- Equisol.py ------- ''' Equinozio= 2021-03-20 09:37:26.12 Solstizio= 2021-06-21 03:32:14.89 Equinozio= 2021-09-22 19:21:08.96 Solstizio= 2021-12-21 15:59:05.43 Equinozio successivo= 2022-03-20 15:33:19.37 '''
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 gradi2005 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
