Wednesday, March 6, 2019

Caricare i kernel - CSPICE data types - Ulteriori funzioni

Proseguiamo l'estensione delle funzioni del toolkit accessibili tramite JNA. Come già visto nei post dedicati al linguaggio Python, il toolkit fornisce due metodi che possiamo utilizzare per il caricamento e per la dismissione dei kernel. I kernel sono i file dati utilizzati dal toolkit, ce ne sono di vario tipo, per una disamina completa potete leggere questa pagina della documentazione. Per ora ci intereressa caricare il kernel delle effemeridi (de421.bsp o simili) e il naif0012.tls per i leapsecond. Quindi aggiungeremo in seguito il kernel per specificare un frame di riferimento customizzato (trueepoch.tf)

Nel piccolo file interfaccia creato nel post precedente era già integrata la dichiarazione delle due funzioni di libreria furnsh_c e unload_c, ora ne scriveremo la definizione.

Partiamo dalle signature:


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

In Java le tradurremo come:


        void furnsh_c(String s);
        void unload_c(String s);

e la corrispondente definizione sarà:


    public void furnsh_c(String s){    
        INSTANCE.furnsh_c(s);        
        return;
    }

    public void unload_c(String s){
        INSTANCE.unload_c(s);
        return;
    }

Rivediamo il file CSpice.java e Main.java modificati:

import com.sun.jna.*;
import com.sun.jna.ptr.*;

public class CSpice {
    private static CLibrary INSTANCE;

    static {
        INSTANCE = Native.load("cspice", CLibrary.class);
    }

    public interface CLibrary extends Library {
        double pi_c();
        double spd_c();
        double dpr_c();
        void furnsh_c(String s);
        void unload_c(String s);

    }
    public double pi_c() {
        return INSTANCE.pi_c();
    }

    public double spd_c() {
        return INSTANCE.spd_c();
    }

    public double dpr_c() {
        return INSTANCE.dpr_c();
    }

    public void furnsh_c(String s){    
        INSTANCE.furnsh_c(s);        
        return;
    }

    public void unload_c(String s){
        INSTANCE.unload_c(s);
        return;
    }

}

e
public class Main {

     public static void main(String[] args) {
        CSpice spice = new CSpice();
        System.out.println("Valore di Pi greco: " + spice.pi_c());
        System.out.println("Secondi in un giorno: " + spice.spd_c());
        System.out.println("Numero di gradi per radiante: " + spice.dpr_c());
        spice.furnsh_c("de421.bsp");
        spice.furnsh_c("naif0012.tls");
        System.out.println("OK Kernel caricati");
        spice.unload_c("de421.bsp");
        spice.unload_c("naif0012.tls");
        System.out.println("OK Kernel dismessi");        
     }
}

Ricompiliamo il tutto ed eseguiamo:

ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ javac -classpath .:jna-5.2.0.jar -g Main.java
ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ java -classpath .:jna-5.2.0.jar Main
Valore di Pi greco: 3.141592653589793
Secondi in un giorno: 86400.0
Numero di gradi per radiante: 57.29577951308232
OK Kernel caricati
OK Kernel dismessi
ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ 

Siamo pra pronti ad inserire la chiamata alla funzione che, data una stringa di formattazione di data e ora in formato ISO 8601, restituisce il numero di secondi TDB (Temps Dynamique Barycentrique) a partire dalla epoch J2000.

La signature della funzione è la seguente:


void str2et_c ( ConstSpiceChar * str,
                SpiceDouble    * et   )

dove str è la data in formato ISO, et è il puntatore al valore double fornito come parametro e che va recuperato come valore.

per i puntatori abbiamo già considerato l'import delle classi del blocco com.sun.jna.ptr (seconda riga del file CSpice). La traduzione in Java sarà la seguente:


       double str2et_c(String s, DoubleByReference p);

per la dichiarazione e


    double str2et_c(String s, DoubleByReference p){
        INSTANCE.str2et_c(s,p);
        return p.getValue();
    }

per la definizione.

In buona sostanza: la funzione str2et_c, che nel toolkit restituisce un void, viene trasformata da noi in modo da restituire un double, che corrisponde al nuovo contenuto della variabile p, modificato dalla funzione ed estratto (dereferenziato, con la terminologia usata dal linguaggio C) con il metodo getValue(). Molto simile, nel procedimento, a quanto già fatto con python ctypes.

Rivediamo i file CSpice.java e Main.java con le ultime modifiche:

import com.sun.jna.*;
import com.sun.jna.ptr.*;

public class CSpice {
    private static CLibrary INSTANCE;

    static {
        INSTANCE = Native.load("cspice", CLibrary.class);
    }

    public interface CLibrary extends Library {
        double pi_c();
        double spd_c();
        double dpr_c();
        void furnsh_c(String s);
        void unload_c(String s);
        double str2et_c(String s, DoubleByReference p);

    }
    public double pi_c() {
        return INSTANCE.pi_c();
    }

    public double spd_c() {
        return INSTANCE.spd_c();
    }

    public double dpr_c() {
        return INSTANCE.dpr_c();
    }

