Monday, February 2, 2015

Verifiche di accuratezza delle posizioni planetarie

Prima di procedere oltre, vorrei fare una prima verifica di accuratezza delle posizioni calcolate con il metodo di Paul Schlyter, utilizzando le effemeridi calcolate via web alla pagina ephemeris.com per la data del 19 aprile 1990 ore 00.00, utilizzando i server NASA/JPL.

Affiancati ai valori ottenuti in tal modo, quelli che ho calcolato con il programma planets.py

Date/Time: 1990.04.19 00:00:00 UTC (GMT - Delta T), JD = 2448000.500000 Sidereal Time: 13:47:20, Delta T = 57.1 seconds Geocentric positions Phase of Moon: 0.773 (0.000=New; 0.250=First 1/4; 0.500=Full; 0.750=Last 1/4) Planet Longitude Latitude Right Asc. Declination longitudine latitudine asc. retta declinazione Sun 28 Ari 40'46" - 0°00'00" 01:46:36 11°00'24" ( 28, 41, 13) ( 0, 0, 0) ( 1, 46, 38) ( 11, 0, 30) Moo 06 Aqr 57'08" - 0°33'08" 20:37:59 -19°04'17" (306, 56, 54) ( 0, -35, - 8) (20, 38, 0) (-19, - 6, -12) Mer 16 Tau 37'08" 2°57'53" 02:53:01 19°38'45" ( 46, 37, 28) ( 2, 57, 52) ( 2, 53, 2) ( 19, 38, 45) Ven 13 Psc 18'20" - 0°21'46" 22:59:01 - 6°53'50" (343, 18, 48) ( 0, -21, -48) (22, 59, 3) (- 6, -53, -39) Mar 28 Aqr 33'42" - 1°25'07" 22:04:53 -13°18'21" (328, 33, 45) (- 1, -25, - 9) (22, 4, 53) (-13, -18, -18) Jup 05 Cnc 02'50" 0°06'16" 06:22:01 23°27'01" ( 95, 2, 31) ( 0, 6, 17) ( 6, 21, 59) ( 23, 26, 57) Sat 25 Cap 07'49" 0°11'01" 19:48:11 -20°55'48" (295, 6, 1) ( 0, 11, 12) (19, 48, 3) (-20, -55, -51) Ura 09 Cap 34'42"R - 0°18'17" 18:41:47 -23°24'01" (279, 34, 17) ( 0, -18, -21) (18, 41, 45) (-23, -24, - 1) Nep 14 Cap 34'22"R 0°51'54" 19:02:53 -21°47'08" (284, 34, 4) ( 0, 51, 54) (19, 2, 51) (-21, -47, - 4) Plu 16 Sco 53'59"R 15°59'26" 15:15:42 - 1°31'06" (226, 53, 31) ( 15, 58, 55) (15, 15, 39) (- 1, -31, -23)

I valori nelle tuple relative a latitudine e declinazione sono negativi anche nei minuti e secondi, per maggior agio nell'utilizzo successivo. La tabella di ephemeris.com riporta il grado di longitudine nel segno zodiacale di appartenenza, come si puo' facilmente intuire.

Pur rilevando una precisione discreta, ai fini della redazione del tema natale, ho qualche dubbio sulla possibilità di utilizzo dove è richiesta una maggiore accuratezza, per esempio nelle tecniche direzionali e nei ritorni planetari. Iniziamo ad esplorare, nei prossimi post, qualche metodo alternativo.

Saturday, January 31, 2015

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.

Tuesday, January 27, 2015

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.

Tuesday, January 20, 2015

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.

Saturday, January 17, 2015

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.

Tuesday, March 25, 2014

Il modulo ctypes e la libreria libswe.so

Siccome la mia class SVG2Tk è ancora un po' indietro di cottura, apro una piccola parentesi sull'accesso diretto alle librerie condivise scritte in linguaggio C.

