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-16 Dicembre 2016). Operazione di restyling del codice asteroidale aster16.c pubblicato sotto. Per la lettura delle stringhe del database non viene usato il laborioso comando fseek ma fscanf che si avvale di un file di deposito esterno per copiarvi (in una riga) i dati del record intercettato. Ciò permette di modificare rapidamente il codice C in caso di variazione del template dei campi (larghezza e numero). L'algoritmo è stato provato sul file AsterSelect.c ed applicato successivamente allo script finale AsterDIC16.c, entrambi reperibili qui, insieme allo stesso dbase di 20,000 pianetini.
(28 Settembre 2016). Una interessante variante del software asteroidale è questa. I due Database degli elementi orbitali, richiamati dai codici AstGeocP.c e AstNEAp.c, contengono rispettivamente 20,000 asteroidi del catalogo MPC e 14,690 NEA, acronimo di Near Earth Asteroids; essi vengono scandagliati scegliendo come parametro di screening la distanza geocentrica; per una migliore consultazione, i risultati ottenuti sono copiati su file testo. Nel caso dei NEA, si percepisce subito quali siano i più vicini alla Terra alla data dell'effemeride, come nell'esempio seguente:
Distanza Geocentrica limite assunta= 0.02500 Data Effemeride= 2016 9 24 2016 RM20 Distanza GEOCENTRICA d= 0.01991001 (UA) Distanza Eliocentr. r= 1.01175656 (UA) Long.Eclitt.Eliocen. l= 1.07918 (gradi) Lati.Eclitt.Eliocen. b= -1.01969 (gradi) Ascensione Retta AR= 36.12525 (gradi) = 2H 24M 30S Declinazione DE= -58.36825 (gradi) psi= -114.98831 (gradi) Elongaz. Ovest Angolo di Fase beta= 63.98968 (gradi) Magnitudine= 20.0 2016 RD34 Distanza GEOCENTRICA d= 0.01327780 (UA) Distanza Eliocentr. r= 1.00928176 (UA) Long.Eclitt.Eliocen. l= 1.82229 (gradi) Lati.Eclitt.Eliocen. b= -0.01367 (gradi) Ascensione Retta AR= 62.36386 (gradi) = 4H 9M 27S Declinazione DE= 19.96857 (gradi) psi= -116.99746 (gradi) Elongaz. Ovest Angolo di Fase beta= 62.33090 (gradi) Magnitudine= 20.1 ============================= N.ro corpi selezionati= 2 =============================
Altro esempio, per gentile concessione dell'autore di AstGeoc.f90 [Fortran 90],
l'astrofilo piemontese Aldo Nicola.
Date and Time UTC
==============================
Input : Year , Month, Day (yyyy,mm,dd)
2016,10,1
Input : Maximum Geocentric distance 1.00 [UA]
(2016, revised by G. Matarazzo)
Time of Ephemeris ID 2457662.5007891473 Time of epoch JD 2457600.5000000000
--------------------------------------------------------------------------------------------------
R.A.2000 Decl.2000 Delta r Phase Mag Elong.
num body hh mm ss.ss dd pp ss.s A.U. A.U. deg V deg
--------------------------------------------------------------------------------------------------
00433 Eros 21 0 12.87 - 2 21 50.50 0.846252 1.657753 28.7 12.75 127.4 E
01863 Antinous 0 4 18.17 + 16 11 57.34 0.919574 1.905601 7.5 17.30 165.6 E
01865 Cerberus 5 55 50.21 + 19 54 23.11 0.685827 1.298868 49.6 18.44 99.0 W
01943 Anteros 1 37 11.49 + 25 41 24.22 0.783340 1.733142 15.6 17.27 152.2 W
02100 Ra-Shalom 3 7 42.15 - 25 51 57.79 0.162216 1.116143 41.8 14.03 131.9 W
03103 Eger 8 12 41.29 - 19 25 44.57 0.520085 0.909234 84.2 16.70 64.6 W
03200 Phaethon 13 58 8.90 + 82 51 22.89 0.403303 1.056688 71.0 15.24 86.6 E
03352 McAuliffe 3 1 54.05 + 6 45 50.91 0.818840 1.721235 21.0 17.58 141.9 W
03753 Cruithne 6 50 56.68 - 13 46 43.52 0.644850 1.133080 61.5 17.11 84.0 W
04179 Toutatis 18 52 11.08 - 23 21 37.69 0.633991 1.221312 54.9 16.86 93.9 E
04341 Poseidon 23 50 16.60 - 31 14 20.35 0.674437 1.598263 21.4 17.21 144.3 E
04581 Asclepius 17 8 41.34 - 19 0 22.81 0.727503 1.015213 67.8 22.43 69.9 E
04953 1990 MU 6 5 16.61 - 3 36 31.94 0.773186 1.325591 48.7 15.98 95.8 W
05381 Sekhmet 9 35 54.04 - 51 52 52.45 0.649299 0.885928 79.8 18.19 60.6 W
05653 Camarillo 4 18 22.87 + 27 1 12.91 0.934055 1.679978 30.9 18.40 120.4 W
05786 Talos 23 16 17.22 - 18 4 25.08 0.956216 1.899822 14.3 19.20 152.1 E
05836 1993 MF 3 10 7.85 + 30 24 38.95 0.545095 1.431187 30.6 15.58 133.3 W
05863 Tara 18 7 28.72 - 12 24 16.33 0.785987 1.206008 55.6 17.71 84.0 E
06239 Minos 16 7 49.78 - 23 12 17.01 0.946415 0.921460 64.8 20.50 56.4 E
06569 Ondaatje 3 18 54.45 - 34 37 49.91 0.794182 1.596334 30.9 18.24 125.1 W
07088 Ishtar 0 51 33.27 - 16 12 50.45 0.708463 1.684057 11.8 17.80 159.9 W
07341 1991 VK 0 18 39.83 + 19 14 29.00 0.497297 1.485093 10.9 16.72 163.7 E
07888 1993 UC 6 29 48.04 - 12 57 4.86 0.814011 1.281518 51.4 17.09 89.2 W
10150 1994 PN 16 20 16.84 - 56 9 18.88 0.904571 1.095014 59.2 17.50 69.9 E
10636 1998 QK56 23 10 39.21 - 10 15 51.21 0.539845 1.510947 15.5 18.01 156.2 E
11500 Tomaiyowit 5 21 31.17 + 33 9 0.18 0.516005 1.245886 50.6 19.22 105.9 W
15745 Yuliya 3 6 15.24 - 11 39 45.48 0.878535 1.757120 22.3 19.31 138.3 W
16834 1997 WU22 19 57 10.94 + 14 32 16.72 0.682906 1.407616 41.3 17.13 112.0 E
END OF PROGRAM
(Settembre 2016).Con il software asteroidale pubblicato qui appresso,
voglio rendere omaggio all'amico, collega ed astrofilo Sergio Foglia
ben conosciuto, in Italia e all'estero, per la sua poliedrica attività
di scopritore di pianetini e cultore di programmi astronomici di rara
eleganza e precisione.
Nel lontano 1996 ho avuto il privilegio di ricevere una personal release
del codice MPEPH.c per il calcolo delle effemeridi degli asteroidi tramite
lettura del Dbase degli elementi orbitali pubblicato dal Minor Planet Center (MPC).
Per avere un'idea più dettagliata basta leggere questo suo esaustivo
paper sull'argomento.
Va da sé, però, che nel corso degli anni il MPC ha cambiato lo standard del
database, aumentando considerevolmente il numero dei campi (fissi) da cui
estrarre i dati e le modalità di lettura di alcuni: stringhe invece di numeri.
Ho quindi rivisto la parte del codice da cambiare e con il permesso dell'autore
procedo alla pubblicazione della sua opera intellettuale, scaricabile da
qui; il file zippato comprende
il più recente dbase limitato ai primi 20mila asteroidi del catalogo MPC
e il nuovo listato aster16.c, per la cui compilazione è sufficiente
usare (gcc) di licenza GNU, gratuito e di facile reperibilità sia per sistemi
Windows che Linux.
Grazie di cuore, Sergio.
+-------------------------------------------+
| MINOR PLANET EPHEMERIS |
| M P E P H v. 2.1 |
| 1996, S.Foglia, Serafino Zani Observatory |
| (2016, revisited by G. Matarazzo) |
+-------------------------------------------+
Minor Planet catalogue number [1 to 20,000]: 1620
Anomalia Media M= 322.32014 (gradi)
Argomento Peri w= 276.86945 (gradi)
Longitud. Nodo O= 337.21625 (gradi)
Inclinazione i= 13.33747 (gradi)
Eccentricita' e= 0.3353926
Moto Medio n= 0.70916069 (gradi/g)
Semiasse Magg. a= 1.2453926 (UA)
Numero Catalogo : (1620)
Nome Pianetino : Geographos
Parametro H= 15.60
Parametro G= 0.15
Sigla Secolo : K [K=20; J=19]
Anni-Secolo = 2000
Sigla Anni : 16 [01-99]
Anni = 16
Sigla Mese : 7 [1-9; A=10; B=11; C=12]
Mese = 7
Sigla Giorni : V [1-9; A=10; B=11; ... U=30; V=31]
Giorno = 31
------------------------------------
Epoca Anomalia Media: anno= 2016
mese= 7
gior= 31
------------------------------------
Giorno Giuliano Epoca Anomalia Media= 2457600.50000
Ephemeris begin
year month day : 2016 09 14
Ephemeris end
year month day : 2016 09 30
Ephemeris step (days) : 2
Numero Catalogo: (1620)
Nome Pianetino: Geographos
MPEPH: Minor Planet Ephemeris
1996, S.Foglia, Serafino Zani Observatory
(2016, revisited by G. Matarazzo)
------------------------------------------------------------------------------
Date R.A.2000 Decl.2000 Delta r Phase Mag Elong.
year mo day hh mm ss.ss dd pp ss.s A.U. A.U. deg V deg
------------------------------------------------------------------------------
2016 9 14 13 10 48.65 -15 25 8.9 1.50147 0.83246 39.3 17.6 31.6E
2016 9 16 13 19 47.21 -16 22 43.0 1.48924 0.83041 39.9 17.6 32.0E
2016 9 18 13 28 56.08 -17 19 8.1 1.47726 0.82893 40.6 17.6 32.5E
2016 9 20 13 38 15.45 -18 14 13.7 1.46556 0.82803 41.2 17.6 32.9E
2016 9 22 13 47 45.46 -19 7 48.8 1.45417 0.82770 41.9 17.6 33.4E
2016 9 24 13 57 26.19 -19 59 42.4 1.44313 0.82795 42.5 17.6 33.9E
2016 9 26 14 7 17.66 -20 49 43.3 1.43247 0.82878 43.1 17.6 34.4E
2016 9 28 14 17 19.79 -21 37 39.9 1.42223 0.83018 43.7 17.6 34.9E
2016 9 30 14 27 32.44 -22 23 21.0 1.41244 0.83216 44.2 17.7 35.4E
------------------------------------------------------------------------------
(Agosto 2016). Rivisitazione del Teorema di Lambert con l'elegante algoritmo di Battin. Codice in linguaggio C e risultati sono riportati sotto.
/* Battin.c :: Teorema di Lambert (metodo di Battin) in C :: 13/15-Ago-2016
* convertito da BestBatt.bas */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.141592653589793238462643 /* in double conserva 15 cifre significative */
#define RAD 0.01745329251994329576923691 /* in double conserva 15 cifre significative */
#define KGAUSS 0.01720209895 /* in double mantiene le 11 cifre designate */
#define MU 1 /* parametro gravitazionale unitario espresso in UA^3/UT0^2 */
#define uv0 29.78469169 /* costante di conversione velocità da UC in km/s */
int main(void)
{
float r1,r2,tetagrad,tgiorni;
// double r1,r2,tetagrad,tgiorni;
double teta,t,c,s,lam,Tadim, cosf,sinf,tgf2,L,rop,m;
double x,y, yit, eta,gs0,gs1,gs2,gs3,gs4, csi;
double hden,h1,h2,U,B, fs0,fs1,fs2,fs3,fs4, ku;
double a,p,ecc, V1,V2, f,g,fpunto,gpunto,sigma1,V1L,V2L;
/*
+-------------------------------------------------------+
. Da BATTIN: An elegant Lambert algorithm - 1984 .
. (+Vaughan) (3^ parte dell'articolo- p.669-670) .
. Rimozione singolarità θ=180° e RAPIDA CONVERGENZA .
. con uso della frazione continua K(u) nonchè dell'ar- .
. ticio del calcolo di csi(x) tramite 'eta' .
+-------------------------------------------------------+
*/
printf("\tProblema di LAMBERT: calcolo della traiettoria di trasferimento\n");
printf("\tdi una sonda noti a) la posizione geometrica -> r1, r2, teta\n");
printf("\t b) il tempo di percorrenza -> t\n");
printf("\t\t\t [Moto Centrale Eliocentrico]\n");
printf("\n Raggio vettore r1 [UA], p.es. 1.0= "); scanf("%f", &r1);
printf(" Raggio vettore r2 [UA], p.es. 1.8= "); scanf("%f", &r2);
printf(" Angolo tra i due vettori teta [gradi], p.es. 109.884= "); scanf("%f", &tetagrad);
printf(" Tempo di percorrenza t [giorni], p.es. 253.86= "); scanf("%f", &tgiorni);
//
printf("\n\t r1= %.3f UA\t r2= %.3f UA\t teta= %.3f gradi\t delta-t= %.3f giorni\n\n", r1,r2,tetagrad,tgiorni);
/* ----------- Inizio calcoli ---------- */
teta= tetagrad*RAD; // angolo tra i due vettori in radianti
t= tgiorni*KGAUSS; // tempo in UC (Unità canoniche)
c= sqrt(r1*r1 + r2*r2 -2*r1*r2*cos(teta)); // c= corda del triangolo dei vettori dati
s= (r1 + r2 + c)/2 ; // semiperimetro del triangolo
//
lam= sqrt(r1*r2)/s * cos(teta/2) ; // 1^ parametro adimensionale
Tadim= sqrt(8*MU/s) * t/s; // 2^ parametro adimensionale
//
cosf= sqrt(r1*r2)*cos(teta/2)/((r1+r2)/2); sinf= sqrt(1-cosf*cosf);
tgf2= (1 - cosf)/sinf; L= tgf2*tgf2;
rop= (r1 + r2 + 2*lam*s)/4;
m= MU*t*t/(8*rop*rop*rop);
printf("\tParametri intermedi\n");
printf("\tc= %.14f UA s= %.14f UA\n\tlam=%.14f adim. Tadim=%1.14f adim. \n\n", c,s,lam,Tadim);
/* Da qui inizia il ciclo iterativo*/
x=L; // valore iniziale di x
int count=0;
goto label11;
for (;;) // ciclo infinito che è interrotto dal break
{
x= sqrt(((1 - L)/2)*((1 - L)/2) + m /(y*y)) - (1 + L)/2; /* 1^ equaz. di GAUSS */
label11:
{count++; // contatore iterazioni
eta= x/(1+sqrt(1+x))/(1+sqrt(1+x)); /* Param. ausiliario {eta} -> -1 < eta < +1 formula (52)*/
/* Calcolo della Funzione Ipergeometrica Csi = f(x,eta) */
gs3= 169.0/25.0/27.0*eta/(1 + 196.0/27.0/29.0*eta/(1 + 225.0/29.0/31.0*eta/(1 + gs4)));
gs2= 100.0/19.0/21.0*eta/(1 + 121.0/21.0/23.0*eta/(1 + 144.0/23.0/25.0*eta/(1 + gs3)));
gs1= 49.0/13.0/15.0*eta/(1 + 64.0/15.0/17.0* eta/(1 + 81.0/17.0/19.0*eta/(1 + gs2)));
gs0= 16.0/7.0/9.0*eta/(1 + 25.0/9.0/11.0* eta/(1 + 36.0/11.0/13.0*eta/(1 + gs1)));
csi= 8*(1 + sqrt(1 + x))/(3 + 1/(5 + eta + 9.0/7.0*eta/(1 + gs0))); // Formula (53)
/* Parametri ausiliari h1,h2 */
hden= (1 + 2*x + L)*(4*x + csi*(3 + x));
h1 = (L + x)*(L + x)*(1 + 3*x + csi)/hden;
h2 = m * (x - L + csi)/hden;
/* Parametri ausiliari B,U */
B= 27*h2/(4*(1+h1)*(1+h1)*(1+h1));
U= -B/2/(1 + sqrt(1 + B));
/*Calcolo della Funzione Ipergeometrica K(u) */
fs3= 928.0/3591.0*U/(1 - 1054.0/4347.0*U/(1 - 266.0/1035.0*U/(1 - fs4)));
fs2= 418.0/1755.0*U/(1 - 598.0/2295.0*U/(1 - 700.0/2907.0*U/(1 - fs3)));
fs1= 22.0/81.0*U/(1 - 208.0/891.0*U/(1 - 340.0/1287.0*U/(1 - fs2)));
fs0= 4.0/27.0*U/(1 - 8.0/27.0*U/(1 - 2.0/9.0*U/ (1 - fs1)));
ku= 1.0/3.0/(1 - fs0); //Formula (58)
/* 2^ Equaz. di Gauss -> Radice positiva y */
yit= (1 + h1)/3*(2 + sqrt(1 + B)/(1 - 2*U*ku*ku));
printf("\t %d) \t yit = %.8f\t Residuo= %.2e \n", count,yit,yit-y);
} // chiude da count
if (fabs(yit-y) < 1E-6)
break;
else
y=yit;
} // fine ciclo infinito di for
printf("\n\t\t y-fin = %.8f\t x-fin= %.8f \n", yit,x);
/* Calcolo dei PARAMETRI ORBITALI {a,p,e} */
a= m*s*(1 + lam)*(1 + lam)/(8 * x * yit*yit);
p= 2*r1*r2*(y*(1 + x)*sin(teta/2))*(y*(1 + x)*sin(teta/2))/((1 + lam)*(1 + lam))/m/s;
ecc= sqrt(1 - p/a);
/* Calcolo delle VELOCITA' {V1, V2} - Metodo CLASSICO */
V1= sqrt(2* MU/r1 - MU/a); V2= sqrt(2* MU/r2 - MU/a);
/* Calcolo delle VELOCITA' {V1, V2} tramite i parametri lagrangiani {f,g,f',g'}
per verificare se i parametri lagrangiani di sotto sono CORRETTI */
f= 1 - r2/p*2*sin(teta/2)*sin(teta/2);
g= r1*r2*sin(teta)/sqrt(MU*p);
fpunto= (2/p*sin(teta/2)*sin(teta/2) - 1/r1 - 1/r2)*tan(teta/2)*sqrt(MU/p);
gpunto= 1 - r1/p*2*sin(teta/2)*sin(teta/2);
//
V1L = sqrt((r2*r2 - 2*f*r1*r2*cos(teta) + f*f * r1*r1)/(g*g));
/* Vettore V2 = fpunto * r1 + gpunto * V1 */
sigma1= (r1*r2*cos(teta) - f*r1*r1)/g;
V2L= sqrt((fpunto*r1)*(fpunto*r1) + (gpunto*V1L)*(gpunto*V1L) + 2*fpunto*gpunto*sigma1);
/* OUTPUT finali */
printf("\n\n\t RISULTATI della CONICA (parametri di forma e velocita')\n");
printf("\t -------------------------------------------------------\n");
printf("\t\t a= semiasse magg.= %.10f UA\n\t\t p= param.focale= %.10f UA\n\t\t e= eccentricita'= %.10f\n", a,p,ecc);
printf("\n\t\tCalcolo delle velocita' con il metodo classico\n");
printf("\t\t V1= velocita' in P1= %.5f km/s\n\t\t V2= velocita' in P2= %.5f km/s\n\n", V1*uv0,V2*uv0);
printf("\t\tCalcolo delle velocita' con i parametri lagrangiani\n");
printf("\t\t V1= velocita' in P1= %.5f km/s\n\t\t V2= velocita' in P2= %.5f km/s\n", V1L*uv0,V2L*uv0);
printf("\t ------------------- Fine Calcoli ----------------------\n\n");
return (EXIT_SUCCESS);
}
Problema di LAMBERT: calcolo della traiettoria di trasferimento
di una sonda noti a) la posizione geometrica -> r1, r2, teta
b) il tempo di percorrenza -> t
[Moto Centrale Eliocentrico]
Raggio vettore r1 [UA], p.es. 1.0= 1
Raggio vettore r2 [UA], p.es. 1.8= 1.8
Angolo tra i due vettori teta [gradi], p.es. 109.884= 109.884
Tempo di percorrenza t [giorni], p.es. 253.86= 253.86
r1= 1.000 UA r2= 1.800 UA teta= 109.884 gradi delta-t= 253.860 giorni
Parametri intermedi
c= 2.33761010582490 UA s= 2.56880502907059 UA
lam=0.30000160296500 adim. Tadim=3.00001894683454 adim.
1) yit = 1.25093289 Residuo= 1.25e+000
2) yit = 1.24381067 Residuo= -7.12e-003
3) yit = 1.24380628 Residuo= -4.39e-006
4) yit = 1.24380628 Residuo= -1.68e-012
y-fin = 1.24380628 x-fin= 0.50885139
RISULTATI della CONICA (parametri di forma e velocita')
-------------------------------------------------------
a= semiasse magg.= 1.2853374112 UA
p= param.focale= 1.0495667071 UA
e= eccentricita'= 0.4282884337
Calcolo delle velocita' con il metodo classico
V1= velocita' in P1= 32.92514 km/s
V2= velocita' in P2= 17.19032 km/s
Calcolo delle velocita' con i parametri lagrangiani
V1= velocita' in P1= 32.92514 km/s
V2= velocita' in P2= 17.19032 km/s
------------------- Fine Calcoli ----------------------
(Luglio 2016). Raccolta dei miei lavori (works) anno per anno.
1997
File PDF da aprire o scaricare
Aggiornamento del codice di MOID in MatLab, prelevabile da questo file zippato e contenente il relativo DataBase. Di sotto sono riportati listato e output di un esempio applicativo.
% ---------- moid_new.m --------- 21 Luglio 2016
% Nota: la creazione della funzione m_iter.m è opera dell'amico
% e astrofilo Aldo Nicola, che ringrazio.
% ---------------------------------------------------------------------------
%
clc; clear all;
rad=pi/180;
format long
%
disp ('----------------------------------------------');
disp (' Calcolo del MOID ');
disp (' DataBase di 20,000 pianetini ');
disp ('----------------------------------------------');
disp ('');
%1^ parte: viene selezionata la riga del dBase, scartando la colonna con %*f
i=input('Numero di catalogo dell''asteroide [da 1 a 20,000]: ');
[perielio nodo inclinazione eccentricita semiassemaggiore nome] = ...
textread('DB20k.txt', '%f %f %f %f %*f %f %s',i);
%
nome(i)
fprintf ('\n omega=%10.5f Omega=%10.5f i=%9.5f e=%10.7f a=%10.7f\n', ...
perielio(i), nodo(i), inclinazione(i), eccentricita(i), semiassemaggiore(i));
disp('');
% ---------------- Fine lettura dei parametri del Dbase -----------
a= semiassemaggiore(i); e= eccentricita(i);
% ------- angoli in radianti --------
lnodo=nodo(i)*rad; om=perielio(i)*rad; incl=inclinazione(i)*rad;
%
% Terra (Fonte: Astronomical Algorithms di J.Meeus)
et = .01670862; lp = 102.937348*rad; at = 1.000001018;
%'****************************************************************************
p = a * (1 - e ^ 2); pt = at * (1 - et ^ 2);
%'****************************************************************************
% Angoli di EULERO {i, Omega, omega} --> Per Calcolo MATRICE di Trasformazione
R11 = cos(lnodo) * cos(om) - sin(lnodo) * cos(incl) * sin(om);
R12 = -cos(lnodo) * sin(om) - sin(lnodo) * cos(incl) * cos(om);
R21 = sin(lnodo) * cos(om) + cos(lnodo) * cos(incl) * sin(om);
R22 = -sin(lnodo) * sin(om) + cos(lnodo) * cos(incl) * cos(om);
R31 = sin(incl) * sin(om);
R32 = sin(incl) * cos(om);
%
% a) scelta della coppia dei valori iniziali (te,el) vicino al nodo ascendente
te = -om; el = lnodo;
% ----------------------------------------------
[x,y,z,xt,yt] = m_iter (R11,R21,R31,R12,R22,R32,pt,e,te,el,et,lp,p);
% ----------------------------------------------
moid1=sqrt((x-xt)^2+(y-yt)^2+z^2);
fprintf('\n d1= minima distanza tra le orbite nello spazio [nodo Asc]= %8.6f UA\n\n',moid1);
%
% b) scelta della coppia dei valori iniziali (te,el) vicino al nodo discendente
te = -om+pi; el = lnodo+pi;
% ----------------------------------------------
[x,y,z,xt,yt] = m_iter (R11,R21,R31,R12,R22,R32,pt,e,te,el,et,lp,p);
% ----------------------------------------------
moid2=sqrt((x-xt)^2+(y-yt)^2+z^2);
fprintf('\n d2= minima distanza tra le orbite nello spazio [nodo Disc]= %8.6f UA\n\n',moid2);
fprintf(' ----------------------------------------------------- \n');
if moid1 > moid2
moid = moid2;
else
moid = moid1;
end
%
fprintf(' MOID= Distanza minima Terra-Pianetino= %8.6f UA\n',moid);
fprintf(' ----------------------------------------------------- \n');
% ---- EOF: moid_new.m
function [x,y,z,xt,yt] = m_iter (R11,R21,R31,R12,R22,R32,pt,e,te,el,et,lp,p)
%
% --- file: m_iter.m da copiare nella stassa cartella di moid_new.m
%
ddte=1; ddel=1; % inizializzazione variabili di controllo
while (abs(ddte) > 1e-8) && (abs(ddel) > 1e-8)
%
x = (R11 * cos(te) + R12 * sin(te)) * p/(1+e*cos(te));
y = (R21 * cos(te) + R22 * sin(te)) * p/(1+e*cos(te));
z = (R31 * cos(te) + R32 * sin(te)) * p/(1+e*cos(te));
xt = pt * cos(el)/(1+et*cos(el-lp));
yt = pt * sin(el)/(1+et*cos(el-lp));
%----------------------------------------------------------------------------
Atil = -R11 * sin(te) + R12 * (cos(te) + e);
Btil = -R21 * sin(te) + R22 * (cos(te) + e);
Ctil = -R31 * sin(te) + R32 * (cos(te) + e);
Dtil = et * (-sin(el) * cos(el - lp) + cos(el) * sin(el - lp)) - sin(el);
Etil = et * (cos(el) * cos(el - lp) + sin(el) * sin(el - lp)) + cos(el);
%
A = Atil*p/(1+e*cos(te))^2;
B = Btil*p/(1+e*cos(te))^2;
C = Ctil*p/(1+e*cos(te))^2;
D = Dtil*pt/(1+et*cos(el-lp))^2;
E = Etil*pt/(1+et*cos(el-lp))^2;
%
Aast = -R11 * cos(te) - R12 * sin(te);
Bast = -R21 * cos(te) - R22 * sin(te);
Cast = -R31 * cos(te) - R32 * sin(te);
Dast = -cos(el);
East = -sin(el);
%-----------------------------------------------------------------------
F1 = Atil*(x-xt)+Btil*(y-yt)+Ctil*z;
F2 = Dtil*(x-xt)+Etil*(y-yt);
%
F11 = Aast*(x-xt) + Atil*A + Bast*(y-yt) + Btil*B + Cast*z + Ctil*C;
F12 = -Atil*D - Btil*E;
F21 = Dtil*A + Etil*B;
F22 = Dast*(x-xt) - Dtil*D + East*(y-yt) - Etil*E;
%
W0 = F11*F22 - F21*F12;
W1 = F1*F22 - F2*F12;
W2 = F11*F2 - F21*F1;
%
ddte= W1/W0; ddel= W2/W0;
te1= te - ddte; el1 = el - ddel; %Iterazione di Newton
%
fprintf ('\n te_i = %12.10f el_i = %12.10f', te1,el1);
te = te1; el=el1;
end
--------------------------------------------------------------------------
PHA's - Potential Hazardous Asteroids nell'intervallo [1-20,000]
--------------------------------------------------------------------------
(1566) Icarus (4034) Vishnu (6491) 1991 OA
(1620) Geographos (4179) Toutatis (7335) 1989 JA
(1862) Apollo (4183) Cuno (7341) 1991 VK
(1981) Midas (4450) Pan (7482) 1994 PC1
(2101) Adonis (4486) Mithra (7753) 1988 XB
(2102) Tantalus (4581) Asclepius (7822) 1991 CS
(2135) Aristaeus (4660) Nereus (8014) 1990 MF
(2201) Oljato (4769) Castalia (8566) 1996 EN
(2340) Hathor (4953) 1990 MU (9856) 1991 EE
(3122) Florence (5011) Ptah (10115) 1992 SK
(3200) Phaethon (5189) 1990 UQ (11500) Tomaiyowit
(3361) Orpheus (5604) 1992 FE (12538) 1998 OH
(3362) Khufu (5693) 1993 EA (12923) Zephyr
(3671) Dionysus (6037) 1988 EG (13651) 1997 BR
(3757) Anagolay (6239) Minos (14827) Hypnos
(4015) Wilson-Harrington (6489) Golevka (16960) 1998 QS52
--------------------------------------------------------------------------
----------------------------------------------
Calcolo del MOID
DataBase di 20,000 pianetini
----------------------------------------------
Numero di catalogo dell'asteroide [da 1 a 20,000]: 4179
ans =
'(4179)-Toutatis'
omega= 278.73541 Omega= 124.37315 i= 0.44711 e= 0.6294515 a= 2.5342052
te_i = -5.3537474518 el_i = 1.6797376225
te_i = -5.6123064511 el_i = 1.4225598679
te_i = -5.7380636521 el_i = 1.2973889389
te_i = -5.7782101638 el_i = 1.2573642509
te_i = -5.7826308079 el_i = 1.2529544594
te_i = -5.7826834982 el_i = 1.2529018907
te_i = -5.7826835057 el_i = 1.2529018833
d1= minima distanza tra le orbite nello spazio [nodo Asc]= 0.006111 UA
te_i = -1.0867007754 el_i = 5.9450770007
te_i = -0.7960267800 el_i = 6.2383106446
te_i = -0.6426328717 el_i = 6.3926840315
te_i = -0.5825931891 el_i = 6.4529380460
te_i = -0.5722515375 el_i = 6.4633071040
te_i = -0.5719498614 el_i = 6.4636095082
te_i = -0.5719496075 el_i = 6.4636097626
te_i = -0.5719496075 el_i = 6.4636097626
d2= minima distanza tra le orbite nello spazio [nodo Disc]= 0.007131 UA
-----------------------------------------------------
MOID= Distanza minima Terra-Pianetino= 0.006111 UA
-----------------------------------------------------
1997
File PDF da aprire o scaricare
Aggiornamento del codice delle Variabili Universali in C. Di sotto sono riportati listato (variabuniv.c) e output di un'orbita parabolica.
/* variabuniv.c, Rif. 1997b --- 07.Ago.2016 ---
* ' Dati iniziali (Elementi UNIVERSALI) nell'ordine: al, q, i, Om, om, tau
* */
#include ⟨stdio.h⟩
#include ⟨stdlib.h⟩
#include ⟨math.h⟩
#define PI 3.141592653589793238462643 /* in double conserva 15 cifre significative */
#define RAD 0.01745329251994329576923691 /* in double conserva 15 cifre significative */
#define KGAUSS 0.01720209895 /* in double mantiene le 11 cifre designate */
#define MU 1
int main(void)
{
float AL, Q, EI, BOM, OM, TAU;
double TOL = 1E-10;
double EIrad,BOMrad,OMrad, PX,PY,PZ, QX,QY,QZ, VP,XP,YP,ZP, VXP,VYP,VZP;
double PSIA,PSI, B,S3PARZ,S3, S2PARZ,S2, S1,S0;
double R,F,G, F1,G1, X,Y,Z, VX,VY,VZ, V;
/* Esempio orbita ELLITTICA
AL = -1.321604534; Q = .455635165; EI = 6.87631;
BOM = 63.07572; OM = 205.6752; TAU = 1.457528167;
*/
// Esempio orbita PARABOLICA
AL = 0; Q = 0.22432; EI = 122.639;
BOM = 188.943; OM = 131.202; TAU = -0.632503976;
/* Esempio orbita IPERBOLICA
AL = 4.98736E-4; Q = 0.555404; EI = 72.5488;
BOM = 237.8971; OM = 276.7690; TAU = 1.157986814;
*/
printf("\t\t ---------- I 6 Dati Orbitali Universali ------------\n");
printf(" alfa=(-mu/a)= %.8f q=dist.peri.= %.8f UA tau= t-Tp= %.8f TU\n",AL,Q,TAU);
printf(" Om= Long.Nodo= %.8f gr. w= Arg.peri.= %.8f gr. incl= %.8f gr.\n",BOM,OM,EI);
printf("\t\t ----------------------------------------------------\n\n");
EIrad= EI*RAD; BOMrad= BOM*RAD; OMrad= OM*RAD;
PX= cos(OMrad)*cos(BOMrad) - sin(OMrad)*cos(EIrad)*sin(BOMrad);
PY= cos(OMrad)*sin(BOMrad) + sin(OMrad)*cos(EIrad)*cos(BOMrad);
PZ= sin(OMrad)*sin(EIrad);
QX= -sin(OMrad)*cos(BOMrad) - cos(OMrad)*cos(EIrad)*sin(BOMrad);
QY= -sin(OMrad)*sin(BOMrad) + cos(OMrad)*cos(EIrad)*cos(BOMrad);
QZ= cos(OMrad)*sin(EIrad);
VP= sqrt(2*MU/Q + AL); /*Velocità al Perielio*/
XP= Q*PX; YP= Q*PY; ZP= Q*PZ;
VXP= VP*QX; VYP= VP*QY; VZP= VP*QZ;
// Risoluzione equazione di Keplero espressa con le V.U.
PSIA= TAU/Q; /* PSIA= valore iniziale dell'iterazione*/
for (;;) // ciclo infinito che è interrotto dal break
{
B= AL*PSIA*PSIA;
S3PARZ= ((((B/342 + 1)*B/272+1)*B/210 + 1)*B/156 + 1);
S3= PSIA*PSIA*PSIA/6*((((S3PARZ*B/110 + 1)*B/72 + 1)*B/42 + 1)*B/20 + 1);
S2PARZ= ((((B/306 + 1)*B/240 + 1)*B/182 + 1)*B/132 + 1);
S2= PSIA*PSIA/2*((((S2PARZ*B/90 + 1)*B/56 + 1)*B/30 + 1)*B/12 + 1);
S1 = PSIA + AL*S3;
S0 = 1 + AL*S2;
PSI= PSIA - (Q*S1 + MU*S3 - TAU)/(Q*S0 + MU*S2);
printf("\t PSI_iterato= %.10f\n",PSI);
if (fabs(PSI - PSIA) ⟨ TOL)
break;
else
PSIA= PSI;
} // fine ciclo infinito di for
printf("\n\t PSI_finale (PSIA)= %.10f\n",PSIA);
printf(" \t PSI_finale ( PSI)= %.10f\n\n",PSI);
/* Calcolo dei Vettori Posizione e Velocita'*/
R= Q*S0 + MU*S2; F= 1 - MU*S2/Q; G= TAU - MU*S3;
F1= -MU*S1/(R*Q); G1= 1 - MU*S2/R;
X= F*XP + G*VXP;
Y= F*YP + G*VYP;
Z= F*ZP + G*VZP;
VX= F1*XP + G1*VXP;
VY= F1*YP + G1*VYP;
VZ= F1*ZP + G1*VZP;
V= sqrt(VX*VX + VY*VY + VZ*VZ);
printf("\t X= %.8f Y= %.8f Z= %.8f --> R= %.8f UA\n",X,Y,Z,R);
printf("\t VX= %.8f VY= %.8f VZ= %.8f --> V= %.8f UA/UT\n",VX,VY,VZ,V);
return (EXIT_SUCCESS);
}
---------- I 6 Dati Orbitali Universali ------------
alfa=(-mu/a)= 0.00000000 q=dist.peri.= 0.22431999 UA tau= t-Tp= -0.63250399 TU
Om= Long.Nodo= 188.94299316 gr. w= Arg.peri.= 131.20199585 gr. incl= 122.63899994 gr.
----------------------------------------------------
PSI_iterato= -1.9299712990
PSI_iterato= -1.4514440116
PSI_iterato= -1.2927892425
PSI_iterato= -1.2761833321
PSI_iterato= -1.2760124511
PSI_iterato= -1.2760124331
PSI_iterato= -1.2760124331
PSI_finale (PSIA)= -1.2760124331
PSI_finale ( PSI)= -1.2760124331
X= -1.02901223 Y= -0.09682577 Z= 0.10041273 --> R= 1.03842386 UA
VX= 1.23710597 VY= 0.46747707 VZ= 0.42074911 --> V= 1.38780251 UA/UT
2002
File PDF da aprire o scaricare
2002
File PDF da aprire o scaricare
2003
File PDF da aprire o scaricare
2003
File PDF da aprire o scaricare
2003
File PDF da aprire o scaricare
2004
File PDF da aprire o scaricare
2004
File PDF da aprire o scaricare
2004
File PDF da aprire o scaricare
File PDF da aprire o scaricare
2004
2004
File PDF da aprire o scaricare
2004
File PDF da aprire o scaricare
2005
File PDF da aprire o scaricare
2005
File PDF da aprire o scaricare
2005
File PDF da aprire o scaricare
Aggiornamento del codice degli elemnti orbitali noti 2 Vettori e l'angolo con la tangente in C. Di sotto sono riportati listato (orbit-r1-r2.c) e output.
/* orbit-r1-r2.c :: Riferimento paper 2005c.pdf ---- 02-Ago-2016 --- */
#include ⟨stdio.h⟩
#include ⟨stdlib.h⟩
#include ⟨math.h⟩
#define PI 3.141592653589793238462643 /* in double conserva 15 cifre significative */
#define RAD 0.01745329251994329576923691 /* in double conserva 15 cifre significative */
#define KGAUSS 0.01720209895 /* in double mantiene le 11 cifre designate */
#define MU 1
int main(void)
{
/* Dati i vettori (r1,r2) con coordinate in float */
float x1 = -0.106418; float y1 = 0.137154; float z1 = 1.637343;
float x2 = -2.60002887; float y2 = 1.62023766; float z2 = 2.21048897;
float be = 63.54333316; /* angolo beta in gradi */
/* Inizio Calcoli*/
double r1= sqrt(x1*x1 + y1*y1 + z1*z1);
double r2= sqrt(x2*x2 + y2*y2 + z2*z2);
double cosalfa = (x1*x2 + y1*y2 + z1*z2)/(r1*r2);
double alfa = acos(cosalfa);
/* Costanti */
double K2 = cosalfa; double K1 = sqrt(1-K2*K2);
double B = r1 * sin(be*RAD); double K0 = 1/(B*B);
double K4 = r1/r2 + K1*sqrt(K0*r1*r1 - 1);
/* Calcolo di x=e*COS(te1), y=e^2, z=p */
double x = (1 - K4)/(K4 - K2);
double y = K0*r1*r1 * (1+x)*(1+x) - 2*x - 1;
double z = r1*(1+x);
double ecc = sqrt(y);
double p = z;
double a = p/(1-ecc*ecc);
/* Calcolo Anomalie (te1,te2) */
double te1 = acos(x/ecc);
double te2 = te1 + alfa;
/* Calcolo dei Versori Perifocali P=(Px,Py,Pz) e Q=(Qx,Qy,Qz) */
double a1= r1*cos(te1); double b1= r1*sin(te1);
double a2= r2*cos(te2); double b2= r2*sin(te2);
double Det= a1*b2 - b1*a2;
double Px= (x1*b2 - x2*b1)/ Det;
double Py= (y1*b2 - y2*b1)/ Det;
double Pz= (z1*b2 - z2*b1)/ Det;
double Qx= (a1*x2 - a2*x1)/ Det;
double Qy= (a1*y2 - a2*y1)/ Det;
double Qz= (a1*z2 - a2*z1)/ Det;
/* Calcolo delle Velocità (V1,V2) */
double V1= sqrt(MU*(1 + 2*ecc*cos(te1) + ecc*ecc)/p);
double V2= sqrt(MU*(1 + 2*ecc*cos(te2) + ecc*ecc)/p);
/* Calcolo degli elementi orbitali (i,Omega,w) */
double om= atan2(Pz,Qz);
if (om ⟨ 0)
om = om + 2*PI;
double num1= Py*cos(om) - Qy*sin(om);
double den1= Px*cos(om) - Qx*sin(om);
double nodo= atan2(num1,den1);
if (nodo ⟨ 0)
nodo = nodo + 2*PI;
double num2= Pz/sin(om);
double den2= -(Px*sin(om) + Qx*cos(om))/sin(nodo);
double incl= atan2(num2,den2);
/* OUTPUTS */
printf("\t Modulo Vettore r1= %.8f UA\n",r1);
printf("\t Modulo Vettore r2= %.8f UA\n\n",r2);
printf("\t Angolo tra i due Vettori= %.8f gradi \n\n",alfa/RAD);
printf("\t a= %.8f UA ecc= %.8f p= %.8f UA\n\n",a,ecc,p);
if (a ⟨ 0)
printf("\t La Traiettoria e' IPERBOLICA\n\n");
else
printf("\t La Traiettoria e' ELLITTICA\n\n");
printf("\t Anomalia-1= %.8f gradi\n",te1/RAD);
printf("\t Anomalia-2= %.8f gradi\n\n",te2/RAD);
printf("\t Arg.Peri(w)= %.8f gradi\n",om/RAD);
printf("\t Long.Nodo(O)= %.8f gradi\n",nodo/RAD);
printf("\t Inclinaz.(i)= %.8f gradi\n\n",incl/RAD);
printf("\t\t Fine Calcoli \n");
return (EXIT_SUCCESS);
}
Modulo Vettore r1= 1.64652005 UA
Modulo Vettore r2= 3.77777468 UA
Angolo tra i due Vettori= 48.54151166 gradi
a= -1.88461051 UA ecc= 1.73559597 p= 3.79238923 UA
La Traiettoria e' IPERBOLICA
Anomalia-1= 41.33077865 gradi
Anomalia-2= 89.87229031 gradi
Arg.Peri(w)= 54.28322750 gradi
Long.Nodo(O)= 329.70534118 gradi
Inclinaz.(i)= 87.73564140 gradi
2005
File PDF da aprire o scaricare
2005
File PDF da aprire o scaricare
2008
File PDF da aprire o scaricare
2012
File PDF da aprire o scaricare
Aggiornamento del codice delle Variabili Universali in C. Di sotto sono riportati listato (variabuniv.c) e output di un'orbita parabolica.
/* variabuniv.c, Rif. 1997b --- 07.Ago.2016 ---
* ' Dati iniziali (Elementi UNIVERSALI) nell'ordine: al, q, i, Om, om, tau
* */
#include ⟨stdio.h⟩
#include ⟨stdlib.h⟩
#include ⟨math.h⟩
#define PI 3.141592653589793238462643 /* in double conserva 15 cifre significative */
#define RAD 0.01745329251994329576923691 /* in double conserva 15 cifre significative */
#define KGAUSS 0.01720209895 /* in double mantiene le 11 cifre designate */
#define MU 1
int main(void)
{
float AL, Q, EI, BOM, OM, TAU;
double TOL = 1E-10;
double EIrad,BOMrad,OMrad, PX,PY,PZ, QX,QY,QZ, VP,XP,YP,ZP, VXP,VYP,VZP;
double PSIA,PSI, B,S3PARZ,S3, S2PARZ,S2, S1,S0;
double R,F,G, F1,G1, X,Y,Z, VX,VY,VZ, V;
/* Esempio orbita ELLITTICA
AL = -1.321604534; Q = .455635165; EI = 6.87631;
BOM = 63.07572; OM = 205.6752; TAU = 1.457528167;
*/
// Esempio orbita PARABOLICA
AL = 0; Q = 0.22432; EI = 122.639;
BOM = 188.943; OM = 131.202; TAU = -0.632503976;
/* Esempio orbita IPERBOLICA
AL = 4.98736E-4; Q = 0.555404; EI = 72.5488;
BOM = 237.8971; OM = 276.7690; TAU = 1.157986814;
*/
printf("\t\t ---------- I 6 Dati Orbitali Universali ------------\n");
printf(" alfa=(-mu/a)= %.8f q=dist.peri.= %.8f UA tau= t-Tp= %.8f TU\n",AL,Q,TAU);
printf(" Om= Long.Nodo= %.8f gr. w= Arg.peri.= %.8f gr. incl= %.8f gr.\n",BOM,OM,EI);
printf("\t\t ----------------------------------------------------\n\n");
EIrad= EI*RAD; BOMrad= BOM*RAD; OMrad= OM*RAD;
PX= cos(OMrad)*cos(BOMrad) - sin(OMrad)*cos(EIrad)*sin(BOMrad);
PY= cos(OMrad)*sin(BOMrad) + sin(OMrad)*cos(EIrad)*cos(BOMrad);
PZ= sin(OMrad)*sin(EIrad);
QX= -sin(OMrad)*cos(BOMrad) - cos(OMrad)*cos(EIrad)*sin(BOMrad);
QY= -sin(OMrad)*sin(BOMrad) + cos(OMrad)*cos(EIrad)*cos(BOMrad);
QZ= cos(OMrad)*sin(EIrad);
VP= sqrt(2*MU/Q + AL); /*Velocità al Perielio*/
XP= Q*PX; YP= Q*PY; ZP= Q*PZ;
VXP= VP*QX; VYP= VP*QY; VZP= VP*QZ;
// Risoluzione equazione di Keplero espressa con le V.U.
PSIA= TAU/Q; /* PSIA= valore iniziale dell'iterazione*/
for (;;) // ciclo infinito che è interrotto dal break
{
B= AL*PSIA*PSIA;
S3PARZ= ((((B/342 + 1)*B/272+1)*B/210 + 1)*B/156 + 1);
S3= PSIA*PSIA*PSIA/6*((((S3PARZ*B/110 + 1)*B/72 + 1)*B/42 + 1)*B/20 + 1);
S2PARZ= ((((B/306 + 1)*B/240 + 1)*B/182 + 1)*B/132 + 1);
S2= PSIA*PSIA/2*((((S2PARZ*B/90 + 1)*B/56 + 1)*B/30 + 1)*B/12 + 1);
S1 = PSIA + AL*S3;
S0 = 1 + AL*S2;
PSI= PSIA - (Q*S1 + MU*S3 - TAU)/(Q*S0 + MU*S2);
printf("\t PSI_iterato= %.10f\n",PSI);
if (fabs(PSI - PSIA) ⟨ TOL)
break;
else
PSIA= PSI;
} // fine ciclo infinito di for
printf("\n\t PSI_finale (PSIA)= %.10f\n",PSIA);
printf(" \t PSI_finale ( PSI)= %.10f\n\n",PSI);
/* Calcolo dei Vettori Posizione e Velocita'*/
R= Q*S0 + MU*S2; F= 1 - MU*S2/Q; G= TAU - MU*S3;
F1= -MU*S1/(R*Q); G1= 1 - MU*S2/R;
X= F*XP + G*VXP;
Y= F*YP + G*VYP;
Z= F*ZP + G*VZP;
VX= F1*XP + G1*VXP;
VY= F1*YP + G1*VYP;
VZ= F1*ZP + G1*VZP;
V= sqrt(VX*VX + VY*VY + VZ*VZ);
printf("\t X= %.8f Y= %.8f Z= %.8f --> R= %.8f UA\n",X,Y,Z,R);
printf("\t VX= %.8f VY= %.8f VZ= %.8f --> V= %.8f UA/UT\n",VX,VY,VZ,V);
return (EXIT_SUCCESS);
}
---------- I 6 Dati Orbitali Universali ------------
alfa=(-mu/a)= 0.00000000 q=dist.peri.= 0.22431999 UA tau= t-Tp= -0.63250399 TU
Om= Long.Nodo= 188.94299316 gr. w= Arg.peri.= 131.20199585 gr. incl= 122.63899994 gr.
----------------------------------------------------
PSI_iterato= -1.9299712990
PSI_iterato= -1.4514440116
PSI_iterato= -1.2927892425
PSI_iterato= -1.2761833321
PSI_iterato= -1.2760124511
PSI_iterato= -1.2760124331
PSI_iterato= -1.2760124331
PSI_finale (PSIA)= -1.2760124331
PSI_finale ( PSI)= -1.2760124331
X= -1.02901223 Y= -0.09682577 Z= 0.10041273 --> R= 1.03842386 UA
VX= 1.23710597 VY= 0.46747707 VZ= 0.42074911 --> V= 1.38780251 UA/UT
(06-Lug-13). L'uploading dei miei lavori finisce qui