    public void furnsh_c(String s){    
        INSTANCE.furnsh_c(s);        
        return;
    }

    public void unload_c(String s){
        INSTANCE.unload_c(s);
        return;
    }

    double str2et_c(String s, DoubleByReference p){
        INSTANCE.str2et_c(s,p);
        return p.getValue();
    }

}
import com.sun.jna.ptr.DoubleByReference;

public class Main {

     public static void main(String[] args) {
        CSpice spice = new CSpice();
        System.out.println("Valore di Pi greco: " + spice.pi_c());
        System.out.println("Secondi in un giorno: " + spice.spd_c());
        System.out.println("Numero di gradi per radiante: " + spice.dpr_c());
        spice.furnsh_c("de421.bsp");
        spice.furnsh_c("naif0012.tls");
        System.out.println("OK Kernel caricati");

        DoubleByReference et = new DoubleByReference();
        System.out.println("Secondi dal J2000: "+spice.str2et_c("2000-01-01T12:00:00", et));
        
        spice.unload_c("de421.bsp");
        spice.unload_c("naif0012.tls");
        System.out.println("OK Kernel dismessi");        
     }
}

E ricompiliamo

ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ javac -classpath .:jna-5.2.0.jar -g Main.java
ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ java -classpath .:jna-5.2.0.jar Main
Valore di Pi greco: 3.141592653589793
Secondi in un giorno: 86400.0
Numero di gradi per radiante: 57.29577951308232
OK Kernel caricati
Secondi dal J2000: 64.18392728473108
OK Kernel dismessi
ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ 

Tuesday, March 5, 2019

Marshalling di alcune funzioni riferite a costanti fondamentali

Come ho fatto per il linguaggio Python, scrivero' per prime alcune funzioni in Java finalizzate al recupero di valori di costanti matematiche o astronomiche di rilievo (pi greco, numero di secondi in un giorno, gradi per radiante.

Le funzioni del toolkit che consentono queste operazioni sono le seguenti, con la loro signature completa:

SpiceDouble pi_c  ( void )
SpiceDouble spd_c ( void )
SpiceDouble dpr_c ( void )

Nel file CSpice.java introdurro' le seguenti modifiche:

import com.sun.jna.*;
import com.sun.jna.ptr.*;

public class CSpice {
    private static CLibrary INSTANCE;

    static {
        INSTANCE = Native.load("cspice", CLibrary.class);
    }

    public interface CLibrary extends Library {
        double pi_c();
        double spd_c();
        double dpr_c();
        void furnsh_c(String s);
        void unload_c(String s);

    }
    public double pi_c() {
        return INSTANCE.pi_c();
    }

    public double spd_c() {
        return INSTANCE.spd_c();
    }

    public double dpr_c() {
        return INSTANCE.dpr_c();
    }

}

Nell'interface ho definito le funzioni traducendo i tipi di variabili SpaceDouble del toolkit nel tipo double di Java e il tipo void del parametro delle funzioni in uno spazio vuoto. In questo caso il marshalling è particolarmente semplice. La stesura effettiva delle funzioni relative è appena poco più estesa: le funzioni prima dichiarate diventano metodi della INSTANCE della libreria nativa cspice.

Per finire modifichiamo opportunamente la classe Main.java per estrarre i valori delle costanti dichiarate tramite i metodi associati:

public class Main {

     public static void main(String[] args) {
        CSpice spice = new CSpice();
        System.out.println("Valore di Pi greco: " + spice.pi_c());
        System.out.println("Secondi in un giorno: " + spice.spd_c());
        System.out.println("Numero di gradi per radiante: " + spice.dpr_c());

     }
}

A questo punto ricompiliamo il tutto:

ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ javac -classpath .:jna-5.2.0.jar -g Main.java
ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ java -classpath .:jna-5.2.0.jar Main
Valore di Pi greco: 3.141592653589793
Secondi in un giorno: 86400.0
Numero di gradi per radiante: 57.29577951308232
ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ 

Come è facile vedere, il risultato è perfettamente corrispondente a quello atteso. Il più è fatto. Da qui in poi estenderemo le funzioni usate fino a riprodurre quanto abbiamo già realizzato con il linguaggio python. Nel post successivo iniziamo a definire funzioni un po' più complesse

JNA e jpl toolkit - creazione dell'interfaccia