La libreria pyswisseph, che ho usato nei primi post, è un lavoro splendido e ci permette di effettuare chiamate alla libreria con istruzioni in Python nativo. Se pero' volessimo usare la libreria originaria, installata dal software center di Ubuntu, dovremmo confrontarci con un problema: non è scritta in Python ma in C e compilata verosimilmente con gcc, il compilatore GNU. La stessa pyswisseph è un insieme di funzioni wrapper (involucro) che fanno da ponte fra i due linguaggi. Anche noi, però, possiamo farlo direttamente dal nostro codice.

Tra le librerie integrate nell'installazione di Python c'è un modulo realizzato appositamente per chiamare una libreria in C, a patto che si conoscano le specifiche della libreria, i nomi delle funzioni e i parametri in input e in output.

Fortunatamente la libreria swisseph è molto ben documentata e con un po' di sforzo possiamo effettuare delle chiamate come se utilizzassimo una libreria pitonica. Vediamo come, prendendo a spunto un post del curatore delle pyswisseph Stanislas Marquis link.

Per prima cosa apriamo IDLE e scriviamo l'import di ctypes e di datetime. A questo punto creiamo una variabile shared_lib che costituirà il link alla libreria, attraverso il metodo CDLL di ctypes (o, in alternativa, il metodo cdll.LoadLibrary, sovrapponibile al precedente). L'installazione di swisseph da ubuntu software center mette la libreria libswe.so in una specifica posizione sul disco. Qualora volessimo spostarla nella stessa directory del sorgente Python sarà necessario cambiare il path della libreria stessa. Aggiungiamo inoltre due dizionari per i nomi dei pianeti e dei segni zodiacali.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# -*- coding: utf-8 -*-
import datetime
import ctypes as ct 

shared_lib = ct.CDLL('/usr/lib/x86_64-linux-gnu/libswe.so')

pianeti = {0:'Sole',
           1:'Luna',
           2:'Mercurio',
           3:'Venere',
           4:'Marte',
           5:'Giove',
           6:'Saturno',
           7:'Urano',
           8:'Nettuno',
           9:'Plutone',
           10:'Nodo lunare medio',
           11:'Nodo lunare vero'}

segni = {0:'Ariete',
         1:'Toro',
         2:'Gemelli',
         3:'Cancro',
         4:'Leone',
         5:'Vergine',
         6:'Bilancia',
         7:'Scorpione',
         8:'Sagittario',
         9:'Capricorno',
         10:'Acquario',
         11:'Pesci'}

Della libreria swisseph, in questo programma, uso solo due funzioni, una per il calcolo della data giuliana, l'altra per il calcolo della posizione dei pianeti. Vediamo la prima.


1
2
3
4
5
6
7
_julday=shared_lib.swe_julday
_julday.argtypes=[ct.c_int, ct.c_int, ct.c_int, ct.c_double, ct.c_int]
_julday.restype=ct.c_double

def julday(year, month, day, hour, gregflag=1):
    return _julday(ct.c_int(year), ct.c_int(month), ct.c_int(day),
                   ct.c_double(hour), ct.c_int(gregflag))

La prima parte è molto semplice. La libreria ctypes consente di accedere alla funzione nascosta nella shared library swe_julday definendo prima gli argtypes, cioè le tipologie di variabili da usare come parametri di input, poi il restype, cioè il tipo di variabile di ritorno dalla funzione (ricordiamo che il linguaggio C non consente la tipizzazione dinamica, cioè l'assegnazione arbitraria in runtime di un tipo all'atto della creazione della variabile. Essendo un linguaggio compilato, richiede, ovviamente, la tipizzazione di tipo statico).

La chiamata da codice Python risulta molto naturale, in questo caso i parametri passati alla funzione julday vengono trasferiti ai tipi definiti per l'interfaccia al linguaggio C. La funzione restituisce una variabile double, in Python corrisponde a una float e fornisce il valore della data giuliana calcolata da data e ora.

Un po' più complessa risulta, invece, la seconda funzione, che calcola le posizioni planetarie. Di seguito un estratto della guida Programming Interface to the Swiss Ephemeris:

