Moto ellittico e determinazione della posizione di un pianeta dati gli elementi orbitali

La seconda area di interesse, nella costruzione del nostro software astrologico, è quella dedicata alla determinazione, quanto più possibile precisa, della posizione di un corpo celeste nello spazio. L'astrologia adotta una rappresentazione estremamente semplificata, di fatto riducendosi a considerare esclusivamente la longitudine geocentrica eclittica. Sarà quindi compito del nostro software determinare questa variabile e, per farlo, utilizzeremo le normali modalità di calcolo astronomico in uso, cercando di contenere l'errore e, quando si renderà utile o necessario, verificando più metodi di calcolo. La precisione richiesta per la determinazione della longitudine eclittica di un pianeta varia secondo le tecniche utilizzate.

Inizieremo con il calcolo della posizione del Sole, passando in rassegna diverse tecniche e cercando di stimare l'errore. Come anticipato in un post precedente, partiamo dalla pagina web di Paul Schlyter.

Dalle leggi di Keplero si ricavano i metodi per la determinazione della posizione di un corpo celeste che si muove in un'orbita ellittica.

Prima di procedere oltre ho deciso di rendere accessibile on line il software che progressivamente produco. Per farlo usero' github. Usate liberamente questo link: https://github/zdomjus60. Il repository di interesse si chiama astrometry. Se sapete usare git potete clonare il repository e lavorarci sopra, altrimenti scaricatelo come file zip. Cerchero' di tenerlo aggiornato in linea con il blog.

Passiamo al moto ellittico. In questo post uso il metodo indicato da Paul Schlyter nella pagina How to compute planetary positions.

Per usare il programma planets.py è sufficiente un'operazione di import come si vede nel file test.py che riporto integralmente di seguito:

# -*- coding: utf-8 -*-
# calculate Sun longitude according to Paul Schlyter formulas
from trigon import *
from planets import Planet
year, month, day = (1990, 4, 19)
for body in ('Sun', 'Moon', 'Mercury', 'Venus', 'Mars',
               'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto'):
    
    planet = Planet(body, year, month, day).position()
    print body, planet

Ho cercato di seguire fedelmente le istruzioni di Paul Schlyter per facilità di confronto del codice e dei risultati di calcolo. Ho lasciato, per esempio, l'ascensione retta in gradi, anche se di solito viene rappresentata in ore (é sufficiente dividere i gradi per 15 per avere la rappresentazione oraria).

La precisione che si raggiunge non è male, ma nel seguito faremo altre prove. Invito sempre chi legge a contribuire, anche su github, mettendo a disposizione qualche suggerimento. Il programma è ancora un po' ruvido e richiede semplificazione e ottimizzazione. Nel seguito faro' qualche aggiornamento. Alla prossima.

Funzioni helper per la gestione dell'ora e del calendario

Nel calcolo astronomico è di primaria importanza la gestione del calendario e dell'ora. Ogni luogo sulla Terra, se basasse il computo del tempo sui fenomeni astronomici quali la rotazione apparente del Sole o delle stelle, avrebbe necessariamente un'ora diversa da tutti gli altri luoghi. Il sorgere e tramonto del Sole avvengono in istanti diversi secondo la longitudine geografica e sono variabili secondo la latitudine. La necessità di far coincidere gli orologi, almeno nell'ambito di una comunità organizzata quale quella di un Paese, ha reso necessaria la standardizzazione dell'ora a scapito dell'allineamento con i fenomeni astronomici reali, introducendo variazioni in eccesso o in difetto per adottare un'ora solare convenzionale. Analogamente, l' introduzione dell'ora legale ha comportato l' artificiosa alterazione del corso delle ore, aggiungendo o togliendo in certi giorni dell'anno un'ora dal corso normale.

A queste cause di alterazione dell'ora si aggiunge l'adozione di un fuso orario che talora rende molto evidente il divario tra l'ora dell'orologio e la posizione del Sole. Per esempio la Spagna usa la stessa ora della Germania quando sarebbe più naturale, per longitudine, l'allineamento all'ora media di Greenwich.

L'uso tradizionale del sistema sessagesimale rispetto a quello decimale per gli angoli e le ore comporta inoltre una difficoltà aggiuntiva.

Per rappresentare l'ora di calcolo della posizione dei corpi celesti in modo unico per ogni luogo della Terra, indipendentemente dagli orologi, occorre considerare almeno le seguenti funzioni:

• Conversione di angoli e ore dal sistema sessagesimale al decimale e viceversa
• Gestione del calendario gregoriano (introdotto nel 1582 per la correzione del ritardo storico tra posizione reale del Sole e la sua rappresentazione nel calendario, ma che ha, di fatto, introdotto una discontinuità, resa ancora più complicata dall' adozione non simultanea del calendario gregoriano stesso nei diversi paesi del mondo). La libreria standard di Python purtroppo non consente di gestire il calendario giuliano, per cui è necessario adottare librerie di terze parti o scrivere da sè le principali funzioni, scelta che preferisco

Inizio quindi con la creazione di un modulo, che chiamero' time_func.py, dove inseriro' alcune routine di gestione del calendario e dell'ora. Allego il listato del codice per commentarlo poi brevemente. Per la stesura delle routine mi sono basato anche su un famoso testo di calcolo astronomico ben noto ai programmatori: Duffett-Smith, Zwart - Practical Astronomy with Your Calculator or Spreadsheet (Cambridge, 4th Ed., 2011)

# -*- coding: utf-8 -*-
""" helper functions for time management
"""
import math

def dms2ddd(hour, minute, second):
    """ from sexagesimal to decimal """
    return hour+minute/60.0+second/3600.0

def ddd2dms(dec_hour):
    """ from decimal to sexagesimal representation of hours and angles."""
    total_seconds = int(dec_hour * 3600+.5)
    seconds = total_seconds % 60
    total_minutes = int((total_seconds - seconds)/60)
    minutes = total_minutes % 60
    hours = int((total_minutes - minutes)/60)
    return (hours, minutes, seconds)

def cal2jul(year, month, day, hour=0, minute=0, second=0):
    """ converts calendar date to julian date
        this routine and the following are built following Duffet Smith /Zwart instructions
        as given in Peter Duffett-Smith- Practical Astronomy with your Calculator or Spreadsheet
        Fourth Edition, Cambridge University Press, Fourth Ed. 2011
        For an easier use of the function, hours minutes and seconds are defaulted to 0, so it's
        not necessary to give them as parameters when the hour is 00:00:00 
    """
    month2 = month
    year2 = year
    if month2 <= 2:
        year2 -= 1
        month2 += 12
    else:
        pass
    if (year*10000 + month*100 + day) > 15821015:
        a = math.trunc(year2/100)
        b = 2 - a + math.trunc(a/4)
    else:
        a = 0
        b = 0
    if year < 0:
        c = math.trunc((365.25 * year2)-0.75)
    else:
        c = math.trunc(365.25 * year2)
    d = math.trunc(30.6001 *(month2 + 1))
    return b + c + d + day + hour / 24.0 + minute / 1440.0 + second / 86400.0 + 1720994.5

def jul2cal(jd):
    """ converts julian date to calendar date """
    jd += 0.5
    i = math.modf(jd)[1]
    f = math.modf(jd)[0]
    if i > 2299160:
        a = math.trunc((i-1867216.25)/36524.25)
        b = i + a - math.trunc(a/4)+1
    else:
        b = i
    c = b + 1524
    d = math.trunc((c-122.1)/365.25)
    e = math.trunc(365.25 * d)
    g = math.trunc((c-e)/30.6001)
    day = c-e+f-math.trunc(30.6001*g)
    if g < 13.5:
        month = g - 1
    else:
        month = g - 13
    if month > 2.5:
        year = d - 4716
    else:
        year = d - 4715
    
    hours_frac = math.modf(day)[0]*24
    day = int(day)
    hour, minute, second = ddd2dms(hours_frac) 
    return (year, month, day, hour, minute, second)

def day_of_the_week(year, month, day):
    """ given a calendar date, the routine returns a tuple with the Day Of The Week in number and in plaintext
        0 for Sunday 1 for Monday and so on up to 6 Saturday
    """
    doth = {0:'Sunday', 1:'Monday', 2:'Tuesday',
            3:'Wednesday', 4:'Thursday', 5:'Friday',
            6:'Saturday'}
    jd = cal2jul(year, month, day, 0, 0, 0)
    a = (jd+1.5)/7
    f = math.trunc((a % 1)*7 +.5)
    return (f,doth[f])

def lt2ut(year, month, day, hour=0, minute=0, second=0, timezone=0, DS=0):
    """ Given, for a location on the Earth,a date, a time, a timezone (East + West - in hours) and the Daylight
        Savings (0 normal time 1 Daylight Savings), this routine gives back a calendar date in Universal Time
        representation (year, month, day, hour, minute, second).
        It aims to restore a common date and time for all places in the Earth. Timezone and
        Daylight Savings can be automized knowing the location using the pytz module (Olson
        database)
    """
    ut = dms2ddd(hour,minute,second) - timezone - DS
    greenwich_calendar_date = day + ut/24
    jd = cal2jul(year, month, greenwich_calendar_date)
    greenwich_calendar_date = jul2cal(jd)
    return greenwich_calendar_date

def ut2lt(year, month, day, hour=0, minute=0, second=0, timezone=0, DS=0):
    """ Given a date, a time for Greenwich in UT format this routine gives back a calendar date
        in local time representation (year, month, day, hour, minute, second).
        It's the inverse function of the previous formula 
    """
    lt = dms2ddd(hour,minute,second) + timezone +DS
    local_calendar_date = day + lt/24
    jd = cal2jul(year, month, local_calendar_date)
    local_calendar_date = jul2cal(jd)
    return local_calendar_date

def ut2gst(year, month, day, hour, minute, second):
    """ Sidereal time is a time-keeping system astronomers use to keep track of the direction to point
        their telescopes to view a given star in the night sky.
        Briefly, sidereal time is a "time scale that is based on the Earth's rate of rotation measured
        relative to the fixed stars." (source Wikipedia)
        This routine converts Universal Time to Sidereal Time for Greenwich (Greenwich Sidereal Time)
    """
    jd = cal2jul(year, month, day)
    S = jd - 2451545.0
    T = S/36525.0
    T0 = (6.697374558 + (2400.051336 * T)+ 0.000025862 *T*T) % 24
    UT = dms2ddd(hour, minute, second)*1.002737909
    GST = ddd2dms((UT + T0) % 24)
    return GST
    
def gst2ut( year, month, day, hour, minute, second):
    """ Inverse of the previous function
    """
    jd = cal2jul(year, month, day, 0,0,0)
    S = jd - 2451545.0
    T = S/36525.0
    T0 = (6.697374558 + 2400.051336 * T + 0.000025862 *T*T) % 24
    GST = (dms2ddd(hour, minute, second) - T0) % 24
    while GST <0:
        GST += 24
    UT = GST * .9972695663
    return ddd2dms(UT)

def gst2lst( hour, minute, second, long_degree, long_minute, long_second=0):
    """ Corrects GST for a different location on the Earth
    """
    GST = dms2ddd(hour,minute,second)
    lg = dms2ddd(long_degree, long_minute, long_second)/15
    lst = ddd2dms((GST + lg) % 24)
    return lst

def lst2gst( hour, minute, second, long_degree, long_minute, long_second=0):
    """ Inverse of the previous method
    """
    lst = dms2ddd(hour,minute,second)
    lg = dms2ddd(long_degree, long_minute, long_second)/15
    GST = ddd2dms((lst + lg) % 24)
    return GST
    
if __name__ == '__main__':
    dash_line = 80 * '-'
    print dash_line
    ddd = dms2ddd(15,24,35)
    print '''from sexagesimal to decimal of 15h 24' 35": ''', ddd
    dms = ddd2dms(ddd)
    print 'and inverse function: ', dms
    print dash_line
    jd = cal2jul(2009,6,19,18,3,21)
    print 'julian date from calendar date 2009/06/19 18:03:21', jd
    cal = jul2cal(jd)
    print 'inverse: calendar date from julian date ', jd, ' :',cal
    print dash_line
    epoch = cal2jul(1999,12,31)
    date_ = cal2jul(1990, 4,19)
    print 'days from 1999/12/31', date_ - epoch
    print dash_line
    doth = day_of_the_week(2009, 6, 19)
    print 'day of the week for 2009/6/19 :',doth
    doth = day_of_the_week(2015, 1, 25)
    print 'day of the week for 2015/1/25 :',doth
    ddd = dms2ddd(19,20,37)
    print dash_line
    print """conversion 19h 20' 37": """, ddd
    print "inverse conversion :", ddd2dms(ddd)
    print dash_line 
    ut = lt2ut(2013,7,1, 3, 37,0, 4, 1)
    print "greenwich calendar date for 2013/7/1 3h37 timezone +4 Daylight Savings)", ut
    lt = ut2lt(ut[0], ut[1], ut[2], ut[3], ut[4], ut[5], 4, 1)
    print "inverse of the previous calculus", lt 
    print dash_line
    print """from Universal Time to Greenwich Sidereal Time for 1980/4/22 14h36'52" """, ut2gst(1980,4,22, 14,36,52)
    print """from Greenwich Sidereal Time to Universal Time for 1980/4/22 4h40'06"  """, gst2ut(1980,4,22, 4,40,06)
    print dash_line
    print """ from GST 4h40'06" to LST for 64°W """, gst2lst(4,40,6,-64,0,0)
    print """ inverse of the previous formula               """, lst2gst(0,24,6,64,0,0)
    

Ho preferito inserire i commenti in lingua inglese per consentire una maggiore diffusione delle routine a chiunque fosse interessato. Dovrebbero risultare sufficientemente chiari anche ai lettori di lingua italiana. Per ogni vostra necessità attendo commenti e suggerimenti. Al prossimo post.

Alcune considerazioni su progettazione del software, moduli e classi

Nel momento in cui si iniziano a definire le linee progettuali di un software, si cerca anche di capire quali strumenti un linguaggio di programmazione mette a disposizione per assistere il programmatore nella costruzione del software stesso.

Gli esempi di codice che ho pubblicato finora sono sufficientemente piccoli da non richiedere una eccessiva strutturazione. Man mano che le dimensioni del progetto aumentano fino a diverse migliaia di linee di codice, se non abbiamo fatto una scelta iniziale corretta degli strumenti, siamo destinati a entrare in crisi. Una buona progettazione puo' proteggerci dal rischio di fallimento iniziale e da quelli futuri, quando si imporrà necessariamente un ammodernamento del software o una sua parziale riscrittura.

La prima cosa che dobbiamo fare è individuare le "aree" in cui il nostro software si muove: per un software astrologico, possiamo tentare di schematizzarle così:

1. Area del tempo, del calendario, delle convenzioni locali
2. Area del movimento dei pianeti, delle posizioni relative degli stessi, dei loro moti apparenti
3. Area dello zodiaco, dell'eclittica e dell'equatore celeste, delle angolarità tra gli elementi del tema
4. Area di ciò che non rientra necessariamente nelle prime tre

Cio' detto, si impone una prima separazione tra queste aree per creare una prima mappa del software e iniziarne lo sviluppo.

Cosa conviene utilizzare per il tempo? Sicuramente delle funzioni del tempo, che abbiano a che fare con la gestione dell'ora e della localizzazione temporale. Forse lo strumento più adatto è un modulo dedicato, fatto di definizioni di funzioni. In python un modulo è un file con estensione .py che si trova in una posizione rintracciabile dal programma principale, di solito nella stessa directory o in una sottodirectory. Ci basterà chiamare il modulo con un import, seguito dal nome del modulo senza estensione .py. Ogni funzione sarà accessibile chiamando il nome del modulo seguito da un punto e dal nome della funzione stessa.

Per i pianeti puo' valere forse la pena di considerarli come oggetti, istanze di un prototipo comune definito attraverso una classe. In questo modo potremo creare tanti oggetti quanti sono i pianeti (e altri componenti del tema) per legare in un insieme unico proprietà e metodi, cioè variabili di classe o di istanza (ne parleremo a breve) e funzioni proprie della classe.

In effetti scrivere una classe ha senso se dobbiamo creare dei prototipi con molte istanze, quindi va bene per oggetti abbastanza elementari ma anche per il tema natale nel suo complesso, qualora volessimo confrontare fra loro più temi (fatti di oggetti più semplici), riferiti a persone diverse o a momenti diversi della vita di una stessa persona.

Dovremo scegliere anche un modo per garantire lo storage definitivo dei dati raccolti, per esempio attraverso un database, che potrebbe essere quindi lo strumento adatto per questa funzione.

Per andare avanti abbiamo bisogno anche di un documento teorico, che consenta di effettuare i calcoli astronomici e/o astrologici con accuratezza e precisione. Per la sua semplicità, relativamente al moto dei pianeti, penso di avvalermi, come già in passato, delle pagine di Paul Schlyter How to compute planetary positions e Computing planetary positions - a tutorial with worked examples.

Da queste iniziamo, a partire dal prossimo post, a costruire un primo nucleo di funzioni dedicate alla misura del tempo e al moto dei pianeti. Ove necessario, attingeremo ad altre fonti. Al prossimo post.

Linguaggio Python. Alcune particolarità interessanti

Riprendo a scrivere il mio blog, anche se con un po' di fatica, dopo molti mesi di assenza. Ho accantonato il progetto di conversione tra i formati SVG e Tk perchè si tratta di un progetto impegnativo e non ho più avuto il tempo di stargli dietro. Lo riprendero' appena avro' di nuovo voglia e tempo di farlo.

Ho invece lasciato un po' in sospeso la parte dedicata ai riferimenti temporali e geografici del l'ora, fondamentali per il corretto calcolo della domificazione, cioè dell'Ascendente e del Medio Cielo. Mi riprometto di continuare quanto già iniziato in precedenza con il database delle località e di parlare un po' del tempo nei suoi aspetti convenzionali, cioè il tempo degli orologi rispetto al tempo terrestre, i fusi orari e i Daylight Savings, che in italiano chiamiamo ora legale, e la libreria tz di Python che consente di accedere velocemente alla storia del tempo nei vari paesi del mondo.

Prima di procedere oltre vorrei accennare ad un problema che ha a che fare con i metodi di computo della posizione dei corpi celesti e i diversi sistemi di coordinate. Pur avendo finora usato le swiss ephemerides nella versione di libreria Python Pyswisseph, credo che sia utile parlare dei metodi astronomici di calcolo e di come si puo' costruire una libreria ex novo.

Le funzioni trigonometriche sono disponibili in Python attraverso la libreria math. Per chiamarla nel nostro codice sarà sufficiente un'operazione di import:

import math

Utilizzando il metodo dir(math) disporremo dell'intero set di funzioni matematiche, alcune delle quali goniometriche, e di due costanti: pi greco e numero e.

>>> dir(math)
['__doc__', '__name__', '__package__', 'acos', 'acosh', 'asin', 'asinh', 'atan', 'atan2', 'atanh', 'ceil', 'copysign', 'cos', 'cosh', 'degrees', 'e', 'erf', 'erfc', 'exp', 'expm1', 'fabs', 'factorial', 'floor', 'fmod', 'frexp', 'fsum', 'gamma', 'hypot', 'isinf', 'isnan', 'ldexp', 'lgamma', 'log', 'log10', 'log1p', 'modf', 'pi', 'pow', 'radians', 'sin', 'sinh', 'sqrt', 'tan', 'tanh', 'trunc']
>>> 

Ogni funzione sarà accessibile richiamando nel codice il package math seguito da un punto e dalla funzione desiderata.

>>> print math.pi
3.14159265359
>>> print math.e
2.71828182846
>>> print math.sin(math.pi/6)
0.5
>>> 

Nel calcolo astronomico si utilizza comunemente il formato gradi/minuti/secondi o gradi e frazioni di grado anzichè i radianti come normalmente usati in trigonometria. La conversione è semplice: 180° corrispondono a pi radianti, per cui un modo per usare nativamente funzioni goniometriche di grado anzichè di radiante consiste nel definire, nel namespace principale o in un modulo, delle nuove funzioni dedicate:

def sin(x):
 return math.sin(math.radians(x))

Quindi se chiamiamo sin(30) anziché math.sin(math.pi/6) otterremo lo stesso risultato ma potremo utilizzare le formule che normalmente usano gli astronomi.

La soluzione proposta è in assoluto la più semplice. Un'alternativa possibile è l'uso dei decoratori:

import math

def deg(func):
    def wrapper(x):
        return func(math.radians(x))
    return wrapper

@deg
def sin(x):
    return math.sin(x)

@deg
def cos(x):
    return math.cos(x)
    
print sin(30)
print cos(30)
>>> 
0.5
0.866025403784
>>> 

In questo secondo caso si usa il decoratore per modificare la variabile argomento prima di applicare la funzione goniometrica, non è l'esempio più felice di uso dei decoratori, ha solo, evidentemente, una funzione dimostrativa.

Terza possibilità: usare i metodi getattr e setattr per generare le nuove funzioni goniometriche nello spazio globale:

import math

def deg(func):
    def wrapper(x):
        return func(math.radians(x))
    return wrapper

for i in ("sin", "cos", "tan"):
    globals()[i]=deg(getattr(math,i))

>>> dir()
['__builtins__', '__doc__', '__file__', '__name__', '__package__', 'cos', 'deg', 'i', 'math', 'sin', 'tan']
>>> globals()
{'cos': <function wrapper at 0x7fbbc1780938>, '__builtins__': <module '__builtin__' (built-in)>, '__file__': '/home/mint/prova.py', '__package__': None, 'i': 'tan', '__name__': '__main__', 'tan': <function wrapper at 0x7fbbc17809b0>, 'sin': <function wrapper at 0x7fbbc17802a8>, '__doc__': None, 'math': <module 'math' (built-in)>, 'deg': <function deg at 0x7fbbc627c6e0>}
>>> sin(30)
0.49999999999999994
>>> cos(30)
0.8660254037844387
>>>   

Va da sè che le possibilità offerte dal linguaggio Python sono tante, ma atteniamoci alla metodologia più semplice, che è quella illustrata per prima. Nei prossimi post inizieremo a scrivere una nuova libreria di calcolo astronomico, in cui proporremo delle formule di calcolo approssimato, ma sufficiente per molti dei nostri scopi.