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

(13 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