2.1. The call parameters
swe_calc_ut() was introduced with Swisseph version 1.60 and makes planetary calculations a bit simpler. For the steps required, see the chapter  The programming steps to get a planet’s position.
swe_calc_ut() and swe_calc() work exactly the same way except that swe_calc() requires Ephemeris Time( more accurate: Dynamical Time ) as a parameter whereas swe_calc_ut() expects Universal Time. For common astrological calculations, you will only need swe_calc_ut() and will not have to think anymore about the conversion between Universal Time and Ephemeris Time.
swe_calc_ut() and swe_calc() compute positions of planets, asteroids, lunar nodes and apogees. They are defined as follows:
 
int swe_calc_ut ( double tjd_ut, int ipl, int iflag, double* xx, char* serr),
where
tjd_ut     =Julian day, Universal Time
ipl       =body number
iflag    =a 32 bit integer containing bit flags that indicate what kind of computation is wanted
xx       =array of 6 doubles for longitude, latitude, distance, speed in long., speed in lat., and speed in dist.
serr[256] =character string to return error messages in case of error.
 
and
int swe_calc(double tjd_et, int ipl, int iflag, double *xx, char *serr),
same but
tjd_et     =     Julian day, Ephemeris time,  where tjd_et = tjd_ut + swe_deltat(tjd_ut)

