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.