Saturday, February 28, 2015

Widget Calendario con Tkinter - Tkinter Calendar Widget part 2

Siamo arrivati alla parte 'intelligente' del codice, cioè la gestione della data conseguente all'uso del mouse sulle spinbox o sui button del giorno.

In un post precedente ho parlato delle Variable Classes, cioè di quelle particolari variabili che sono collegate ai widget e determinano, in relazione a specifici eventi, di chiamare delle funzioni di callback.

Allego di seguito un'altra sezione del codice:

        # some Variable Classes, when updated a callback function is activated
        self.year = IntVar(self.Calendar)
        self.year.set(self.dtime.year)
        self.year.trace('w', self.check_and_reconfigure)
        self.month = IntVar(self.Calendar)
        self.month.set(self.dtime.month)
        self.month.trace('w', self.check_and_reconfigure)
        self.day = IntVar(self.Calendar)
        self.day.set(self.dtime.day)
        self.day.trace('w', self.check_and_reconfigure)
        self.hour = IntVar(self.Calendar)
        self.hour.set(self.dtime.hour)
        self.hour.trace('w', self.check_and_reconfigure)
        self.minute = IntVar(self.Calendar)
        self.minute.set(self.dtime.minute)
        self.minute.trace('w', self.check_and_reconfigure)
        # bottom Frame and call to functions
        self.month_length = 0
        self.spin_year.config(textvariable = self.year)
        self.spin_month.config(textvariable = self.month)
        self.spin_hour.config(textvariable = self.hour, justify = RIGHT)
        self.spin_minute.config(textvariable = self.minute, justify = RIGHT)
        self.jd = self.cal2jul(self.year.get(),
                                                 self.month.get(),
                                                 self.day.get(),
                                                 self.hour.get(),
                                                 self.minute.get())
        self.time_lbl.config(text = "{0}:{1}     JD {2}".format(
                                    str(self.hour.get()).zfill(2),
                                    str(self.minute.get()).zfill(2),
                                    self.jd))                   
        self.check_and_reconfigure()

Le variable classes riproducono i tipi fondamentali di variabili elementari di Python: boolean, string, integer, float. I tipi corrispondenti sono BooleanVar, StringVar, IntVar e DoubleVar. Si creano come istanze, assegnandole ad una nuova variabile, e ricevono un valore dato dal metodo .set della variable class. Per leggere il contenuto si usa invece il metodo .get(). L'assegnazione al widget comporta l'iscrizione della nuova variabile tra i parametri del widget stesso, attraverso l'istruzione 'textvariable', nelle modalità che abbiamo già incontrato:

label = Label(self.master, textvariable = <class variable>)
oppure con il metodo config: label.config(textvariable = <class variable>
oppure ancora con la selezione diretta della key: label['textvariable'] = <class variable> 

Il comportamento dinamico di queste variabili è fornito dal metodo .trace: <class variable>.trace('w', <funzione di callback> quando la variabile viene modificata, per esempio agendo sulle frecce delle spinbox, oppure inserendo direttamente il valore in un widget di tipo entry. Esiste anche l'attivazione del callback qualora si legga il contenuto della variabile, in tal caso il trace si fa con il parametro 'r' anziche 'w'.

Nel codice sopra riportato la funzione di callback è self.check_and_reconfigure, che vediamo immediatamente:

    def check_and_reconfigure(self, *args):
        year = self.year.get()
        month = self.month.get()
        day = self.day.get()
        hour = self.hour.get()
        minute = self.minute.get()
        self.act_time=""
        # recalculate month's length in days
        if month in (1,3,5,7,8,10,12):
            self.max_days = 31
        elif month in (4,6,9,11):
            self.max_days = 30
        elif month == 2:
            if year < 1582:
                if year%4 == 0:
                    self.max_days = 29
                else:
                    max_day = 28
            else:
                if (year%4 == 0) and (year%400 == 0):
                    self.max_days = 29
                else:
                    self.max_days = 28
        self.month_length = self.max_days
        # fill days' buttons
        for i in self.table:
            i.config(text = '')
            i.config(bg   = 'grey')
            i.config(state = DISABLED)
        for day in range (1, self.month_length+1):
            dotw = self.day_of_the_week(self.year.get(), self.month.get(), day)
            posx = dotw
            if day == 1:
                lag = posx
            posy = int((day+lag-1)/7)
            self.table[posx + posy * 7].config(text = str(day))
            self.table[posx + posy * 7].config(state = NORMAL)
            if day == self.day.get():
                self.table[posx + posy * 7].config(bg='white')
                self.act_time = "{0}:{1}     JD {2}".format(
                        str(self.hour.get()).zfill(2), str(self.minute.get()).zfill(2),
                        self.cal2jul(self.year.get(), self.month.get(), self.day.get(),
                        self.hour.get(), self.minute.get()))
        self.time_lbl['text']=self.act_time

Nell'ordine vengono eseguite le seguenti operazioni:

  1. Si leggono, nell'ordine, i valori della variabile self.dtime riferita all'anno, al mese, all'ora e al minuto
  2. con un semplice algoritmo si determina la lunghezza, in giorni, del mese selezionato, sia per le date gregoriane che per quelle più antiche (qui ci sarà un problema, il calendario gregoriano non è stato adottato da tutti i paesi nel 1582, bisognerà tenerne conto)
  3. si ripulisce la table dei giorni, disabilitando, temporaneamente, la possibilità di cliccare sui button
  4. con un altro semplice algoritmo si calcola, per i giorni da 1 alla lunghezza del mese, la posizione del giorno nella table in riga e colonna e si ripristina la possibilità di interazione per il mouse per i soli button selezionati
  5. si crea la stringa che compare in fondo al widget con l'ora e minuto e la data giuliana

Rimane la parte finale del codice, che segue immediatamente:

    def change_day(self, button):
        self.day.set(button.cget('text'))
        
        
    def cal2jul(self, year, month, day,
                hour=0, minute =0):
        month2 = month
        year2 = year
        if month2 <= 2:
            year2 -= 1
            month2 += 12
        else:
            pass
        if (year*10000 + month*100 + day) > 15821015:
            a = int(year2/100)
            b = 2 - a + int(a/4)
        else:
            a = 0
            b = 0
        if year < 0:
            c = int((365.25 * year2)-0.75)
        else:
            c = int(365.25 * year2)
        d = int(30.6001 *(month2 + 1))
        return b + c + d + day + hour / 24.0 + minute / 1440.0 + 1720994.5

   
    def day_of_the_week(self, year, month, day):
                
        jd = self.cal2jul(year, month, day)
        a = (jd+1.5)/7
        f = int((a % 1)*7 +.5)
        return f        

Chiudiamo con alcune funzioni helper:

self.change_day(self) banalmente effettua l'aggiornamento della variabile self.day, abbiamo visto come la funzione lambda consenta di aggiornare il valore della variabile associata ad uno specifico button

cal2jul è una della funzioni che abbiamo usato nel calcolo di posizione dei pianeti. Ho preferito integrarla anzichè richiamarla da un modulo esterno per dare indipendenza al nuovo widget

infine day_of_the_week serve per dare un valore numerico ad ogni giorno del mese, è una funzione che utilizza la data giuliana.

Con questo ho finito, appuntamento ad un prossimo post

Widget Calendario con Tkinter - Tkinter Calendar Widget

Ho abbandonato Tkinter per quasi un anno, è venuto il momento di riprenderne l'uso e di scrivere qualche applicazione utile per il futuro utilizzo nel software astrologico autoprodotto.

Se curiosate tra i widget di Tkinter, anche nella ricca collezione dell'estensione ttk, non trovate un widget calendario. Qualcosa si trova sul web, ma ho pensato che, a questo punto, mi conviene scrivermene uno da solo.

In github trovate il codice di quello fatto da me, è un vero calendario perpetuo, tarato sull'intervallo di anni -3000 : 3000, utile per i nostri calcoli astrologici, facilmente incorporabile all'interno di altri programmi grafici in Python. Di seguito allego un'immagine del widget:

Per spiegarne l'utilizzo dividero' il sorgente in alcuni blocchi che commentero' brevemente. Nei post del 2014 trovate molte utili spiegazioni ulteriori, che possono aiutare a capire meglio come si lavora graficamente con il linguaggio python.

Iniziamo con le prime e le ultime righe di codice:

#/usr/bin/env python
# -*- coding: utf-8 -*-
from Tkinter import *


class tkCalendar:
    import datetime

.............................

if __name__ == '__main__':

    root = Tk()
    cal  = tkCalendar(root)

    root.mainloop()

Iniziamo importando la libreria Tkinter. In fondo al modulo, per testare il widget, avviamo un'istanza della classe tkCalendar e lanciamo il metodo mainloop sull'istanza, che serve a mettere il widget in ciclo continuo, in attesa di eventi.

Subito sotto la dichiarazione di classe, iniziamo a costruire i widget che compongono il calendario.

    def __init__(self, master):
        "this section initializes the graphic objects"
        # main window
        self.Calendar = master
        # top frame (date)
        self.date_labels = Frame(self.Calendar)
        self.date_labels.pack(fill = BOTH, expand = 1)
        for i in ('Year', 'Month', 'Hour', 'Minute'):
            label = Label(self.date_labels, text=i, bg = 'tomato')
            label.pack(side = LEFT, fill = X, expand = 1) 
        
        self.dtime = self.datetime.datetime.now()
        self.date = Frame(self.Calendar)
        self.date.pack(side = 'top', fill = X, expand = 1)

Come abbiamo visto nei post dedicati a Tkinter, la creazione dei widget è abbastanza semplice: si definiscono delle variabili, cui si associa la creazione degli oggetti grafici secondo lo schema:

<variabile> = <Widget>(self.<Nome dell'istanza>, <parametro1>=<valore1>, <parametro2>=<valore2>, ...) 

Dopo la dichiarazione di classe e l'importazione delle librerie che servono, si apre il costruttore con

     def __init__(self, master)

master è un termine convenzionale (potete usare quello che preferite) quindi si crea una finestra contenitore, cui si dà un nome preceduto da self, in questo caso ho usato self.Calendar, che acquisisce il riferimento all'istanza dal parametro master.

Subito dopo creo un frame, cioè un contenitore o, se preferite, una cornice, in cui inserisco quattro label che contengono le etichette in chiaro dell'anno, del mese, dell'ora e del minuto. Per i giorni uso un altro sistema, che vediamo nel seguito. Da notare che perchè i widget siano visibili, devono chiamare il metodo pack() - o grid() o place(). integrati nel pack() ci sono i parametri fill e expand, che determinano il comportamento delle label in caso di modifica della grandezza della finestra generale. Se preferite modificare la direzione di allungamento o di espansione, o i colori che ho preimpostato, potete liberamente intervenire sul codice.

Le ultime due righe del codice sopra riportato creano un frame che conterrà dei widget particolari, chiamati Spinbox.

        # spinbox for year
        self.spin_year = Spinbox(self.date, from_ = -3000, to = 3000, width = 3)
        self.spin_year.pack(side = 'left', fill = X, expand = 1)
        self.spin_year.config(state = 'readonly')
        # spinbox for month
        self.spin_month = Spinbox(self.date,
                                  from_ = 1, to = 12, 
                                  width = 3)
        self.spin_month.pack(side = 'left', fill = X, expand = 1)
        self.spin_month.config(state = 'readonly')# spinbox for hour
        self.spin_hour = Spinbox(self.date,
                                 from_ = 0, to = 23,
                                 width = 3)
        self.spin_hour.pack(side = LEFT, fill = X, expand = 1)
        self.spin_hour.config(state = 'readonly')
        # spinbox for minute
        self.spin_minute = Spinbox(self.date,
                                 from_ = 0, to = 59,
                                 width = 3)
        self.spin_minute.pack(side = LEFT, fill = X, expand = 1)
        self.spin_minute.config(state = 'readonly')

Gli spinbox sono widget molto comodi per selezionare un range ben definito di valori, in questo caso per l'anno fra 3000 e 3000, il mese tra 1 e 12, l'ora tra 0 e 23, il minuto tra 0 e 59. La scelta si fa con le frecce in alto e in basso. Questo puo' risultare scomodo per la selezione di anni estremi nell'intervallo, ma ho preferito disabilitare l'accesso diretto con tastiera per evitare errori grossolani nella gestione del dato.

Arriviamo alla gestione del giorno: i normali calendari di windows o linux consentono la selezione del giorno direttamente sul calendario, con un click del mouse. Ho cercato di realizzare questo comportamento, riproducendo una tabella di valori con i giorni del mese. Il codice seguente, dopo avere creato un frame contenente le labels per i sette giorni della settimana, partendo dalla domenica, realizza appunto una table in cui per ogni giorno del mese, in corripondenza della label indicante il giorno della settimana, si puo' selezionare uno specifico giorno con il click sinistro del mouse.

        # frame for days of the week labels
        self.dotw_labels= Frame(self.Calendar)
        self.dotw_labels.pack(fill = BOTH, expand = 1)
        for i in ('Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat'):
            label = Label(self.dotw_labels, text=i, bg = 'light green',
                          relief = GROOVE)
            label.pack(side = LEFT, fill = BOTH, expand = 1) 
        # table (list of buttons) for days
        self.table = []
        for row in range(0,6):
            # here 6 horizontal frames are created
            # each one receives a block of 7 buttons
            self.row = Frame(self.Calendar)
            self.row.pack(side = TOP, fill = BOTH, expand = 1)
            for col in range(0,7):
                button = Button(self.row,
                                width = 2,
                                bg = 'grey',
                                text = ''
                                )
                # just a little trick: because of lambda usage the callback function
                # receives an updated parameter, referred to the specific button
                # widget, not the one available during
                # initialization
                button.config(command = lambda button = button: self.change_day(button))
                button.pack(side = LEFT, fill = BOTH, expand = 1)
                self.table.append(button)
        self.time_lbl = Label(self.Calendar, text='')
        self.time_lbl.pack()

Per realizzare la tabella dei giorni ho nuovamente usato il trucco che avevo già sperimentato in precedenza: la variabile self.table è una lista, cui associo 6 frame in verticale, abbastanza per contenere i 31 giorni del mese più lungo quando il mese inizia per sabato. Ogni frame viene riempito con sette button, allineati in modo tale da corrispondere alle etichette dei giorni della settimana.

Faccio notare solo il trucchetto usato per recuperare i 'quanti' del giorno direttamente dal testo del button: l'uso della funzione inline lambda consente di lanciare una funzione di callback (self.change_day) con il valore attuale del button selezionato, anzichè quello (vuoto) disponibile al momento della creazione del widget stesso. Per i dettagli sull'uso di queste funzioni potete consultare l'eccellente volume di Mark Lutz "Programming in Python - 4th edition" che dedica ampi capitoli a Tkinter, oppure scrivetemi o inserite un commento a questo post.

Siccome il post è diventato lunghissimo, continuo il commento al codice nel post successivo.

Wednesday, February 18, 2015

Metodi di calcolo delle posizioni planetarie e verifiche di accuratezza e precisione

Con il precedente post ho mostrato come si puo' ottenere un metodo di calcolo delle posizioni planetarie basato sulle VSOP87. Rimane da fare la verifica di quanto questa metodologia sia accurata in relazione ad altri metodi. Se i risultati sono soddisfacenti, possiamo ritenere il metodo utile per le finalità di calcolo astrologico che sono l'obiettivo fondamentale di questo lavoro.

Nel frattempo ho inserito anche il Sole tra i corpi celesti misurabili e ho modificato opportunamente il programma planets.py nella directory vsop87c rinominandolo planets_sun.py. Come già detto trovate tutto in github/zdomjus60/astrometry. Nella directory principale ho inserito uno script per mettere a confronto tre diverse librerie in linguaggio Python: pyephem di Brandon Rhodes, la mia libreria vsop87c e le pyswisseph di Stanislas Marquis.

Riporto di seguito il testo integrale dello script di verifica, di cui trovate copia in github sotto il nome di test_vsop87.py, nella directory principale.

from time_func import *
import math
import random
import ephem
import numpy
import time, datetime
import swisseph
from vsop87c.planets import Planet 
planets = {'Earth': ephem.Sun, 'Mercury':ephem.Mercury, 'Venus':ephem.Venus,
           'Mars':ephem.Mars, 'Jupiter':ephem.Jupiter, 'Saturn':ephem.Saturn,
           'Uranus':ephem.Uranus, 'Neptune':ephem.Neptune}
swe_plan = {'Earth': 0, 'Mercury':2, 'Venus':3,
           'Mars':4, 'Jupiter':5, 'Saturn':6,
           'Uranus':7, 'Neptune':8}
c = []; d = []; e = []
for n in range(1000):
    year = random.randrange(-3000,3000)
    doty = random.randrange(1,365)
    date = cal2jul(year,1,1)+doty
    cal_date = jul2cal(date)
    year = cal_date[0]
    month = cal_date[1]
    day = cal_date[2]
#    print "\njulian date: %d calendar_date %s" % (date , str(cal_date))
#    print "%20s %12s %12s %12s %12s %12s %12s" % ('body','ephem','vsop87c','swisseph','diff1-2','diff2-3', 'diff1-3') 
    for body in planets:
        body_selected = planets[body](str(year)+"/"+str(month)+"/"+str(day), epoch = str(year)+"/"+str(month)+"/"+str(day))
        pos1 = math.degrees(ephem.Ecliptic(body_selected).lon)
        vsop87c_body = Planet(body, year, month, day)
        pos2 = vsop87c_body.calc()[1]
        diff12 = (pos2 - pos1)*3600
        pos3 = swisseph.calc(date, swe_plan[body])[0]
        diff23 = (pos3 - pos2)*3600
        diff13 = (pos3 - pos1)*3600
#        print "%20s %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f" % (body, pos1, pos2, pos3, diff12, diff23, diff13)
        c.append(diff12)
        d.append(diff23)
        e.append(diff13)
print "\n" * 2
print "mean(stdev) ephem-vsop87c      %12.4f(%12.4f)" % (numpy.mean(c), numpy.std(c))
print "mean(stdev) vsop87c-swisseph   %12.4f(%12.4f)" % (numpy.mean(d), numpy.std(d))
print "mean(stdev) ephem-swisseph     %12.4f(%12.4f)" % (numpy.mean(e), numpy.std(e))

Il test procede con l'estrazione casuale di 1000 date tra gli anni 3000 a.c. e 3000 d.c. e con il calcolo di posizione in longitudine di 8 corpi celesti, con l'esclusione di Luna e Plutone. Le differenze di risultato tra i tre metodi sono riportate in linea con le posizioni planetarie e al termine del computo sono calcolate medie e deviazioni standard delle differenze

Noterete che le righe di stampa dei risultati sono commentate con #, per escluderle dall'esecuzione. Questo perchè, ai fini della verifica, la stampa di 1000*8 risultati puo' risultare eccessivamente lenta. Al termine vengono presentate le statistiche relative all'elaborazione, distintamente per le differenze tra coppie di metodi

Un estratto di stampa completo delle longitudini è riportato di seguito:

julian date: 2024877 calendar_date (831, 10, 25, 0, 0, 0)
                body        ephem      vsop87c     swisseph      diff1-2      diff2-3      diff1-3
             Mercury     207.3747     207.3312     207.3156    -156.6550     -56.0983    -212.7533
             Jupiter     205.3468     205.3423     205.3313     -16.1494     -39.5695     -55.7188
              Uranus     329.9765     329.9781     329.9698       5.5322     -29.8444     -24.3122
                Mars     273.3969     273.3770     273.3673     -71.5360     -35.0396    -106.5756
               Earth     215.4045     215.3728     215.3642    -114.1573     -30.7234    -144.8807
               Venus     221.3751     221.3423     221.3272    -118.2600     -54.3486    -172.6085
              Saturn     167.9982     167.9971     167.9892      -4.0064     -28.2423     -32.2487
             Neptune     254.1742     254.1743     254.1742       0.4144      -0.2400       0.1744

julian date: 978435 calendar_date (-2034, 10, 24, 0, 0, 0)
                body        ephem      vsop87c     swisseph      diff1-2      diff2-3      diff1-3
             Mercury     209.0923     213.0602     213.0497   14284.4955     -37.6671   14246.8284
             Jupiter       9.8841     330.6850     330.6922 1154883.3028      25.9853 1154909.2881
              Uranus     261.2739     257.3849     257.3680  -14000.5191     -60.9418  -14061.4609
                Mars      35.8908     180.2785     180.2819  519795.8213      11.9315  519807.7528
               Earth     194.0403     193.7273     193.7197   -1126.8160     -27.3218   -1154.1377
               Venus     152.5836     185.1803     185.1658  117347.8506     -51.8912  117295.9594
              Saturn      42.1960      26.6255      26.6088  -56053.7974     -60.0905  -56113.8879
             Neptune      81.9239      79.7036      79.7322   -7993.3666     103.2981   -7890.0685

julian date: 2730409 calendar_date (2763, 7, 5, 0, 0, 0)
                body        ephem      vsop87c     swisseph      diff1-2      diff2-3      diff1-3
             Mercury     113.3284     113.3347     113.3312      22.7658     -12.7297      10.0361
             Jupiter     169.3535     169.3514     169.3422      -7.6423     -32.8086     -40.4509
              Uranus     359.6431     359.6444     359.6370       4.5214     -26.5551     -22.0336
                Mars      24.3008      24.2813      24.2726     -70.0481     -31.3619    -101.4100
               Earth     103.0054     102.9741     102.9637    -112.4822     -37.6763    -150.1585
               Venus     143.7745     143.7400     143.7283    -124.2814     -42.2936    -166.5750
              Saturn      30.6450      30.6448      30.6358      -0.8895     -32.2358     -33.1253
             Neptune     182.0906     182.0913     182.0795       2.4104     -42.4343     -40.0239

julian date: 2294341 calendar_date (1569, 7, 26, 0, 0, 0)
                body        ephem      vsop87c     swisseph      diff1-2      diff2-3      diff1-3
             Mercury     159.0758     159.0753     159.0686      -1.9116     -23.8803     -25.7919
             Jupiter     273.5696     273.5723     273.5738       9.5099       5.3463      14.8562
              Uranus     266.0435     266.0449     266.0466       4.8892       6.0567      10.9460
                Mars      87.9991      88.0019      87.9939       9.9529     -28.8656     -18.9127
               Earth     132.0953     132.0938     132.0879      -5.6827     -21.1337     -26.8164
               Venus     124.1023     124.0967     124.0973     -20.2555       2.2533     -18.0022
              Saturn     185.7830     185.7847     185.7794       5.9762     -18.9236     -12.9475
             Neptune      80.6680      80.6690      80.6667       3.6657      -8.1534      -4.4877

E infine una delle possibili statistiche dopo mille estrazioni casuali:

mean(stdev) ephem-vsop87c        -1672.3936( 231441.6349)
mean(stdev) vsop87c-swisseph        -5.0075(     84.1484)
mean(stdev) ephem-swisseph       -1677.4011( 231440.8876)

Per quanto ci riguarda, è interessante l'accuratezza relativa del metodo vsop87c rispetto al metodo swisseph.

La libreria pyephem fornisce risultati un po' differenti, sia pure contenuti (per i secoli più vicini al nostro tempo). D'altronde tutti i metodi forniscono un'accuratezza a priori che è relativa a intervalli di date ben definite. Quindi non siamo in alcun modo autorizzati a ritenere un metodo più o meno valido dell'altro, tanto più che c'è, nella rappresentazione delle posizioni planetarie, da considerare il sistema di riferimento e l'apparenza della posizione rispetto alla geometria pura e semplice. Per l'astronomo che deve puntare un telescopio è importante sapere dove trovare il pianeta, la cui posizione apparente risente del tempo di percorrenza della luce nello spazio e di vari fattori modificanti, relativistici e di aberrazione luminosa. L'astrologo normalmente non guarda il cielo ma i suoi grafici astrologici, per cui è interessato a collocare le longitudini dei pianeti dove si presume che le leggi del moto li portino fisicamente. Lo scarto medio di soli 5 arcsecondi con elevata precisione tra le longitudini calcolate con la nostra libreria VSOP87C rispetto alle swiss ephemerides è incoraggiante per le finalità che ci siamo proposti. Rimane il problema plutone, che affronteremo in un post successivo

Friday, February 13, 2015

VSOP87 - metodi di alta precisione per la posizione dei pianeti

Se l'arsenale delle librerie di Python è ricco di soluzioni per i problemi più vari, quello dei metodi matematici per stabilire moto e posizione dei corpi celesti non è da meno. Oggi voglio intrattenervi sui metodi semianalitici sviluppati dal Bureau des Longitudes - Paris - France, che vanno sotto il nome di Variations Séculaires des Orbites Planétaires (VSOP) di cui esistono fondamentalmente due versioni, le VSOP82 (1982 è l'anno di pubblicazione di questo lavoro a cura di Pierre Bretagnon), che fornivano un metodo di grande accuratezza per la definizione dei parametri orbitali, considerati gli effetti di perturbazione legati principalmente alle masse planetarie e alla loro reciproca interazione, e le VSOP87.

Le VSOP87, pubblicate cinque anni dopo le VSOP82, forniscono non solo i correttivi per i parametri orbitali ma permettono di ottenere direttamente le coordinate eliocentriche dei pianeti da Mercurio a Nettuno nella forma rettangolare (X, Y, Z) e sferica (longitudine e latitudine eliocentrica). Come abbiamo visto usando il metodo di Paul Schlyter, il passaggio dalle coordinate eliocentriche a quelle geocentriche è abbastanza immediato, richiedendo poche linee di codice.

Non trovando nulla di pronto in linguaggio Python ho effettuato io stesso la conversione da C a Python con l'uso dell'editor IDLE e della funzione di ricerca e sostituzione, sia alfabetica diretta che con l'uso delle regular expressions. Le formule in linguaggio C sono state ottenute usando il Multi-Language VSOP87 Source Code Generator Tool, straordinario lavoro di Jay Tanner - PHP Science Labs che mette a disposizione degli sviluppatori l'intera serie di correzioni (da poche unità a oltre 1000) previste per le coordinate eliocentriche in codice BASIC, CPP / C++, Java, PHP e VB.NET. Python non è tra i linguaggi considerati ma la conversione è molto semplice, basta eliminare un po' di parentesi graffe, punti e virgola e altri ammenicoli tipici dei linguaggi C-like e riformattare il tutto con un po' di attenzione e rispetto della PEP 8. Per praticità ho utilizzato la tabella C delle VSOP87 che fornisce le coordinate rettangolari (X, Y, Z) per l'equinozio della data

Con le VSOP87 possiamo raggiungere un'accuratezza elevata che risulta stabile per migliaia di anni intorno alla J2000. L'unico pianeta di interesse astrologico non considerato dalle VSOP87 è Plutone, per cui useremo un altro metodo in un post successivo.

La dimensione del codice è tale che non posso includerlo direttamente nel post ma ne faccio un breve estratto. L'integrale, a disposizione di chiunque ritenga di volerlo utilizzare, è in github, come tutto il codice che umilmente sviluppo nel corso di queste conversazioni.

Credo pero' che sia utile pubblicare un piccolo estratto del codice, giusto per mostrare come funziona:

import math
class Mercury:
    
    """

       MERCURY - VSOP87 Series Version C
       HELIOCENTRIC DYNAMICAL ECLIPTIC AND EQUINOX OF THE DATE
       Rectangular (X,Y,Z) Coordinates in AU (Astronomical Units)

       Series Validity Span: 2000 BC < Date < 6000 AD
       Theoretical accuracy over span: +-1 arc sec

       R*R = X*X + Y*Y + Z*Z

       t = (JD - 2451545) / 365250


       C++ Programming Language

       VSOP87 Functions Source Code
       Generated By The VSOP87 Source Code Generator Tool
       (c) Jay Tanner 2015

       Ref:
       Planetary Theories in Rectangular and Spherical Variables
       VSOP87 Solutions
       Pierre Bretagnon, Gerard Francou
       Journal of Astronomy & Astrophysics
       vol. 202, p309-p315
       1988

       Source code provided under the provisions of the
       GNU General Public License (GPL), version 3.
       http://www.gnu.org/licenses/gpl.html

    """
        
    def __init__(self, t):
        self.t = t

    def calculate(self):

        # Mercury_X0 (t) // 1853 terms of order 0

        X0 = 0
        X0 += 0.37749277893 * math.cos(4.40259139579 + 26088.1469590577 * self.t)
        X0 += 0.11918926148 * math.cos(4.49027758439 + 0.2438174835 * self.t)
        X0 += 0.03840153904 * math.cos(1.17015646101 + 52176.0501006319 * self.t)
        X0 += 0.00585979278 * math.cos(4.22090402969 + 78263.95324220609 * self.t)
        X0 += 0.00305833424 * math.cos(2.10298673336 + 26087.65932409069 * self.t)
        X0 += 0.00105974941 * math.cos(0.98846517420 + 104351.85638378029 * self.t)
        X0 += 0.00024906132 * math.cos(5.26305668971 + 52175.56246566489 * self.t)

.........

        # Mercury_X1 (t) // 1023 terms of order 1

        X1 = 0
        X1 += 0.00328639517 * math.cos(6.04028758995 + 0.2438174835 * self.t)
        X1 += 0.00106107047 * math.cos(5.91538469937 + 52176.0501006319 * self.t)
        X1 += 0.00032448440 * math.cos(2.68404164136 + 78263.95324220609 * self.t)
        X1 += 0.00009699418 * math.cos(5.42935843059 + 26087.65932409069 * self.t)

..........

        # Mercury_X2 (t) // 413 terms of order 2

        X2 = 0
        X2 += 0.00020000263 * math.cos(5.96893489541 + 26088.1469590577 * self.t)
        X2 += 0.00008268782 * math.cos(0.41593027178 + 0.2438174835 * self.t)

..........

        # Mercury_X3 (t) // 135 terms of order 3

        X3 = 0
        X3 += 0.00000180660 * math.cos(1.53524873089 + 0.2438174835 * self.t)
        X3 += 0.00000056109 * math.cos(1.50368825372 + 52176.0501006319 * self.t)
        X3 += 0.00000023909 * math.cos(5.09285389011 + 78263.95324220609 * self.t)
        X3 += 0.00000011317 * math.cos(2.22514154728 + 104351.85638378029 * self.t)
        X3 += 0.00000011202 * math.cos(6.20552374893 + 26088.1469590577 * self.t)

..........

        # Mercury_X4 (t) // 42 terms of order 4

        X4 = 0
        X4 += 0.00000043303 * math.cos(2.70854317703 + 26088.1469590577 * self.t)
        X4 += 0.00000016746 * math.cos(2.85109602051 + 0.2438174835 * self.t)
        X4 += 0.00000005097 * math.cos(5.82035608585 + 52176.0501006319 * self.t)
        X4 += 0.00000001110 * math.cos(2.85416500039 + 78263.95324220609 * self.t)


..........

        # Mercury_X5 (t) // 16 terms of order 5

        X5 = 0
        X5 += 0.00000000414 * math.cos(4.09017660105 + 0.2438174835 * self.t)
        X5 += 0.00000000327 * math.cos(2.83894329980 + 26088.1469590577 * self.t)
        X5 += 0.00000000134 * math.cos(4.51536199764 + 52176.0501006319 * self.t)
        X5 += 0.00000000046 * math.cos(1.23267980717 + 78263.95324220609 * self.t)
        X5 += 0.00000000016 * math.cos(4.45794773259 + 104351.85638378029 * self.t)

.........


        X = (X0+
            X1*self.t+
            X2*self.t*self.t+
            X3*self.t*self.t*self.t+
            X4*self.t*self.t*self.t*self.t+
            X5*self.t*self.t*self.t*self.t*self.t)

.........

        Z = ( Z0+
            Z1*self.t+
            Z2*self.t*self.t+
            Z3*self.t*self.t*self.t+
            Z4*self.t*self.t*self.t*self.t+
            Z5*self.t*self.t*self.t*self.t*self.t)

        return (X, Y, Z)

Come si vede dalle ultime righe del codice, i termini in coseno vengono sommati serie per serie, quindi moltiplicati per la distanza temporale in millenni dalla J2000 alla data considerata, elevata a potenza dalla prima alla quinta per ogni signola serie. Lo stesso procedimento si usa per le variabili Y e Z. Il programma restituisce alla fine una tupla contenente le tre coordinate rettangolari per il pianeta e la data considerata.

Per ragioni di ordine e praticità ho inserito i moduli planetari in una directory chiamata vsop87c che trovate in github. la richiesta di coordinate avviene attraverso il modulo planets.py, dentro la stessa directory, che ha lo stesso nome di quello già utilizzato ma codice molto diverso. Di questo riporto l'integrale di seguito:

import math
from time_fn import *
from trigon import *
from earth import Earth
from mercury import Mercury
from venus import Venus
from mars import Mars
from jupiter import Jupiter
from saturn import Saturn
from uranus import Uranus
from neptune import Neptune

class Planet:

    def __init__(self, planet, year, month, day, hour=0, minute =0, second=0):
        self.planet = planet
        self.year = year
        self.month = month
        self.day = day
        self.hour = hour
        self.minute = minute
        self.second = second
        self.planet_list = {'Earth':Earth, 'Mercury':Mercury, 'Venus':Venus, 'Mars':Mars,
                            'Jupiter':Jupiter, 'Saturn':Saturn, 'Uranus':Uranus, 'Neptune':Neptune}


    def calc(self):
        t = julian_millennia(self.year, self.month, self.day,
                             self.hour, self.minute, self.second)
        if self.planet in self.planet_list:
            body = self.planet_list[self.planet](t).calculate()
            earth = self.planet_list['Earth'](t).calculate()
            x = body[0]-earth[0]
            y = body[1]-earth[1]
            z = body[2]-earth[2]
            r = math.sqrt(x*x + y*y + z*z)
            longitude = atan2(y,x) % 360 
            latitude = atan2(z, math.sqrt(x*x + y*y))
            obl_ecl = obl_ecl_Laskar(self.year, self.month, self.day, self.hour, self.minute, self.second)
            f0 = sin(obl_ecl)*sin(longitude)*cos(latitude) + cos(obl_ecl)*sin(latitude)
            f1 = cos(longitude)*cos(latitude)
            f2 = cos(obl_ecl)*sin(longitude)*cos(latitude) - sin(obl_ecl)*sin(latitude)
            RA = atan2(f2,f1) % 360
            r0 = math.sqrt(f1*f1 + f2*f2)
            Decl = atan2(f0,r0)
            return (r, longitude, latitude, RA, Decl)
        
        
if __name__ == '__main__':
    for i in ('Mercury', 'Venus', 'Mars', 'Jupiter',
              'Saturn', 'Uranus', 'Neptune'):
        year, month, day = 2011,4,19
        planet = Planet(i, year, month, day).calc()
        print i, ddd2dms(planet[3]/15), ddd2dms(planet[4])

Come si vede dal codice sopra riportato, è il modulo planets.py che richiama le routine relative ad ogni singolo pianeta, evitando così una gestione eccessivamente polverizzata delle chiamate di classe. Alcune funzioni helper sono contenute in un modulo a parte, time_fn.py, che trovate sempre in github

Un piccolo sfizio da programmatori è l'uso di un dizionario per raccogliere ordinatamente i nomi delle classi, in modo da poter utilizzare i nomi dei pianeti (in Inglese e con l'iniziale maiuscola) per creare l'istanza di classe relativa, nella funzione calc.

Nel prossimo post faremo una piccola verifica di accuratezza, come abbiamo fatto precedentemente per la posizione lunare.

Tuesday, February 3, 2015

Nuove verifiche di accuratezza. Moto lunare

Ora disponiamo di due metodi distinti per calcolare la posizione della Luna, e abbiamo bisogno di sapere quale fiducia possiamo riporre nell'uno e nell'altro. Il metodo di Paul Schlyter è dichiaratamente approssimato, mentre il metodo di Jean Meeus viene indicato come affetto da un errore non superiore a 10 secondi d'arco in longitudine. Proviamo allora a scrivere una breve routine per stimare l'accuratezza dei due metodi. Per semplicità usero' le swiss ephemerides nella modalità Moshier-Ephemeris, cioè senza le effemeridi già tabulate in file di supporto (JPL o Swiss), ma utilizzando un approccio semi-analitico sviluppato da Steve Moshier.

Nel listato che segue è contenuta la routine che chiama sequenzialmente il programma planets.py, quindi il programma moon_meeus.py e infine le swiss ephemerides per calcolare la longitudine lunare per i 31 giorni del mese di marzo 1960, calcola la differenza tra i risultati con i tre metodi, li converte in secondi d'arco e li tabula:

import swisseph
from planets import *
from time_func import *
from moon_meeus import moon_meeus

year = 1960
month = 3
print "Date(Y/M/A)\tlong_planets\tlong_moon_meeus\tlong_swisseph\t   diff1:3\t   diff2:3"
for day in range(1,32):
    g1 = Planet('Moon', year, month, day).position()[1]
    g2 = moon_meeus(year, month, day)[1]
    g3 = swisseph.calc(swisseph._julday(year, month, day, 0, 0, 0),1)[0]
    d13 = (g3-g1)*3600
    d23 = (g3-g2)*3600
    print "%4d %2d %2d %15.4f %15.4f %15.4f %15.4f %15.4f" % (year, month, day, g1, g2, g3, d13, d23)
    

        

Questa è la stampa che si ottiene lanciando la routine:

Date(Y/M/A) long_planets long_moon_meeus long_swisseph    diff1:3    diff2:3
1960  3  1         20.2982         20.3317         20.3322        122.4425          1.7232
1960  3  2         32.9255         32.9644         32.9643        139.8962         -0.3071
1960  3  3         45.2724         45.3107         45.3099        134.9965         -2.9113
1960  3  4         57.3999         57.4302         57.4288        104.1592         -4.7730
1960  3  5         69.3775         69.3939         69.3925         53.6586         -5.1180
1960  3  6         81.2786         81.2793         81.2781         -1.7273         -4.3011
1960  3  7         93.1769         93.1654         93.1645        -44.7431         -3.3203
1960  3  8        105.1452        105.1283        105.1275        -63.6183         -2.8243
1960  3  9        117.2526        117.2372        117.2365        -57.8233         -2.6135
1960  3 10        129.5611        129.5510        129.5505        -38.3935         -2.0565
1960  3 11        142.1210        142.1150        142.1148        -22.3835         -0.9185
1960  3 12        154.9648        154.9580        154.9581        -24.1414          0.2718
1960  3 13        168.1039        168.0904        168.0906        -47.8573          0.7046
1960  3 14        181.5274        181.5036        181.5037        -85.2734          0.1422
1960  3 15        195.2044        195.1714        195.1712       -119.7482         -0.7517
1960  3 16        209.0904        209.0533        209.0531       -134.3417         -1.0077
1960  3 17        223.1333        223.1002        223.1001       -119.2792         -0.3058
1960  3 18        237.2800        237.2590        237.2592        -74.7989          0.7129
1960  3 19        251.4806        251.4778        251.4782         -8.8668          1.1368
1960  3 20        265.6901        265.7086        265.7088         67.1196          0.7483
1960  3 21        279.8692        279.9083        279.9084        141.0218          0.2109
1960  3 22        293.9830        294.0386        294.0387        200.5008          0.2882
1960  3 23        307.9996        308.0644        308.0646        234.0084          0.9793
1960  3 24        321.8873        321.9520        321.9524        234.3308          1.5220
1960  3 25        335.6130        335.6691        335.6694        203.1741          1.2192
1960  3 26        349.1427        349.1850        349.1851        152.7543          0.1567
1960  3 27          2.4447          2.4732          2.4730        101.7851         -0.9317
1960  3 28         15.4954         15.5145         15.5141         67.2924         -1.4241
1960  3 29         28.2846         28.3008         28.3004         56.7443         -1.3436
1960  3 30         40.8192         40.8375         40.8372         64.8680         -1.1482
1960  3 31         53.1239         53.1455         53.1452         76.7926         -1.1081

Anche se le Swiss Ephemerides non sono il metodo di riferimento per la posizione dei corpi celesti, sono tuttavia accreditate come metodo molto accurato, in particolare se usano i file di effemeridi.

Le differenze in longitudine sono espresse in secondi d'arco. L'ultima colonna rappresenta la differenza in longitudine tra i risultati del metodo Meeus rispetto alle Swiss Ephemerides. Mi sembra di poter dire che l'accuratezza è notevole, per cui, nella costruzione del nuovo software, per la posizione della Luna userò certamente il metodo Meeus.

Ricordo che i file che sviluppo si trovano tutti in github. La routine che tabula i dati si chiama test_moon.py.

Con i prossimi post verificheremo la possibilità di migliorare l'accuratezza della posizione del Sole e dei pianeti utilizzando altri metodi. In particolare, nel prossimo post vedremo la routine suggerita da E. M. Standish del Solar System Dinamic Group - JPL Caltech.

Posizione della Luna. Metodo di Jean Meeus

Il metodo di calcolo della posizione della Luna, che presento in questo post, è tratto da Jean Meeus - Astronomical Algorithms, già citato in un post precedente. Questo metodo, una semplificazione dell'ELP2000-82 di Chapront, è basato sul computo di longitudine, latitudine e distanza della Luna utilizzando dei fattori correttivi.

La mia versione replica fedelmente quella riportata da Keith Burnett in forma di foglio elettronico. Per maggiore corrispondenza ho mantenuto quanto più possibile i nomi delle variabili così come vengono chiamati sul foglio elettronico, anche se non seguono le buone prassi indicate dalla PEP 8 - Style Guide for Python Code.

Viste le piccole dimensioni del file, lo trascrivo qui integralmente:

import math
from time_func import *
from trigon import *
def moon_meeus(year, month, day, hour=0, minute =0, second=0):
    coeffs=(
        (0,0,1,0,6288774,-20905355,0,0,0,1,5128122),
        (2,0,-1,0,1274027,-3699111,0,0,1,1,280602),
        (2,0,0,0,658314,-2955968,0,0,1,-1,277693),
        (0,0,2,0,213618,-569925,2,0,0,-1,173237),
        (0,1,0,0,-185116,48888,2,0,-1,1,55413),
        (0,0,0,2,-114332,-3149,2,0,-1,-1,46271),
        (2,0,-2,0,58793,246158,2,0,0,1,32573),
        (2,-1,-1,0,57066,-152138,0,0,2,1,17198),
        (2,0,1,0,53322,-170733,2,0,1,-1,9266),
        (2,-1,0,0,45758,-204586,0,0,2,-1,8822),
        (0,1,-1,0,-40923,-129620,2,-1,0,-1,8216),
        (1,0,0,0,-34720,108743,2,0,-2,-1,4324),
        (0,1,1,0,-30383,104755,2,0,1,1,4200),
        (2,0,0,-2,15327,10321,2,1,0,-1,-3359),
        (0,0,1,2,-12528,0,2,-1,-1,1,2463),
        (0,0,1,-2,10980,79661,2,-1,0,1,2211),
        (4,0,-1,0,10675,-34782,2,-1,-1,-1,2065),
        (0,0,3,0,10034,-23210,0,1,-1,-1,-1870),
        (4,0,-2,0,8548,-21636,4,0,-1,-1,1828),
        (2,1,-1,0,-7888,24208,0,1,0,1,-1794),
        (2,1,0,0,-6766,30824,0,0,0,3,-1749),
        (1,0,-1,0,-5163,-8379,0,1,-1,1,-1565),
        (1,1,0,0,4987,-16675,1,0,0,1,-1491),
        (2,-1,1,0,4036,-12831,0,1,1,1,-1475),
        (2,0,2,0,3994,-10445,0,1,1,-1,-1410),
        (4,0,0,0,3861,-11650,0,1,0,-1,-1344),
        (2,0,-3,0,3665,14403,1,0,0,-1,-1335),
        (0,1,-2,0,-2689,-7003,0,0,3,1,1107),
        (2,0,-1,2,-2602,0,4,0,0,-1,1021),
        (2,-1,-2,0,2390,10056,4,0,-1,1,833),
        (1,0,1,0,-2348,6322,0,0,1,-3,777),
        (2,-2,0,0,2236,-9884,4,0,-2,1,671),
        (0,1,2,0,-2120,5751,2,0,0,-3,607),
        (0,2,0,0,-2069,0,2,0,2,-1,596),
        (2,-2,-1,0,2048,-4950,2,-1,1,-1,491),
        (2,0,1,-2,-1773,4130,2,0,-2,1,-451),
        (2,0,0,2,-1595,0,0,0,3,-1,439),
        (4,-1,-1,0,1215,-3958,2,0,2,1,422),
        (0,0,2,2,-1110,0,2,0,-3,-1,421),
        (3,0,-1,0,-892,3258,2,1,-1,1,-366),
        (2,1,1,0,-810,2616,2,1,0,1,-351),
        (4,-1,-2,0,759,-1897,4,0,0,1,331),
        (0,2,-1,0,-713,-2117,2,-1,1,1,315),
        (2,2,-1,0,-700,2354,2,-2,0,-1,302),
        (2,1,-2,0,691,0,0,0,1,3,-283),
        (2,-1,0,-2,596,0,2,1,1,-1,-229),
        (4,0,1,0,549,-1423,1,1,0,-1,223),
        (0,0,4,0,537,-1117,1,1,0,1,223),
        (4,-1,0,0,520,-1571,0,1,-2,-1,-220),
        (1,0,-2,0,-487,-1739,2,1,-1,-1,-220),
        (2,1,0,-2,-399,0,1,0,1,1,-185),
        (0,0,2,-2,-381,-4421,2,-1,-2,-1,181),
        (1,1,1,0,351,0,0,1,2,1,-177),
        (3,0,-2,0,-340,0,4,0,-2,-1,176),
        (4,0,-3,0,330,0,4,-1,-1,-1,166),
        (2,-1,2,0,327,0,1,0,1,-1,-164),
        (0,2,1,0,-323,1165,4,0,1,-1,132),
        (1,1,-1,0,299,0,1,0,-1,-1,-119),
        (2,0,3,0,294,0,4,-1,0,-1,115),
        (2,0,-1,-2,0,8752,2,-2,0,1,107)
        )
    jd = cal2jul(year, month, day, hour, minute, second)
    epoch = cal2jul(2000,1,1, 12)
    t = (jd - epoch)/36525
    obl_ecl = (84381.448-46.815*t-0.00059*t*t+0.001813*t*t*t)/3600
    L1 = (218.3164591
          + 481267.88134236*t
          - 0.0013268*t*t
          + t*t*t/538841
          - t*t*t*t/65194000) % 360
    D = (297.8502042 +
            445267.1115168*t
          - 0.00163*t*t
          + t*t*t/545868
          - t*t*t*t/113065000) %360
    M = (357.5291092
          + 35999.0502909*t
          - 0.0001536*t*t
          + t*t*t/24490000) % 360
    M1 = (134.9634114
          + 477198.8676313*t
          + 0.008997*t*t
          + t*t*t/69699
          - t*t*t*t/14712000) %360
    F = (93.2720993
          + 483202.0175273*t
          - 0.0034029*t*t
          + t*t*t/3526000
          - t*t*t*t/863310000) % 360
    A1 = (119.75 + 131.849*t) % 360
    A2 = (53.09+479264.29*t) % 360
    A3 = (313.45+481266.484*t) % 360
    E = 1 - 0.002516*t - 0.0000074*t*t
    E2 = E*E


    l_eccen = []
    for i in range(0,60):
        if abs(coeffs[i][1])==1:
            l_eccen.append(coeffs[i][4]*E)
        elif abs(coeffs[i][1])==2:
            l_eccen.append(coeffs[i][4]*E2)
        else:
            l_eccen.append(coeffs[i][4]*1.0)

    r_eccen = []
    for i in range(0,60):
        if abs(coeffs[i][1])==1:
            r_eccen.append(coeffs[i][5]*E)
        elif abs(coeffs[i][1])==2:
            r_eccen.append(coeffs[i][5]*E2)
        else:
            r_eccen.append(coeffs[i][5]*1.0)

    b_eccen = []
    for i in range(0,60):
        if abs(coeffs[i][7])==1:
            b_eccen.append(coeffs[i][10]*E)
        elif abs(coeffs[i][7])==2:
            b_eccen.append(coeffs[i][10]*E2)
        else:
            b_eccen.append(coeffs[i][10]*1.0)

    l_term = []
    for i in range(0,60):
        l_term.append(l_eccen[i]*sin(coeffs[i][0] * D +
                                     coeffs[i][1] * M +
                                     coeffs[i][2] * M1 +
                                     coeffs[i][3] * F
                                     )
        )
    r_term = []
    for i in range(0,60):
        r_term.append(r_eccen[i]*cos(coeffs[i][0] * D +
                                     coeffs[i][1] * M +
                                     coeffs[i][2] * M1 +
                                     coeffs[i][3] * F
                                     )
        )

    b_term = []
    for i in range(0,60):
        b_term.append(b_eccen[i]*sin(coeffs[i][6] * D +
                                     coeffs[i][7] * M +
                                     coeffs[i][8] * M1 +
                                     coeffs[i][9] * F
                                     )
        )



    l_add = sum(l_term) + 3958 * sin(A1) + 1962 * sin(L1-F) + 318 * sin(A2)
    r_add = sum(r_term)
    b_add = ( sum(b_term) - 2235 * sin(L1) + 382 * sin(A3)
              + 175 * sin(A1-F) + 175 * sin(A1+F) + 127 * sin(L1-M1)
              + 115 * sin(L1+M1))

    mean_longitude = L1 + l_add/1000000.0
    mean_latitude  = b_add / 1000000.0
    r = 385000.56 + r_add/1000.0
    
    long_asc_lunar_node = (125.04452-1934.136261*t) % 360
    mean_long_sun = (280.4665+36000.7698*t) % 360
    mean_long_moon = (218.3165+481267.8813*t) % 360

    delta_phi = ( -17.2*sin(long_asc_lunar_node)
                  -1.32*sin(2*mean_long_sun)
                  -0.23*sin(2*mean_long_moon)
                  +0.21*sin(2*long_asc_lunar_node))
    delta_e = ( 9.2*cos(long_asc_lunar_node)
                +0.57*cos(2*mean_long_sun)
                +0.1*cos(2*mean_long_moon)
                -0.09*cos(2*long_asc_lunar_node))

    longitude = mean_longitude + delta_phi / 3600.0
    latitude = mean_latitude + delta_e / 3600.0
    
    x = r * cos(longitude) * cos(latitude)
    y = r * sin(longitude) * cos(latitude)
    z = r                  * sin(latitude)

    x_eq = x
    y_eq = y * cos(obl_ecl) - z * sin(obl_ecl)
    z_eq = y * sin(obl_ecl) + z * cos(obl_ecl)

    RA  = atan2( y_eq, x_eq )
    Dec = atan2( z_eq, math.sqrt(x_eq*x_eq+y_eq*y_eq) )

    return r, reduce360(longitude), latitude, reduce360(RA), Dec

if __name__ == '__main__':
    print moon_meeus(1992,4,12)
    

Come per il file planets.py, la routine restituisce una tupla contenente distanza in Raggi Terrestri, longitudine, latitudine, ascensione retta e declinazione, tutto in gradi.

How to create a virtual linux machine with qemu under Debian or Ubuntu with near native graphics performance

It's been a long time since my latest post. I know, I'm lazy. But every now and then I like to publish something that other people c...