Nella dichiarazione della funzione swe_calc compaiono due puntatori (in Python non esistono, o meglio le variabili Python sono dei puntatori, ma non dilunghiamoci), cioè gli indirizzi di memoria nei quali possiamo trovare le variabili xx (che contiene un array di 6 doubles, in Python diremmo una tupla di 6 float) e serr (per l'eventuale messaggio di errore).

I tipi di variabili ctypes possono essere inizializzati nel codice Python con modalità differenti.

Per i tipi più semplici, ctypes.c_int(42) ctypes.c_float(3.1415), ad esempio, il valore tra parentesi viene conservato in una locazione di memoria da cui possono essere recuperati con l'espressione .value. Per esempio:

>>> a = ctypes.c_int(42)
>>> print a
c_int(42)
>>> print a.value
42
>>> b = ctypes.c_float(3.1415)
>>> print b
c_float(3.1414999961853027)
>>> print b.value
3.14149999619

Vediamo ora come si puo' effettuare la traduzione in Python della funzione swe_calc di swisseph utilizzando ctypes.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
_calc = shared_lib.swe_calc
_calc.argtypes = [ct.c_double, ct.c_int, ct.c_int, ct.POINTER(ct.c_double*6), ct.c_char_p]
_calc.restype = ct.c_int

def calc(julday, planet, flag = 0):

    xx = (ct.c_double * 6)()
    err = ct.c_char_p('')
    if _calc(ct.c_double(julday), ct.c_int(planet), ct.c_int(flag),
             ct.byref(xx), err) == 0:
        return [x for x in xx]
    return err.value

Come potete vedere, il tipo POINTER viene usato per tradurre il puntatore ad una variabile (preceduta da *) di C. Per la variabile serr si usa il tipo predefinito ctypes.c_char_p che è esso stesso un puntatore. La funzione accetta come parametri un double che riprende il risultato del calcolo della data giuliana, una variabile c_int per il numero di posizione del pianeta nella serie e infine il tipo c_int(flag), inizializzato a 0 perchè "If no bits are set, i.e. if iflag == 0, swe_calc() computes what common astrological ephemerides (as available in book shops) supply". Se il valore di ritorno di _calc è zero, cioè non viene resituito un errore, l'array di 6 double va a popolare una list di 6 float contenente vari elementi relativi al pianeta considerato. Di questi ci interessa solo l'elemento con indice 0, riferito alla longitudine in gradi. La lista viene generata con una list comprehension [x for x in xx].

L'ultima parte del programma contiene istruzioni per l'uso delle funzioni C che abbiamo creato in precedenza.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
def d2hm(x):
    lon30 = x % 30
    segno = int(x) / 30
    grado = int(lon30)
    minuto = int((lon30-grado)*60)
    return (grado, segno, minuto)


adesso = datetime.datetime.now()
in_questo_momento = julday(adesso.year,
       adesso.month,
       adesso.day,
       adesso.hour + adesso.minute/60.0 + adesso.second/3600.0)

for i in range(0,12):
    longitudine = calc(in_questo_momento,i)[0]
    pos = d2hm(longitudine)
    print "%18s %2d° %10s %2d'" % (pianeti[i], pos[0], segni[pos[1]], pos[2])

Lanciando il programma si ottiene il seguente risultato:

              Sole  4°     Ariete 49'
              Luna 24° Capricorno 33'
          Mercurio  9°      Pesci 42'
            Venere 18°   Acquario 18'
             Marte 23°   Bilancia 56'
             Giove 11°     Cancro  1'
           Saturno 22°  Scorpione 53'
             Urano 12°     Ariete  1'
           Nettuno  6°      Pesci  8'
           Plutone 13° Capricorno 28'
 Nodo lunare medio 29°   Bilancia 50'
  Nodo lunare vero 28°   Bilancia 36'

Riferito ad oggi 25 marzo alle 13.37 tempo di Greenwich

Spero che questo post sia stato utile per capire come interfacciarsi ad una shared library compilata in linguaggio C. Ulteriori informazioni sull'uso del modulo ctypes sono rintracciabili alla pagina http://docs.python.org/2/library/ctypes.html#

Se volete creare voi stessi la libreria libswe.so, scaricate il file a questo link. Decomprimetelo e lanciate l'istruzione make, che crea i file oggetto. Quindi create la libreria condivisa, nella stessa directory, con l'istruzione gcc -g -O9 -Wall -fPIC -lm -shared -o libswe.so swedate.o swehouse.o swejpl.o swemmoon.o swemplan.o swepcalc.o sweph.o swepdate.o swephlib.o swecl.o swehel.o. La libreria libswe.so comparirà nella directory. Copiatela dove avete collocato il sorgente .py di cui a questo post. Modificate il riferimento alla libreria nel vostro codice in shared_lib = ct.CDLL('libswe.so') e lanciate in esecuzione il codice python, funzionerà perfettamente.

Sunday, March 16, 2014

Quarto intermezzo - L'indice di concentrazione planetaria

Si deve ad André Barbault, insigne astrologo francese, l'introduzione di un indice sintetico che rappresenta il grado di maggiore o minore allontanamento reciproco dei pianeti, che, prescindendo dalle metodologie tradizionali dell'astrologia, basate sulla posizione dei pianeti nei segni e sugli aspetti, cioè distanze angolari discrete dotate di tolleranze contenute nell'ordine di qualche grado, consentirebbe la lettura degli eventi collettivi (astrologia mondiale) secondo cicli complessi che variano tra minimi e massimi nell'arco di giorni, mesi, anni o addirittura secoli. Per i necessari approfondimenti rimando all'opera dell'autore citato, in particolare a:

  1. A. Barbault. Il pronostico sperimentale in astrologia. Mursia ed. 1979
  2. A. Barbault. Astrologia mondiale (trad. Clara Negri). Armenia 1980

Per il calcolo dell'indice ho pensato di fornire un piccolo contributo personale ai lettori esperti di linguaggio Python, appoggiandomi ad alcune librerie esterne per Python:

La prima è pyephem di Brandon Rhodes, una libreria astronomica (non astrologica come swisseph di Astrodienst)

la seconda è pylab, dedicata alla grafica matematica, che fornisce una sintesi di numpy, scipy e matplotlib per fornire un ambiente simile all'ambiente MATLAB, realizzato per il linguaggio C.

Come sempre fornisco un listato (quick and dirty) per il calcolo dell'indice. Ho dato un orizzonte temporale di 50 anni, ma, ovviamente, chi ha già familiarità con Python puo' divertirsi a cambiare data di inizio, lunghezza del periodo, e altri parametri. Le librerie datetime e time sono incluse nell'installazione Python di base, quindi vanno solo importate. Per pylab nel repository di Ubuntu trovate il pacchetto python-matplotlib. Numpy può essere installata attraverso pip o dai repository di Ubuntu senza particolari difficoltà.

Ok, vi passo il codice per l'indice di concentrazione planetaria, considerando i pianeti lenti da Giove a Plutone. Nel grafico compare anche (quasi incomprensibile, lo so) la posizione all'inizio di ogni anno dei cinque pianeti considerati, nell'intervallo 0-360 gradi. Puo' aiutare a capire almeno quando si realizzano delle congiunzioni importanti.


import  pylab ,  numpy 
import  ephem 
import  time ,  datetime

planets = { 0 :ephem.Sun, 1 :ephem.Moon, 2 :ephem.Mercury, 3 :ephem.Venus,
            4 :ephem.Mars, 5 :ephem.Jupiter, 6 :ephem.Saturn, 7 :ephem.Uranus,
            8 :ephem.Neptune, 9 :ephem.Pluto}

def  longitude (x):
     return ephem.Ecliptic (x)

year =  1990 
month =  1 
day =  1 
hour =  0 
minute =  0 
time_span =  50

initial_date = datetime.datetime (year, month, day, hour, minute)

x = []
y = []
jup = []
sat = []
ura = []
nep = []
plu = []

for i in  range ( 0 , time_span):
     date = datetime.datetime (year + i, month, day, hour, minute)
     b = []
     icp =  0 
     for key in  range ( 5 , 10 ):
          a = planets[key](date, epoch = date)
          b.append(ephem.degrees(ephem.Ecliptic(a).lon) / ephem.pi * 180 )
     for j in  range ( 0 , len (b)):
          for k in  range (j, len (b)):
               _icp = (b [k] - b [j])
               if _icp <  0 :
                    _icp +=  360 
               if _icp >  180 :
                    _icp =  360 - _icp
               icp += _icp
     x.append (i + year)
     y.append (icp)
     jup.append (b[0])
     sat.append (b[1])
     ura.append (b[2])
     nep.append (b[3])
     plu.append (b[4])

pylab.figure(1,figsize =(15 , 10 ))
pylab.title ( 'Index planetary concentration' )
pylab.xlabel ( 'year' )
pylab.ylabel ( 'ICP' )
pylab.xticks (numpy.linspace (year, year + time_span, 50 , endpoint = True ), rotation = - 90 )
pylab.plot (x, y, 'r' , label = 'ICP' )
pylab.plot (x, jup, 'bo' , label = 'jupiter' )
pylab.plot (x, sat, 'g ^' , label = 'saturn' )
pylab.plot (x, ura, 'c.' , label = 'uranus' )
pylab.plot (x, nep, 'm,' , label = 'neptune' )
pylab.plot (x, plu, 'k.' , label = 'bar' )

pylab.grid (axis = 'Both' )
pylab.legend (loc = 'best' )
pylab.show ()

A fra qualche giorno con l'ultima puntata di SVG2Tk.

Friday, March 14, 2014

XML - ricerca ricorsiva rivisitata

Le cose non vanno sempre lisce quando si esplora un protocollo. Nel caso dei file SVG il procedimento seguito finora puo' non funzionare (anzi non funziona quasi mai) perchè il protocollo SVG è lontano dall'avere sempre gli stessi contenuti e la stessa disposizione degli stessi. L'unica cosa certa è che il file SVG è un file XML, quindi composto da tag nidificate, che i dati sono prevalentemente negli attributi piuttosto che in testo libero e che dobbiamo quindi trovare una modalità di esplorazione, sempre ricorsiva, che in primo luogo ci dia la certezza di estrarre tutti i tag e relativi attributi e quindi ci consenta la comparazione dei tag con modelli tipo e relativi pattern di estrazione con le regexp.

Ho quindi riscritto (e allego di seguito) il prototipo della funzione xmlrecur per adattarsi a qualsiasi file xml o SVG di cui vogliamo fare successive conversioni.


#|/usr/bin/env python
from xml.etree import ElementTree as ET

def xmlrecur(x):
    print '<elemento',x,'>'
    print '<tag>', x.tag
    print '<attributo>',x.attrib
    for key in x.attrib:
        print 'singolo attributo', key, '-', x.attrib[key]
    print '<testo>',x.text
    print '-'*80
        
    for i in x:
            xmlrecur(i)

tree = ET.parse('xxx.xml')
root = tree.getroot()
xmlrecur(root)

Vi prego di notare che la proprietà attrib è un dizionario, di cui si puo' fare la scansione. Su questa particolarità concentreremo l'analisi per l'identificazione delle word caratteristiche di nostro interesse.

Ancora qualche giorno per finire il tutto. A presto.

Friday, March 7, 2014

Ulteriore miglioramento del codice

Avevo annunciato l'intenzione di eliminare, dalla classe SVG2Tk, l'uso di file di appoggio esterni.

Nel file che allego in questo post mostro come si puo' fare. Un difetto fondamentale della classe, come è stata progettata finora, è la pesantezza dell'analisi del file SVG, la necessità di riversare le informazioni ottenuta in file esterni, che vengono scritti e riletti senza che ce ne sia un'effettiva necessità. Un ulteriore miglioramento puo' essere realizzato con l'uso dei generatori, argomento tutt'altro che facile da capire, ma che cercherò di rendere nel modo più semplice possibile.

In primo luogo, disponendo di una metodo ricorsivo per l'identificazione degli attributi fondamentali del file SVG, utili al nostro lavoro di traduzione, possiamo cercare di fare in un'unica passata la lettura e la trascrizione dei punti di tracciato, senza salvataggi intermedi e senza creare lunghe liste in memoria, che, nonostante le notevoli migliorie delle tecnologie attualmente disponibili, rischierebbe di saturare la memoria RAM e allungare i tempi.

I generatori sono funzioni speciali, in cui i risultati ottenuti dallo scorrimento di una sequenza vengono restituiti on demand, cioè man mano che il programma ne richiede la disponibilità.

Ciò significa che il file viene letto in modo da assicurare la immediata disponibilità di un'informazione, che puo' essere trattata ignorando le informazioni successive.

Non essendo più necessari i file intermedi di appoggio, forniamo solo definizione del file SVG, e la forniamo alla classe attraverso il suo inizializzatore o costruttore.

Per agevolare la comprensione, spezzo il codice nelle sezioni corrispondenti ai metodi, che pubblico una alla volta per commentarle immediatamente dopo. Il lavoro non è ancora concluso, usero' ancora una volta (l'ultima) un file per riversare i dati sul disco fisso, ma sarà evidente il guadagno di tempo e memoria e la maggiore snellezza del codice.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
# -*- coding: utf-8 -*-
#|/usr/bin/env python
import re
import xml.etree.ElementTree as ET

