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

Astronomia da Almanacco

(16-Mag-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) Equinozi/Solstizi, ii) Fasi Lunari, iii) 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, quando per gli equinozi la longitudine del Sole Apparente sia nulla o pari a 180 gradi e per i solstizi 90/270 gradi; nel secondo 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 terzo 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.

Primo Esempio

# --------- equisol.py --------- Equinoxes and Solstices - Nov.2019 ----
# ----- Istanti in Temi Universali (UT) ------
import ephem
d1 = ephem.next_equinox('2020')
print ('\n Equinozio=', d1)
d2 = ephem.next_solstice(d1)
print (' Solstizio=', d2)
d3 = ephem.next_equinox(d2)
print (' Equinozio=', d3)
d4 = ephem.next_solstice(d3)
print (' Solstizio=', d4)
#
d5= ephem.next_vernal_equinox('2021')
print ('\n Equinozio successivo= ', d5)

# EOF: -------- Equisol.py -------
'''
---------- OUTPUT --------------
 Equinozio= 2020/3/20 03:49:33
 Solstizio= 2020/6/20 21:43:47
 Equinozio= 2020/9/22 13:30:38
 Solstizio= 2020/12/21 10:02:09

 Equinozio successivo=  2021/3/20 09:37:26

'''	 

Secondo esempio

# fasi_luna.py    -------- Jan.2020 ----------
# ----- Istanti in Temi Universali (UT) ------
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