Per chiamare da Java le funzioni della shared library, dobbiamo imparare le fondamentali tecniche di interfacciamento e di marshalling (cioè di definizione di formati dati di interscambio tra codice gestito (Java) e non gestito (libreria C).

Prima di tutto scarichiamo la libreria in formato .jar dell'ultima versione disponibile di JNA, che troviamo QUI.


Cominciando dall'interfaccia, utilizzeremo un template, cioè un po' di codice boilerplate, come lo chiamano gli Americani, per renderci la vita più semplice.

import com.sun.jna.*;
import com.sun.jna.ptr.*;

public class CSpice {
    private static CLibrary INSTANCE;

    static {
        INSTANCE = Native.load("cspice", CLibrary.class);
    }

    public interface CLibrary extends Library {
//<-----------------definizione dei prototipi delle funzioni------------>
    }
//<-----------------chiamata alle funzioni, una alla volta-------------->
}

Spiegazione: l'import riguarda il core della libreria JNA e il complesso delle classi per la gestione dei puntatori, lo vediamo nel corso dello sviluppo. La classe pubblica CSpice conterrà tutto il codice per la creazione di un oggetto che sarà l'interfaccia da usare. Per il momento dò per scontato il resto del codice, fatta eccezione per la definizione statica dell'INSTANCE come
Native.Load("cspice", CLibrary.class).
Qui la libreria libcspice.so è richiamata esplicitamente, togliendo dal nome il prefisso lib e l'estensione .so. Questo vale per qualsiasi altra situazione in cui una generica libreria .so è utilizzata: libHelloWorld.so -> INSTANCE = Native.load("HelloWorld", CLibrary.class);.

Buttiamo giù un veloce file di classe Main.java, giusto per ospitare un metodo main:

public class Main {

     public static void main(String[] args) {
        CSpice spice = new CSpice();
     }
}

Facciamo una compilazione di prova: ricordiamoci di richiamare esplicitamente la libreria jna con l'opzione classpath:

javac -classpath .:jna-5.2.0.jar -g Main.java

Come potete apprezzare, la classe CSpice che contiene l'interfaccia viene richiamata implicitamente durante la compilazione di Main.java e il file di classe corrispondente viene creato come CSpice.class nella directory di progetto. Se facciamo il listing della directory vediamo i file seguenti:

ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ ls
 backup.java           CSpice.java     libcspice.so   naif0012.tls
'CSpice$CLibrary.class'   de421.bsp       Main.class     trueepoch.tf
 CSpice.class          jna-5.2.0.jar   Main.java
ubuntu@ubuntu-desktop:~/Scrivania/CSpice$

Proviamo a lanciare il programma con il comando seguente:

java -classpath .:jna-5.2.0.jar Main

Se il comando non produce output vuol dire che va tutto bene

Mi fermo qui, continuiamo nel prossimo post.

JAVA, JNA, codice unmanaged e marshalling

Lavorare su codice nativo non gestito da un ambiente operativo diverso da quello di produzione non è sempre facile, ma spesso non avete alternative. Pensate al nostro JPL toolkit, scritto in C e Fortran. Avreste voglia di riscrivere tutto in Python o in Java? sono centinaia di funzioni, non è il caso. Questo vale per molte situazioni diverse, pensate alle librerie scientifiche compilate in Fortran negli anni dei pionieri, o a tutto il codice generato dalle vecchie versioni Visual Studio (non so se interoperabilità COM vi ricorda qualcosa). Anche la nostra piccola shared library compilata qualche post fa, libcspice.so o la sua equivalente in Windows libcspice.dll rischia di restare un bel giocattolo se non troviamo il modo di riutilizzarla con i compilatori o interpreti scritti in altri linguaggi. Abbiamo esplorato un po' ctypes di Python, ora voglio iniziare, come accennato nel post precedente, a JNA di Java.

Per questo piccolo progetto non usero', almeno prima di un certo numero di post, un IDE grafico come Eclipse o Netbeans, preferisco, per semplicità e beneficio didattico, fare tutte le compilazioni da riga di comando. Userò, come compilatore C per JPL SPICE il GNU gcc (ma trovate tutto già spiegato nei post precedenti) e come compilatore java openjdk (una qualsiasi versione successiva alla 7 va bene, attualmente nella mia Ubuntu 18.10 uso la versione 11.0.2). Se operate in ambiente windows, vista l'elevata portabilità di Java, potremo facilmente usare la linea di comando in una finestra terminale, se volete un compilatore C per questa piattaforma avete ampia libertà di scelta, c'è Visual Studio, oppure Mingw che riproduce un ambiente POSIX sotto Windows, vedete voi cosa preferite, alla fine di tutto faro' io un post apposito).

La prima cosa che faremo sarà creare una directory dedicata al progettino Java sul desktop o in qualsiasi posizione comoda per lavorare. All'interno collocheremo alcuni file creati precedentemente, come vi mostra lo schemino seguente.

ubuntu@ubuntu-desktop:~/Scrivania/CSpice$ ls -l
totale 20596
-rw-rw-r-- 1 ubuntu ubuntu 16788480 mar  2 10:56 de421.bsp
-rwxrwxr-x 1 ubuntu ubuntu  4281824 mar  2 10:45 libcspice.so
-rw-rw-r-- 1 ubuntu ubuntu     5257 mar  2 11:00 naif0012.tls
-rwxrwxrwx 1 ubuntu ubuntu      598 ago 12  2017 trueepoch.tf
ubuntu@ubuntu-desktop:~/Scrivania/CSpice$

Riepilogando:

libcspice.so è la libreria condivisa che contiene le funzioni del JPL toolkit (per un accesso veloce ai metodi di compilazione partite da Questo post in avanti)

de421.bsp è il file delle effemeridi fornito dal sito ftp della NASA. Se volete un file più recente c'è ampia scelta, il 421 ha il pregio di essere abbastanza leggero e portabile anche in ambiente Android.

