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.