class SVG2Tk:
    """La classe è finalizzata ad assicurare la traduzione, per quanto possibile
    fedele alle specifiche SVG, negli oggetti caratteristici del widget Canvas
    di Tkinter"""

    def __init__(self, file_in):
        
        self.file_in = open(file_in, 'r')

La prima parte non è di difficile comprensione, riprende più o meno i precedenti contenuti. Ho eliminato i parametri riferiti ai nomi dei file di appoggio e ho inserito delle righe di commento che vanno a comporre la documentazione della classe così come si puo' ottenere con i normali strumenti di documentazione di Python: lo vedremo più chiaramente alla fine.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
    def s_match(self, pattern, d):
        """restituisce una stringa contenente la substringa che soddisfa
            il pattern fornito alla funzione, e la stringa risultante dopo
            l'estrazione della substringa"""
        
        matchObj = re.match(pattern, d)
        if matchObj:
            string = matchObj.group()
            d = d.replace(string, "")
            s_string = re.split(r'[;\s]*',string)
            string = ','.join(s_string[:-1])
            return (string, d)
        return None

Anche questa sezione, che utilizza le regular expressions per il match tra stringhe di SVG e modelli, non differisce dalla precedente.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
    def xmlrecur(self,x):
        "restituisce l'insieme degli attributi che soddisfano le stringhe di ricerca "
        for i in x:
            if set(['id','style','d']).issubset(i.attrib):
                yield (i.attrib['id'], i.attrib['style'], i.attrib['d'])
            elif set(['width','height']).issubset(i.attrib):
                yield ('width - height', i.attrib['width'], i.attrib['height'])
            else:
                for j in self.xmlrecur(i):
                    yield j