naif0012.tls è il kernel generico che contiene l'elenco aggiornato dei leapsecond (secondo intervalare in italiano), ricordate che si tratta della ".. correzione del tempo applicata al tempo coordinato universale (UTC) per mantenerlo allineato al giorno solare medio, l'ente responsabile della regolazione è l'International Earth Rotation and Reference Systems Service (IERS) (sorgente Wikipedia). Trovate l'elenco dei kernel disponibili in questa pagina.

Infine trueepoch.tf è il kernel creato da me seguendo le istruzioni fornite dal NAIF in questa pagina , che servirà più avanti per svincolarsi dal frame di riferimento inerziale e spostarsi su un frame dinamico come quello che usano gli astrologi. Trovate i dettagli per la sua costruzione nei miei post precedenti già citati.

A questo punto possiamo iniziare a scrivere un po' di codice Java. Per semplicità metto tutto nella stessa directory, rinunciando a creare un package, nulla vieta di strutturare meglio il progetto, lo faremo in seguito.
Dopo questa lunga premessa, passo al post successivo per la scrittura delle prime classi Java che ci servono.

Monday, March 4, 2019

JPL Spice Toolkit, Java, JNA - come utilizzare la libreria nativa libcspice.so in ambiente Java

Sono due anni che prometto di occuparmi di domificazione e ogni volta trovo una scusa per rinviare. Sarà che non sento molto il problema, sarà che ce n'è fin troppi di sistemi di domificazione, sarà che il più usato (Placidus) ha un comportamento degenere alle latitudini estreme, per un motivo o per l'altro trovo sempre qualcos'altro di cui occuparmi.

Da qualche giorno sto pasticciando con JAVA (vi ricordo che non sono un professionista né dell'informatica né dell'astrologia), l'ho usato un po' in passato per ambienti desktop e un po' per creare qualcosa per Android. Mi chiedevo se era possibile usare una libreria nativa come la libreria dinamica libcspice.so che avevo compilato qualche post addietro, in piattaforma JAVA, un po' come abbiamo fatto con la libreria ctypes di Python.

