Primo collaudo di CSPICE

La directory test che abbiamo creato dovrebbe avere questo contenuto:

usefultest$ ls -la
totale 21740
drwxr-xr-x  2 ubuntu ubuntu     4096 ago  5 16:31 .
drwxr-xr-x 10 ubuntu ubuntu     4096 ago  5 16:30 ..
-rwxr-xr-x  1 ubuntu ubuntu  1140778 apr  8 03:47 brief
-rwxrwxrwx  1 ubuntu ubuntu 16788480 gen 10  2016 de421.bsp
-rw-r--r--  1 ubuntu ubuntu        0 ago  5 08:19 __init__.py
-rw-r--r--  1 ubuntu ubuntu     6780 ago  4 19:04 jpltoolkit.py
-rwxr-xr-x  1 ubuntu ubuntu  4286040 ago  3 08:41 libcspice.so
-rw-r--r--  1 ubuntu ubuntu     5257 lug 15  2016 naif0012.tls
-rw-r--r--  1 ubuntu ubuntu     1732 ago  5 12:36 testjpl.py
-rwxrwxrwx  1 ubuntu ubuntu    10146 giu 26 09:55 tools.py

In sostanza:

• un file con estensione bsp (file DAF/SPK, contenente i coefficienti del polinomio interpolatore per il calcolo delle effemeridi in un certo un arco temporale)
• l'utilità brief, che ci fornisce le informazioni relative al file di cui sopra
• il kernel generico naif0012.tls, che contiene i leapsecond stabiliti dall'IERS (va sostituito ogni volta che ne viene emesso uno nuovo)
• Il file di libreria shared libcspice.so
• il file che sto costruendo jpltoolkit.py (chiamatelo come volete) che contiene le funzioni ctypes
• il file tools.py (funzioni di calendario, trigonometriche e di supporto alla conversione radianti - gradi)
• il file __init__.py che serve esclusivamente a dichiarare che tools.py e jpltoolkit.py sono moduli
• il file testjpl.py che conterrà il programma di collaudo vero e proprio

cominciamo a scrivere il file di collaudo (testjpl.py, nel mio caso). Importiamo i moduli e lanciamo in esecuzione.

from tools import *
import jpltoolkit as jpl

print "valore di pi greco pi_c :", jpl.pi_c()
print "coefficiente di conversione radianti->gradi dpr_c :", jpl.dpr_c()
print "secondi in un giorno spd_c :", jpl.spd_c()

Python 2.7.12 (default, Nov 19 2016, 06:48:10) 
[GCC 5.4.0 20160609] on linux2
Type "copyright", "credits" or "license()" for more information.
>>> 
 RESTART: /home/ubuntu/Scrivania/ASTRONOMY/JPL toolkit/cspice/usefultest/prova.py 
valore di pi greco pi_c : 3.14159265359
coefficiente di conversione radianti->gradi dpr_c : 57.2957795131
secondi in un giorno spd_c : 86400.0
>>> 

Fin qui tutto bene. Le funzioni relative alle costanti restituiscono effettivamente i valori attesi. Ora carichiamo i nostri kernel usando la funzione furnsh_c di CSPICE.

import ctypes as ct

ct.shared_lib = ct.CDLL('./libcspice.so')

# definitions
_spd_c = ct.shared_lib.spd_c
_pi_c = ct.shared_lib.pi_c
_dpr_c = ct.shared_lib.dpr_c
_furnsh_c = ct.shared_lib.furnsh_c
_unload_c = ct.shared_lib.unload_c

def spd_c():
    _spd_c.restype = ct.c_double
    return _spd_c()

def pi_c():
    _pi_c.restype = ct.c_double
    return _pi_c()

def dpr_c():
    _dpr_c.restype = ct.c_double
    return _dpr_c()

def furnsh_c(k):
    _furnsh_c.argtypes=[ct.c_char_p]
    _furnsh_c(k)

def unload_c(k):
    _unload_c.argtype=[ct.c_char_p]
    _unload_c(k)

Come ho fatto per le costanti, creo una funzione ctypes per linkare la funzione CSPICE furnsh_c (che carica un kernel) e la funzione unload_c (che invece lo rimuove dalla memoria, quindi liberando spazio).

Molte istruzioni d'uso sono nella documentazione in particolare nella pagina delle funzioni più frequentemente utilizzate che è una vera miniera.

In particolare per l'uso di furnsh_c e unload_c ci guida alla signature delle stesse:

void furnsh_c ( ConstSpiceChar  * file ) 
void unload_c ( ConstSpiceChar  * file ) 

che significa che queste funzioni accettano un parametro di input (il nome del file nel formato ConstantSpiceChar, che corrisponde al char* di C, puntatore all'inizio della stringa)

In ctype basta banalmente dichiarare la variabile stringa come ctypes.c_char_p, o, nel mio caso, ct.c_char_p, visto che ho usato un'alias. Queste funzioni non restituiscono valore (sono void). La gestione di eventuali errori va gestita separatamente, ma ne parliamo più avanti.

Siamo pronti per inserire una funzione importante, str2et_c, che assume una stringa rappresentante una data e ora precise e le converte in secondi trascorsi dal J2000 (ore 12 del 1 gennaio 2000) secondo la convenzione TDB - Barycentric Dynamic Time, istituita nel 1976 a sostituzione del precedente Ephemeris Time. La conversione è automatica.

Integro la nuova funzione nel file jpltoolkit.py, evito di riprodurre le funzioni precedenti, spero che sia abbastanza chiaro il procedimento.

...
_str2et_c = ct.shared_lib.str2et_c
....
def str2et_c(s, e):
    _str2et_c.argtypes=[ct.c_char_p, ct.POINTER(ct.c_double)]
    _str2et_c(s, e)
    return e.value

Stavolta la signature è la seguente:

void str2et_c ( ConstSpiceChar * str, SpiceDouble * et ) 

Questo significa che la funzione accetta una stringa (che in C viene dichiarata come char *) e un puntatore a variabile double. Il che spiega perchè negli argomenti della funzione python equivalente _str2et uso un ct.c_char_p come prima e un ct.POINTER(ct.c_double). La funzione mi restituirà void, cioè None in linguaggio Python, ma io prevedo di restituire al chiamante il value della variabile, cioè il contenuto puntato dal puntatore a double. Sembra complicato, ma vi assicuro che ci si fa presto l'abitudine.

Una curiosità: tra UTC e TDB c'è una differenza temporale significativa. Se lanciamo la chiamata alla funzione str2et_c per il J2000 espresso in UTC, otteniamo questa risposta:

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"

(UTC 2000_01_01 12:00) 64.1839272847 seconds from J2000

che è esattamente quanto atteso, vedi questa tabella

Il formato della stringa è un formato standard secondo ISO. La pagina che documenta la funzione str2et_c ne propone diversi formati nel paragrafo Examples

Chiudo qui per il momento, continuero' nel prossimo post.