Le prove fatte finora mostrano che l'implementazione e l'uso di SPICE versione C per Linux è abbastanza semplice. Ci mancano le funzioni di calcolo di posizione che andiamo a scrivere subito.
Iniziamo dalla funzione di calcolo di posizione spkezr_c, vediamo la signature.
void spkezr_c ( ConstSpiceChar *targ,
SpiceDouble et,
ConstSpiceChar *ref,
ConstSpiceChar *abcorr,
ConstSpiceChar *obs,
SpiceDouble starg[6],
SpiceDouble *lt )
Abbiamo quindi bisogno, in input, di una stringa, di un double, di tre stringhe, di un'array di 6 double e di un puntatore a double
la pagina di documentazione ci chiarisce quali sono gli input e quali gli output:
Brief_I/O Variable I/O Description -------- --- -------------------------------------------------- targ I Target body name. et I Observer epoch. ref I Reference frame of output state vector. abcorr I Aberration correction flag. obs I Observing body name. starg O State of target. lt O One way light time between observer and target.
Dobbiamo quindi fornire il nome dell'oggetto target (un pianeta, il sole, la luna o altro, usando la dicitura adatta, che vediamo fra un po'. quindi un'epoch, cioè un istante temporale riferito all'osservatore, il nome del frame di riferimento, l'eventuale correzione per l'aberrazione e infine il punto da cui facciamo l'osservazione (per un astrologo è, di regola, il centro della Terra).
Nel nostro file jpltoolkit.py aggiungeremo le seguenti righe:
... _spkezr_c = ct.shared_lib.spkezr_c ... def spkezr_c(targ, ref, timestring, abcorr, obs): _spkezr_c.argtypes = [ct.c_char_p, ct.c_double, ct.c_char_p, ct.c_char_p, ct.c_char_p, ct.c_double * 6, ct.POINTER(ct.c_double)] et = str2et_c(timestring, ct.c_double(0)) _starg = ct.c_double * 6; starg = _starg(0,0,0,0,0,0) lt = ct.c_double(0) _spkezr_c(targ, et, ref, abcorr, obs, starg, lt) return starg
Notiamo che gli output del linguaggio C sono a loro volta degli input (inizializzati a zero, in ogni caso). Anche questa funzione è void, non restituisce un valore. Nella funzione wrapper di Python decidiamo cosa e come restituiamo al chiamante.
Diamo un'occhiata ai parametri della funzione python: dobbiamo fornire il target, che puo' essere "SUN" "MOON" o altro pianeta, il ref, che in prima istanza sarà "J2000", il sistema di riferimento equatoriale di default, che ci fornirà ascensione retta e declinazione, una timestring come abbiamo già visto formattata secondo la ISO 8061, "NONE" per abcorr, la variabile relativa alla correzione per l'aberrazione (diamo per scontato che non ci interessano i fenomeni apparenti e la visione diretta degli oggetti con il telescopio) e il punto di osservazione ("EARTH"). La definizione della funzione _spkezr_c conterrà anche le variabili di output, starg[6] e lt (quest'ultima verrà ignorata).
Dopo la definiziaone degli argtypes, creiamo le variabili di input vere e e proprie. La timestring che forniamo come parametro è consumata dalla funzione già vista str2et_c. Le variabili stringa sono già definite. La variabile di output starg, la variabile che contiene 6 risultati, 3 riferiti alla posizione secondo le coordinate equatoriali rettangolari e tre riferiti al vettore velocità nelle tre componenti spaziali, è definita come array di 6 elementi e inizializzata con valori fittizi pari a 0. Essendo una variabile ctypes va definita come ct.c_double * 6 e riempita con una variabile di appoggio _starg (array di 6 zeri). La chiamata interna alla funzione _spkezr_c a questo punto è semplice. Il vettore di 6 elementi starg, dopo l'uscita da C, conterrà i valori che ci interessano, e saranno questi che con il return restituiremo al chiamante. Facciamo una prova: SATURN_BARYCENTER (NAIF ID usato dal kernel SPK) come target, EARTH come observer, J2000 come frame, "2017-08-10T18:53:22" come timestring (la T maiuscola è prevista dal formato ISO 8061), "NONE" per l'aberrazione.
Il codice e il risultato della chiamata a funzione sono i seguenti:
timestring = "2000-01-01T12:00:00" et = jpl.ct.c_double(0) seconds = jpl.str2et_c(timestring, et) print "(UTC 2000_01_01 12:00)", seconds, "seconds from J2000\n" timestring = "2017-08-10T18:53:22" ref = "J2000" target = "SATURN_BARYCENTER" observer = "EARTH" abcorr = "NONE" xyz = jpl.spkezr_c(target, ref, timestring, abcorr, observer) elements = ["x :", "y :", "z :", "vx :", "vy :", "vz :"] for i in range(6): print elements[i], xyz[i] recvector = (xyz[0], xyz[1], xyz[2]) r, RA, dec = jpl.recrad_c(recvector) print r, ddd2dms(RA * jpl.dpr_c()/15.), ddd2dms(dec * jpl.dpr_c()) >>> x : -218912917.201 y : -1296275066.93 z : -528970291.665 vx : -10.2998651768 vy : -20.7463383194 vz : -9.37634837706
I risultati sono espressi in km (per le posizioni, essendo distanze dalla Terra) e km al secondo (per le componenti vettoriali della velocità). Ora dobbiamo verificare il risultato. La JPL/NASA mette a disposizione un' interfaccia web per il calcolo online, usando proprio il toolkit che stiamo implementando. Nel form web non trovo, però, la possibilità di impostare le coordinate rettangolari, per cui devo precedere oltre e implementare la conversione in ascensione retta e declinazione (coordinate sferiche equatoriali)
Aggiungiamo la funzione recrad_c che effettua tale conversione. La signature è la seguente:
void recrad_c ( ConstSpiceDouble rectan[3],
SpiceDouble * range,
SpiceDouble * ra,
SpiceDouble * dec )
e l'implentazione possibile in Python potrebbe essere:
... _recrad_c = ct.shared_lib.recrad_c ... def recrad_c(rec_vector): _recrad_c.argtypes = [ct.c_double *3, ct.POINTER(ct.c_double), ct.POINTER(ct.c_double), ct.POINTER(ct.c_double)] _vec = ct.c_double * 3 vec = _vec(rec_vector[0], rec_vector[1], rec_vector[2]) r = ct.c_double(0) right_ascension = ct.c_double(0) declination = ct.c_double(0) _recrad_c(vec, r, right_ascension, declination) return (r.value, right_ascension.value, declination.value)
che non dovrebbe necessitare di spiegazioni, abbiamo già visto il procedimento
Proviamo a trasformare le coordinate rettangolari in sferiche. La funzione recrad_c restituisce valori espressi in radianti, quindi dobbiamo convertire l'ascensione retta in ore ( cioè in gradi / 15 ) e la declinazione in gradi. Ci torna utile la costante dpr_c() che è il fattore di conversione gradi/radianti e la funzione ddd2dms del mio modulo tool.py.
Ecco codice e risultato:
timestring = "2017-08-10T18:53:22" ref = "J2000" target = "SATURN_BARYCENTER" observer = "EARTH" abcorr = "NONE" xyz = jpl.spkezr_c(target, ref, timestring, abcorr, observer) elements = ["x :", "y :", "z :", "vx :", "vy :", "vz :"] for i in range(6): print elements[i], xyz[i] recvector = (xyz[0], xyz[1], xyz[2]) r, RA, dec = jpl.recrad_c(recvector) print string = "distanza {0} asc.retta {1} declinazione {2}" print string.format(r, ddd2dms(RA * jpl.dpr_c()/15.), ddd2dms(dec * jpl.dpr_c()) ) ... x : -218912917.201 y : -1296275066.93 z : -528970291.665 vx : -10.2998651768 vy : -20.7463383194 vz : -9.37634837706 distanza 1417060861.05 asc.retta (17, 21, 39) declinazione (-21, 55, 7)
Il sito Horizon mi dà questo risultato:
*******************************************************************************
Ephemeris / WWW_USER Sun Aug 6 01:03:28 2017 Pasadena, USA / Horizons
*******************************************************************************
Target body name: Saturn (699) {source: sat389}
Center body name: Earth (399) {source: sat389}
Center-site name: GEOCENTRIC
*******************************************************************************
Start time : A.D. 2017-Aug-10 18:53:22.0000 UT
Stop time : A.D. 2017-Aug-11 18:53:22.0000 UT
Step-size : 1440 minutes
*******************************************************************************
Target pole/equ : IAU_SATURN {East-longitude -}
Target radii : 60268.0 x 60268.0 x 54364.0 km {Equator, meridian, pole}
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : High-precision EOP model {East-longitude +}
Center radii : 6378.1 x 6378.1 x 6356.8 km {Equator, meridian, pole}
Target primary : Sun
Vis. interferer : MOON (R_eq= 1737.400) km {source: sat389}
Rel. light bend : Sun, EARTH {source: sat389}
Rel. lght bnd GM: 1.3271E+11, 3.9860E+05 km^3/s^2
Atmos refraction: NO (AIRLESS)
RA format : HMS
Time format : CAL
EOP file : eop.170803.p171022
EOP coverage : DATA-BASED 1962-JAN-20 TO 2017-AUG-03. PREDICTS-> 2017-OCT-21
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s
Table cut-offs 1: Elevation (-90.0deg=NO ),Airmass (>38.000=NO), Daylight (NO )
Table cut-offs 2: Solar elongation ( 0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
Table cut-offs 3: RA/DEC angular rate ( 0.0=NO )
*******************************************************************************
Date__(UT)__HR:MN:SS R.A._(ICRF/J2000.0)_DEC
*************************************************
$$SOE
2017-Aug-10 18:53:22 17 21 39.01 -21 55 05.9
2017-Aug-11 18:53:22 17 21 33.15 -21 55 11.4
$$EOE
*******************************************************************************
La corrispondenza dei risultati è molto buona, la minima differenza nella declinazione (circa 1 secondo d'arco) va riferita forse a limiti del calcolo o più probabilmente alla scelta di un differente kernel o di un diverso NAIF ID per Saturno. Tuttavia finora abbiamo fatto tutto per avere risultati sul piano equatoriale, che è il più usato dagli astronomi e nulla sul piano eclittico, che invece è prediletto dagli astrologi. Il passo successivo, sempre con gli strumenti del jpl toolkit, sarà cambiare frame di riferimento e passare dall'ascensione retta e declinazione a longitudine e latitudine eclittica. A fra un po'.
No comments:
Post a Comment