La routine xmlrecur è stata, invece, pesantemente rimaneggiata. Ogni volta che l'identificazione dell'attributo è positiva, uso la word yield che è simile al return, cioè alla restituzione dei risultati delle attività previste dal metodo, ma avviene ogni volta che ho necessità di fornire al codice chiamante dei valori senza interrompere la scansione dell'iterable, in questo caso del file di origine dei dati, che viene solo temporaneamente sospesa. Ad ogni successiva chiamata verrà ripresa la scansione e verranno forniti ulteriori risultati.

Ciò che rende di difficile comprensione il metodo è il fatto che apparentemente viene restituito al codice chiamante un dato, mentre in realtà viene restituito un metodo, che si chiama generatore.

Vediamo la sezione finale e ne discutiamo immediatamente.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
    def re_load(self):

        tree = ET.parse(self.file_in)
        root = tree.getroot()
        generatore = self.xmlrecur(root)
        tk_out = open('tk_out','w')
        for i in generatore:
            if i[0] == 'width - height':
                self.width  = float(i[1])
                self.height = float(i[2])
            else:
                self.id    = i[0]
                self.style = i[1]
                self.d     = i[2]
                try:
                    pattern=[
                    r'([M]{1}[\s]*[-\d]*[\.]?[\d]*[\,\s]?[-\d]*[\.]?[\d]*[\s]*)',
                    r'([L]{1}[\s]*[-\d]*[\.]?[\d]*[\,\s]{1}[-\d]*[\.]?[\d]*[\s]*)',
                    r'([C]{1}[\s]*)(([-\d]*[\.]?[\d]*[\,\s]{1}[-\d]*[\.]?[\d]*[\s]*){3})',
                    r'([z]{1}[\s]*)'
                    ]
                    
                    while len(self.d) > 0:
                        for p in pattern:
                            ret = self.s_match(p, self.d)
                            if ret != None:
                                string = ret[0]
                                self.d = ret[1]
                                tk_out.write(string+'\n')                            
                except Exception:
                    print Exception
                    tk_out.close()
                    break

if __name__ == '__main__':
    app = SVG2Tk('EU-Italy.svg')
    app.re_load()

La funzione re_load è simile alla precedente, ma ne differisce per il fatto che i dati ottenuti dalla scansione vengono caricati, non come dati, ma come metodo, sull'item generatore da cui vengono estratti uno alla volta e sottoposti al matching. Non è intuitivo, il generatore non contiene alcun dato, non è una lista, bensì una chiamata iterativa al metodo recurxml che estrae i dati on demand. Ci ho messo parecchio a capire come funziona questo metodo e rinuncio a tentare di spiegarlo ulteriormente. Se avete voglia di discuterlo mandatemi un commento e ne parliamo insieme.

Nella sezione finale inizializzo la class col solo nome del file SVG e poi faccio una chiamata al metodo re_load che fa girare il tutto. Nel prossimo post eliminerò anche l'ultimo puntello, cioè il file di appoggio, e farò in modo che l'estrazione avvenga in streaming con la creazione degli oggetti di Tkinter.Canvas. Alla prossima

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