Il metodo effettivamente esiste, anzi ne esiste più di uno, c'è SWIG, c'è JNI, ma poi ho cominciato a usare JNA e ha funzionato abbastanza presto. Ho, pertanto, deciso di scrivere un piccolo tutorial su JNA, visto che quello che si trova in giro, anche in ambienti seri come Stack Overflow è poco più di Hello World. (per inciso, se un giorno scoppierà una rivoluzione, sarà per l'esasperazione dei programmatori che si trovano sempre un HelloWorld inutile ogni volta che cercano qualche informazione sul web).

Se andate indietro di qualche post e rileggete quello che ho scritto su ctypes, vi renderete conto che è una libreria non di immediata comprensione, ma che, una volta imparate le basi, diventa estremamente soddisfacente. La stessa cosa vale per JNA, interfaccia in puro linguaggio JAVA per il binding di funzioni nascoste in una libreria di sistema, sia essa una DLL di Windows o una SO di Linux. Se avrete la pazienza di seguirmi, ho qualche piccolo trucco da esporre in chiaro.

Per rendere le cose più semplici, partirò, come ho fatto per Python ctypes, con il call di costanti presenti nel JPL toolkit e aggiungerò progressivamente delle chiamate a funzioni, facendovi partecipi di questo accostamento iniziale a JNA. Mi perdonerete se l'approccio potrà sembrare poco professionale, ma, ripeto, sono un programmatore amatoriale.

Al prossimo post, penso di scriverlo a giorni, le mie routine sono sufficientemente mature per una breve esposizione. A presto.

Wednesday, May 16, 2018

Piccola deviazione dal tema principale: preparare la carne vegetale

A furia di tentativi ed errori sono riuscito a produrre, con mezzi semplici, un sostituto salutare ed economico della carne, che non sa di molto ma è abbastanza facile da realizzare, poco costosa, molto digeribile. Ricetta adatta all' ingresso di Urano nel Toro, che esprime, in ogni caso, un cambiamento nelle abitudini alimentari.

Senza ulteriori indugi, vi passo la mia ricetta personale.

Ricetta della carne vegetale

Ingredienti

  1. lenticchie secche (varietà economica da discount) g.400
  2. farina manitoba (o comunque ad alto contenuto proteico) Kg 2
  3. Semi di varia estrazione (girasole, lino, sesamo, papavero) g. 200
  4. acqua (300ml per 500g di farina): 1200ml
  5. sale, aromi a piacere, due cucchiai di olio d'oliva, peperoncino macinato, pepe, un cucchiaio di curcuma

Attrezzatura

  1. una pentola grande
  2. un recipiente in plastica di grandi dimensioni (per due kg di pasta)
  3. uno scolapasta in plastica o acciaio di misure corrispondenti al recipiente di cui sopra
  4. un telino bianco pulito e deodorato (se è fresco di bucato, farlo bollire dieci minuti in acqua)
  5. un miscelatore o robot da cucina o impastatrice
  6. recipienti più piccoli, piatti e stoviglie al bisogno

Preparazione

Formare il seitan (glutine)

Impastare tutta la farina con l'acqua, non aggiungere nè lievito nè sale

Formare una palla di pasta e metterla in un recipiente largo, totalmente immersa in acqua, per almeno 4 ore (serve a formare il glutine)

Al termine scolarla, inserirla nello scolapasta e inserire quest'ultimo nel recipiente di cui sopra, dopo averlo completamente scolato.

Con un filo d'acqua, alternando acqua calda e fredda, lavorare l'impasto per eliminare tutto l'amido.

Dopo quindici-venti minuti di lavorazione la pasta, che è diventata molto più morbida, diventa più densa e mantiene la forma e l'acqua di lavaggio da bianca opaca diventa madreperlacea e semitrasparente.

Quando non si riesce più a estrarre amido, mettere il seitan così ottenuto a scolare dentro lo scolapasta.

Preparare le farine di lenticchie e semi oleosi

Le farine possono essere preparate anche in un altro momento e conservate. Per macinare le lenticchie io uso un comune macinacaffè elettrico, ma va bene anche un frullatore. La farina che si ottiene è abbastanza grossa, tipo semola. I semi oleosi (girasole, lino, sesamo) vanno macinati con un po' di lenticchie (che sono più secche).

Preparazione dell'impasto

Far girare le farine con gli aromi, unire l'olio di oliva, fino ad ottenere un impasto granuloso e omogeneo. Continuare a miscelare e unire, poco alla volta, il seitan tagliato a pezzi grossi come un mandarino.

Non serve omogeneizzare troppo, il seitan fa un po' le veci del grasso animale, per cui è meglio se rimane un po' visibile.

Dopo la miscelazione (in base alle dimensioni della pentola si puo' fare una porzione o anche due) dare una forma grossolana all'impasto e chiuderlo dentro il telino bianco.

Legare con spago bianco alimentare lasciando un po' lasco perchè con la cottura cresce di volume.

Immergere il tutto in acqua salata, non servono altri aromi. Cuocere per un'ora. Al termine scolare, mettere su un tagliere di legno ad asciugare, dopo una o due ore tagliare a fette e conservare in congelatore quello che non si consuma subito.

Per non fare attaccare le fette tra loro durante il congelamento usare strisce di carta forno tra una fetta e l'altra.

Dopo scongelamento basta riscaldare, non serve ulteriore cottura. Non ha un gran sapore, tende a prendere quello del condimento, quindi si puo' usare del sugo di pomodoro o quello che la fantasia ispira.

Una o due fette vanno bene per un normale secondo, è molto digeribile.

Valore nutrizionale

Qualcosa si perde con la cottura, va considerata fisiologica una diminuzione della disponibilità finale di vitamine e sali minerali.

L'apporto di carboidrati è minimo, ma non è la fonte principale di approvvigionamento per questo componente

Dal punto di vista proteico, possiamo calcolare:

grammi/100g peso in grammi tipo alimento totale grammi
23 400 lenticchie secche 92
13 2000 farina manitoba 260
21 100 semi di girasole 21
18 100 altri semi 18

totale = 391g (19.55g / 100g)


le proteine della carne vegetale contengono tutti gli aminoacidi essenziali (per la combinazione di legumi e cereali). Una normale fetta di manzo del peso di circa 100g contiene mediamente 18g di proteine, per cui, almeno in peso, c'è una discreta equivalenza. Il costo di preparazione della carne vegetale, al lordo di tutti i costi, non supera i 2 euro per Kg.

il peso totale della carne vegetale a fine cottura è di circa 2 kg, quindi il contenuto finale proteico è 195.5g/1000g = circa il 19.55% del peso al consumo.

Come dicevo sopra, la miscela di amminoacidi essenziali è buona, puo' sostituire integralmente la carne animale. Bisogna ricordare pero' che la dieta rimane scoperta di vitamina B12, che va quindi integrata a parte, in caso di alimentazione vegana pura.

Non mi resta che augurarvi buona sperimentazione e buon appetito.

Wednesday, August 9, 2017

JPL Toolkit in Python ctypes per utenti Windows

Molto brevemente, non utilizzando abitualmente Windows sono poco esperto nell'uso dei compilatori per C e C++ in quest'ambiente operativo. Tuttavia, siccome la maggior parte degli utenti desktop usa Windows, mi dispiace non poter fornire anche a loro qualche utile suggerimento.

Per avere una dll compilata direttamente dai sorgenti del toolkit, consiglio l'installazione preventiva di:

  1. Python 32 bit, direttamente dal sito python.org
  2. il compilatore VC per Python,
  3. che si puo' ottenere seguendo le istruzioni di questo sito. Io ho installato il Python 2.7.13 per cui la versione più adatta è la 9.0. Dopo l'installazione rendete accessibile il path al compilatore modificando le variabili di ambiente per puntare il compilatore in linea di comando cl.exe. Questo si trova nella directory "C:\Users\\AppData\Local\Programs\Common\Microsoft\Visual C++ for Python\9.0\VC\bin", che bisognerà aggiungere alla variabile di ambiente PATH.

  4. Suggerisco caldamente di installare Spicepy di Andrew Annex,
  5. già citato qualche post fa, perchè consente di compilare automaticamente la dll desiderata durante l'installazione. Per come fare leggete con attenzione le istruzioni sul sito GitHub. Io ho preferito fare il download del file zip, decomprimerlo sul mio desktop, lanciare python setup.py install e attendere la compilazione. Trovate la dll desiderata nella sottodirectory SpicyPy-master\spicypy\utils, con il nome cspice.dll, copiatela nella directory di lavoro in cui, se avete seguito le mie precedenti istruzioni, dovreste avere tutti i file necessari. Ad ogni buon conto provo a riepilogarli:

    • de421.bsp o analoghi
    • naif0012.tls
    • trueepoch.tf
    • __init__.py (file vuoto)
    • jpltoolkit.py o come avete preferito chiamarlo
    • tools.py (piccola raccoltà di utilità di calendario e trigonometriche
    • testjpl.py (anche qui, il nome che preferite) che contiene i test della libreria

Naturalmente le chiamate saranno fatte alla .dll anzichè alla .so, ma la sintassi, per il resto, è perfettamente identica a quella degli script che ho scritto in linux mint

Una piccola raccomandazione, nella fretta ho dimenticato che al termine del lavoro bisogna fare l'unload dei kernel, altrimenti rimangono sospesi in memoria. Per facilitare le cose, vi ripubblico il file di test completo, cosi' come l'ho scritto in precedenza

Per prima cosa, nel file che gestisce la libreria (jpltoolkit.py) scrivete:

import ctypes as ct
ct.shared_lib = ct.CDLL('.\cspice.dll')

al posto di :

import ctypes as ct
ct.shared_lib = ct.CDLL('./libcspice.so')

A seguire il file di test completo:

from tools import *
import jpltoolkit as jpl

lista = ["SUN", "MOON", "MERCURY", "VENUS", "MARS",
         "JUPITER_BARYCENTER", "SATURN_BARYCENTER", "URANUS_BARYCENTER",
         "NEPTUNE_BARYCENTER", "PLUTO_BARYCENTER"]

# kernels load
jpl.furnsh_c('./de421.bsp')
jpl.furnsh_c('./naif0012.tls')
jpl.furnsh_c('./trueepoch.tf')

# calculate geocentric equatorial and equatorial coordinates
timestring = "2017-08-08T10:13:00"
ref = "MYTRUEEPOCH"
abcorr = "NONE"

print "equatorial coordinates (earthtrueeopch)"
for i in lista:
    target = i
    observer = "EARTH"
    # first equatorial
    xyz = jpl.spkezr_c(target, ref, timestring, abcorr, observer)
    x = xyz[0]; y = xyz[1]; z = xyz[2]
    rec_vector = (x,y,z)
    result = jpl.recrad_c(rec_vector)
    print "recrad", i, result[0], ddd2dms(result[1] * jpl.dpr_c()), ddd2dms(result[2] * jpl.dpr_c())
print
jpl.unload_c('./trueepoch.tf')
jpl.unload_c('./naif0012.tls')
jpl.unload_c('./de421.bsp')

Notate che ho corretto alcune cose, in particolare l'unload finale dei kernel

Come ho già affermato, non sono un esperto di windows, ma se doveste avere problemi, cercherò di esservi d'aiuto. Alla prossima

Aggiornamento

Ho fatto diverse prove con la dll compilata da spiceypy in vari ambienti operativi. Riesco a farla funzionare, al momento, solo con Python 32 bit. In alcuni casi non sono riuscito a compilare spiceypy in modo efficace con la procedura python setup.py install ma ho avuto successo con pip install spiceypy. Se avete problemi provate questo secondo metodo. La cspice.dll si trova dentro la directory C:\Python27\Lib\site-packages\spiceypy-2.0.1.dev0-py2.7.egg\spiceypy\utils

.

Tuesday, August 8, 2017

Coordinate eclittiche con JPL Toolkit - reference frame e file FK

NB: Questa pagina è stata interamente riscritta in data 25-03-2019, alcune delle considerazioni fatte all'epoca della scrittura erano inesatte e fuorvianti, me ne scuso con i lettori

L'astrologia non usa un sistema di riferimento inerziale. Inerziale vuol dire che gli oggetti, se non sono sottoposti a forze esterne, mantengono la loro velocità iniziale. Se compaiono accelerazioni, quindi forze impresse agli oggetti, il sistema di riferimento NON è inerziale. Il più inerziale fra quelli disponibili è il International Celestial Reference Frame (ICRF), adottato dall'International Astroniomical Union dal 1 gennaio 1998, in sostituzione del precedente FK5. E' quasi inerziale perchè basato sulla misura di più di 200 sorgenti extragalattiche nello spazio profondo.

Precisiamo inoltre che gli astrologi usano un sistema di riferimento che è assolutamente non inerziale, perchè non è basato su punti fissi dello spazio profondo, ma su un punto di origine nel cerchio dell'eclittica, il punto vernale o 0° di Ariete, dato dall'intersezione di equatore ed eclittica, nell'istante in cui la declinazione solare (in termini di coordinate equatoriali) torna ad essere positiva. Questo vale per tutta la superficie terrestre, anche se i cicli stagionali variano da emisfero nord a emisfero sud, creando espressioni climatiche differenti.

Il punto vernale è dato dall' intersezione dell'equatore (che esprime il moto di rotazione terrestre) e dell'eclittica (cioè dal moto di rivoluzione della Terra intorno al Sole), ma risente di almeno altre due componenti del moto. Le quattro principali sono la rotazione, la rivoluzione, la precessione (che è un moto periodico di circa 26000 anni in cui il Polo Nord celeste descrive un cerchio nel riferimento inerziale) e la nutazione, che è un'ondulazione di quest'ultimo cerchio, come appare come nella figura seguente:

Nel post precedente abbiamo utilizzato il frame J2000, che la documentazione di SPICE descrive come "Earth mean equator, dynamical equinox of J2000. The root reference frame for SPICE". Ai fini astrologici dobbiamo passare ad un sistema di riferimento eclittico, dinamico, cioè tempo dipendente, che integri tutte le componenti principali del moto. Il JPL toolkit ci dà la possibilità di costruire un kernel su misura, che costruiremo secondo le regole indicate nella documentazione. In primo luogo renderemo dinamico il frame, riferito pero' al sistema di riferimento eclittico.

Usiamo un template che ci fornisce la documentazione:

   FRAME_             =  
   FRAME__NAME          =  
   FRAME__CLASS         =  5
   FRAME__CLASS_ID      =  
   FRAME__CENTER        =  399
   FRAME__RELATIVE      = 'ECLIPJ2000'
   FRAME__DEF_STYLE     = 'PARAMETERIZED'
   FRAME__FAMILY        = 'MEAN_ECLIPTIC_AND_EQUINOX_OF_DATE'
   FRAME__PREC_MODEL    = 'EARTH_IAU_1976'
   FRAME__OBLIQ_MODEL   = 'EARTH_IAU_1980'
   FRAME__ROTATION_STATE= 'ROTATING'

In primo luogo, essendo un frame kernel, sarà un file con estensione .tf:

NAIF recommends kernel names use only lower case letters. 
NAIF further recommends one follows the conventions established for kernel 
name extensions, shown below.

            .bc    binary CK
            .bes   binary Sequence Component EK
            .bpc   binary PCK
            .bsp   binary SPK
            .tf    text FK
            .ti    text IK
            .tls   text LSK
            .tm    text meta-kernel (FURNSH kernel)
            .tpc   text PCK
            .tsc   text SCLK
 

quindi lo salveremo, per esempio, come "trueepoch.tf".

Quindi assegneremo un frame_ID, che dovrà essere tra quelli consentiti da NAIF: "If the frame is a TK frame, the class ID must match the frame ID. For both ID codes you should use a positive integer in the range from 1400000 to 2000000 (unless you are working in an official project capacity in which case you should ask NAIF to provide a CLASS_ID for you)". Potremo scegliere un numero come 1987654, e sostituirlo nel template dove previsto.

   FRAME_             =  1987654
   FRAME_1987654_NAME          =  
   FRAME_1987654_CLASS         =  5
   FRAME_1987654_CLASS_ID      =  1987654
   FRAME_1987654_CENTER        =  399
   FRAME_1987654_RELATIVE      = 'ECLIPJ2000'
   FRAME_1987654_DEF_STYLE     = 'PARAMETERIZED'
   FRAME_1987654_FAMILY        = 'MEAN_ECLIPTIC_AND_EQUINOX_OF_DATE'
   FRAME_1987654_PREC_MODEL    = 'EARTH_IAU_1976'
   FRAME_1987654_OBLIQ_MODEL   = 'EARTH_IAU_1980'
   FRAME_1987654_ROTATION_STATE= 'ROTATING'

Per il frame name sceglieremo qualcosa di indicativo, da richiamare nella funzione spkezr_c già tradotta in Python in un post precedente, per esempio "MYTRUEEPOCH"

   FRAME_MYTRUEEPOCH             =  1987654
   FRAME_1987654_NAME          =  'MYTRUEEPOCH'
   FRAME_1987654_CLASS         =  5
   FRAME_1987654_CLASS_ID      =  1987654
   FRAME_1987654_CENTER        =  399
   FRAME_1987654_RELATIVE      = 'ECLIPJ2000'
   FRAME_1987654_DEF_STYLE     = 'PARAMETERIZED'
   FRAME_1987654_FAMILY        = 'MEAN_ECLIPTIC_AND_EQUINOX_OF_DATE'
   FRAME_1987654_PREC_MODEL    = 'EARTH_IAU_1976'
   FRAME_1987654_OBLIQ_MODEL   = 'EARTH_IAU_1980'
   FRAME_1987654_ROTATION_STATE= 'ROTATING'

I primi 6-8 byte del file sono dedicati all'identificazione del tipo di kernel, per un kernel FK la dicitura prevista è KPL/FK, in una riga a se stante, seguito da una breve descrizione. Inoltre il blocco di definizioni deve essere compreso fra le istruzioni escaped \begindata e \begintext. Il risultato finale potrebbe essere:

KPL/FK

private Kernel FK,  dynamic, geocentric, time-based. It integrates precession, nutation and ecliptic obliquity models


\begindata

   FRAME_MYTRUEEPOCH           =  1987654
   FRAME_1987654_NAME          =  'MYTRUEEPOCH'
   FRAME_1987654_CLASS         =  5
   FRAME_1987654_CLASS_ID      =  1987654
   FRAME_1987654_CENTER        =  399
   FRAME_1987654_RELATIVE      = 'ECLIPJ2000'
   FRAME_1987654_DEF_STYLE     = 'PARAMETERIZED'
   FRAME_1987654_FAMILY        = 'MEAN_ECLIPTIC_AND_EQUINOX_OF_DATE'
   FRAME_1987654_PREC_MODEL    = 'EARTH_IAU_1976'
   FRAME_1987654_OBLIQ_MODEL   = 'EARTH_IAU_1980'
   FRAME_1987654_ROTATION_STATE= 'ROTATING'

\begintext

La classe 5 è prevista per i modelli dinamici, 'PARAMETERIZED' è una definizione di default, 'ROTATING' significa non inerziale. Non prevedo ulteriori informazioni testuali, per cui lascio vuoto lo spazio oltre \begintext.

Salviamo tutto come trueepoch.tf e cominciamo a usare il kernel.

Scriviamo un breve programma python che calcola le posizioni eclittiche di Sole, Luna e pianeti usando il nuovo modello per la data di oggi e per il mezzogiorno della data in cui sto scrivendo: 8-Ago-2017.

Codice e risultati riferiti a tale data:

from tools import *
import jpltoolkit as jpl

lista = [b"SUN", b"MOON", b"MERCURY", b"VENUS", b"MARS",
         b"JUPITER_BARYCENTER", b"SATURN_BARYCENTER", b"URANUS_BARYCENTER",
         b"NEPTUNE_BARYCENTER", b"PLUTO_BARYCENTER"]
# kernels load
jpl.furnsh_c(b'./de421.bsp')
jpl.furnsh_c(b'./naif0012.tls')

# calculate geocentric ecliptic coordinates
jpl.furnsh_c(b'./trueepoch.tf')
timestring = b"2013-08-08T12:00:00"
ref = b"MYTRUEEPOCH"
abcorr = b"NONE"

print ("ecliptic coordinates (earthtrueeopch) 2013-08-08 12:00:00 UTC")
for i in lista:
    target = i
    observer = b"EARTH"
    # first equatorial
    xyz = jpl.spkezr_c(target, ref, timestring, abcorr, observer)
    x = xyz[0]; y = xyz[1]; z = xyz[2]
    rec_vector = (x,y,z)
    result = jpl.reclat_c(rec_vector)
    a = result[1]
    while a < 0:
        a += math.pi*2
    print ("reclat", i, result[0], ddd2dms(a * jpl.dpr_c()), ddd2dms(result[2] * jpl.dpr_c()))
print("\n")
print("\n")
print("\n")
---------------------------
ecliptic coordinates (earthtrueeopch) 2013-08-08 12:00:00 UTC
('reclat', 'SUN', 151683411.43760258, (136, 6, 25), (0, 0, 1))
('reclat', 'MOON', 396667.3758218271, (154, 9, 47), (-4, 37, 28))
('reclat', 'MERCURY', 165672145.9315892, (119, 59, 32), (0, 14, 9))
('reclat', 'VENUS', 191961492.69616985, (170, 19, 16), (1, 4, 57))
('reclat', 'MARS', 355674573.1827142, (107, 18, 15), (0, 46, 30))
('reclat', 'JUPITER_BARYCENTER', 885254140.6884873, (99, 32, 37), (0, 9, 16))
('reclat', 'SATURN_BARYCENTER', 1493132563.8426943, (215, 36, 25), (2, 20, 15))
('reclat', 'URANUS_BARYCENTER', 2911577405.9178605, (12, 19, 43), (0, 42, 57))
('reclat', 'NEPTUNE_BARYCENTER', 4341288355.286687, (334, 27, 18), (0, 40, 55))
('reclat', 'PLUTO_BARYCENTER', 4737631830.616889, (279, 25, 10), (3, 8, 39))

Per confronto allego una copia della pagina corrispondente delle Rosycrucian Ephemeris 2000-2100 12h TDT (noon)- International Edition The Rosicrucian Felloship - Publisher - Oceanside). Si tenga comunque conto delle differenze di impostazione (TDT anzichè UTC, JPL DE102 anzichè JPL DE421, rotazione alla J2000 standard epoch, riduzione alle coordinate apparenti, deflessione della luce nel campo gravitazionale del Sole ecc.):

Sunday, August 6, 2017

Calcolo di posizione con CSPICE - JPL Toolkit

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'.

Saturday, August 5, 2017

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.

<h1>Magic: Will, Reality and Science</h1>

Conversation between me and Gemini Me: Gemini, let's hypothesize that magic, that is, the art of causing changes in accordance w...