'''	

Terzo esempio

# --------- crepus2.py ----------  Apr/May.2020 ----------
# Calcolo dei 3 tipi di crepuscolo: civile, nautico e astronomico
# oltre a sorgere, transito, tramonto Sole
# Eliminate le parentesi () di print per far girare lo script su Python 2.7
# Avviare con: py -2.7 crepus2.py
import ephem

#Make an observer
paris      = ephem.Observer()

#PyEphem takes and returns only UTC times. Paris (France)
paris.date = "1961-08-22"

#Location of paris
paris.lon  = str(2.330833) #Note that lon should be in string format
paris.lat  = str(48.868889) #Note that lat should be in string format

#Elevation of paris in metres
paris.elev = 32

#To get U.S. Naval Astronomical Almanac values, use these settings
paris.pressure= 0
paris.horizon = '-0:34'

sunrise=paris.previous_rising(ephem.Sun()) #Sunrise
noon   =paris.next_transit   (ephem.Sun(), start=sunrise) #Solar noon
sunset =paris.next_setting   (ephem.Sun()) #Sunset

#We relocate the horizon to get twilight times
paris.horizon = '-18' #-6=civil twilight, -12=nautical, -18=astronomical
beg_twilight=paris.previous_rising(ephem.Sun(), use_center=True) #Begin civil twilight
end_twilight=paris.next_setting   (ephem.Sun(), use_center=True) #End civil twilight
end1=end_twilight
print' Crepuscolo astronomico mattino = ',beg_twilight, 'UT'
#
paris.horizon = '-12' #-6=civil twilight, -12=nautical, -18=astronomical
beg_twilight=paris.previous_rising(ephem.Sun(), use_center=True) #Begin civil twilight
end_twilight=paris.next_setting   (ephem.Sun(), use_center=True) #End civil twilight
print' Crepuscolo     nautico mattino = ',beg_twilight, 'UT'
end2=end_twilight
#
paris.horizon = '-6' #-6=civil twilight, -12=nautical, -18=astronomical
beg_twilight=paris.previous_rising(ephem.Sun(), use_center=True) #Begin civil twilight
end_twilight=paris.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: crepus2.py
# -------- OUTPUT ----------
# 
# Crepuscolo astronomico mattino =  1961/8/21 02:46:42 UT
# Crepuscolo     nautico mattino =  1961/8/21 03:34:29 UT
# Crepuscolo      civile mattino =  1961/8/21 04:16:58 UT
#
#                  Sole    sorge =  1961/8/21 04:51:10 UT
#                  Sole transita =  1961/8/21 11:53:48 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.  confronto calcoli eclisse

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:  GitHub link 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.

 star altaz graph

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  1997a-work

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  1997b-work

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  2002a-work 2002 File PDF da aprire o scaricare  2002b-work

2003 File PDF da aprire o scaricare  2003a-work 2003 File PDF da aprire o scaricare  2003b-work 2003 File PDF da aprire o scaricare  2003c-work

2004 File PDF da aprire o scaricare  2004a-work 2004 File PDF da aprire o scaricare  2004b-work 2004 File PDF da aprire o scaricare  2004c-work File PDF da aprire o scaricare 2004  2004d-work 2004 File PDF da aprire o scaricare  2004e-work 2004 File PDF da aprire o scaricare  2004f-work

2005 File PDF da aprire o scaricare  2005a-work 2005 File PDF da aprire o scaricare  2005b-work 2005 File PDF da aprire o scaricare  2005c-work

Aggiornamento del codice degli elemnti orbitali noti 2 Vettori e l'angolo con la tangente in C. Di sotto sono riportati listato (orbit-r1-r2.c) e output.

/* orbit-r1-r2.c ::  Riferimento paper 2005c.pdf ---- 02-Ago-2016 --- */
#include ⟨stdio.h⟩
#include ⟨stdlib.h⟩
#include ⟨math.h⟩

#define PI 3.141592653589793238462643     /* in double conserva 15 cifre significative */
#define RAD 0.01745329251994329576923691  /* in double conserva 15 cifre significative */
#define KGAUSS 0.01720209895              /* in double mantiene le 11 cifre designate */
#define MU 1

int main(void)
{
	/* Dati i vettori (r1,r2) con coordinate in float */
    float x1 = -0.106418; float y1 = 0.137154; float z1 = 1.637343;
    float x2 = -2.60002887; float y2 = 1.62023766; float z2 = 2.21048897;
	float be = 63.54333316;  /* angolo beta in gradi */
	
	/* Inizio Calcoli*/
	double r1= sqrt(x1*x1 + y1*y1 + z1*z1);
    double r2= sqrt(x2*x2 + y2*y2 + z2*z2);
  
    double cosalfa = (x1*x2 + y1*y2 + z1*z2)/(r1*r2);
    double alfa = acos(cosalfa);
    
    /* Costanti */
    double K2 = cosalfa;  double K1 = sqrt(1-K2*K2);
    double B = r1 * sin(be*RAD); double K0 = 1/(B*B);
    double K4 = r1/r2 + K1*sqrt(K0*r1*r1 - 1);
    
    /* Calcolo di x=e*COS(te1), y=e^2, z=p */
    double x = (1 - K4)/(K4 - K2);
    double y = K0*r1*r1 * (1+x)*(1+x) - 2*x - 1;
    double z = r1*(1+x);
    
    double ecc = sqrt(y);
    double p = z;
    double a = p/(1-ecc*ecc);
    
    /* Calcolo Anomalie (te1,te2) */
    double te1 = acos(x/ecc);
    double te2 = te1 + alfa; 

    /* Calcolo dei Versori Perifocali P=(Px,Py,Pz) e Q=(Qx,Qy,Qz) */
    double a1= r1*cos(te1); double b1= r1*sin(te1);
    double a2= r2*cos(te2); double b2= r2*sin(te2);
    
    double Det= a1*b2 - b1*a2;
    
    double Px= (x1*b2 - x2*b1)/ Det;
    double Py= (y1*b2 - y2*b1)/ Det;
    double Pz= (z1*b2 - z2*b1)/ Det;
    
    double Qx= (a1*x2 - a2*x1)/ Det;
    double Qy= (a1*y2 - a2*y1)/ Det;
    double Qz= (a1*z2 - a2*z1)/ Det; 
    
    /* Calcolo delle Velocità (V1,V2) */
    double V1= sqrt(MU*(1 + 2*ecc*cos(te1) + ecc*ecc)/p);
    double V2= sqrt(MU*(1 + 2*ecc*cos(te2) + ecc*ecc)/p);

    /* Calcolo degli elementi orbitali (i,Omega,w) */
    double om= atan2(Pz,Qz);
	if (om ⟨ 0)
		om = om + 2*PI; 
    
    double num1= Py*cos(om) - Qy*sin(om);
    double den1= Px*cos(om) - Qx*sin(om);
    double nodo= atan2(num1,den1);    
    if (nodo ⟨ 0)
		nodo = nodo + 2*PI; 
		
	double num2= Pz/sin(om);
	double den2= -(Px*sin(om) + Qx*cos(om))/sin(nodo);	
    double incl= atan2(num2,den2);
    
/*   OUTPUTS */       
  
    printf("\t Modulo Vettore r1= %.8f UA\n",r1);   
    printf("\t Modulo Vettore r2= %.8f UA\n\n",r2);   
    printf("\t Angolo tra i due Vettori= %.8f gradi \n\n",alfa/RAD);
   
    printf("\t a= %.8f UA  ecc= %.8f   p= %.8f UA\n\n",a,ecc,p);
    
	if (a ⟨ 0)
		printf("\t La Traiettoria e' IPERBOLICA\n\n"); 
	else
		printf("\t La Traiettoria e' ELLITTICA\n\n"); 
    
    printf("\t Anomalia-1= %.8f gradi\n",te1/RAD);
    printf("\t Anomalia-2= %.8f gradi\n\n",te2/RAD);
    
    printf("\t  Arg.Peri(w)= %.8f gradi\n",om/RAD);
    printf("\t Long.Nodo(O)= %.8f gradi\n",nodo/RAD);
    printf("\t Inclinaz.(i)= %.8f gradi\n\n",incl/RAD);
    
    printf("\t\t  Fine Calcoli \n");
    return (EXIT_SUCCESS);

}
	  
	
         Modulo Vettore r1= 1.64652005 UA
         Modulo Vettore r2= 3.77777468 UA

         Angolo tra i due Vettori= 48.54151166 gradi

         a= -1.88461051 UA  ecc= 1.73559597   p= 3.79238923 UA

         La Traiettoria e' IPERBOLICA

         Anomalia-1= 41.33077865 gradi
         Anomalia-2= 89.87229031 gradi

          Arg.Peri(w)= 54.28322750 gradi
         Long.Nodo(O)= 329.70534118 gradi
         Inclinaz.(i)= 87.73564140 gradi

	
2005 File PDF da aprire o scaricare  2005d-work 2005 File PDF da aprire o scaricare  2005e-work

2008 File PDF da aprire o scaricare  2008a-work

2012 File PDF da aprire o scaricare  2012a-work

1997 File PDF da aprire o scaricare  1997b-work

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


Valid HTML 4.01 Strict