VdS-Journal 88: Berechnung von Lichterscheinungen auf dem Mond

Die Wiederkehr von Lichterscheinungen erfordert, dass die Licht/Schatten-Grenze auf demselben Längenkreis verläuft, wie zur Referenzbeobachtung. Die Grenze des Morgenterminators wird als Colongitude bezeichnet. Das nachfolgende Programm berechnet diesen Wert und sucht nach Zeitpunkten, wo dieselbe Konstellation auftritt.
# Hesiodus.py - Berechnung von Lichterscheinungen au dem Mond
# Wie in Meeus: AA werde alle Winkel in Grad geführt und erst in den 
# Winkelfunktionen #umgewandelt. Dies erleichtert die Tests.

from math import *
import math
toRad = math.pi / 180

# Tab.  47.A, 339
LDterms = [
[0,  0,  1,  0,    6288774,    -20905355],
[2,  0, -1,  0,    1274027,     -3699111],
[2,  0,  0,  0,     658314,     -2955968],
[0,  0,  2,  0,     213618,      -569925],
[0,  1,  0,  0,    -185116,        48888],
[0,  0,  0,  2,    -114332,        -3149],
[2,  0, -2,  0,      58793,       246158],
[2, -1, -1,  0,      57066,      -152138],
[2,  0,  1,  0,      53322,      -170733],
[2, -1,  0,  0,      45758,      -204586],
[0,  1, -1,  0,     -40923,      -129620],
[1,  0,  0,  0,     -34720,       108743],
[0,  1,  1,  0,     -30383,       104755],
[2,  0,  0, -2,      15327,        10321],
[0,  0,  1,  2,     -12528,            0],
[0,  0,  1, -2,      10980,        79661],
[4,  0, -1,  0,      10675,       -34782],
[0,  0,  3,  0,      10034,       -23210],
[4,  0, -2,  0,       8548,       -21636],
[2,  1, -1,  0,      -7888,        24208],
[2,  1,  0,  0,      -6766,        30824],
[1,  0, -1,  0,      -5163,        -8379],
[1,  1,  0,  0,       4987,       -16675],
[2, -1,  1,  0,       4036,       -12831],
[2,  0,  2,  0,       3994,       -10445],
[4,  0,  0,  0,       3861,       -11650],
[2,  0, -3,  0,       3665,        14403],
[0,  1, -2,  0,      -2689,        -7003],
[2,  0, -1,  2,      -2602,            0],
[2, -1, -2,  0,       2390,        10056],
[1,  0,  1,  0,      -2348,         6322],
[2, -2,  0,  0,       2236,        -9884],
[0,  1,  2,  0,      -2120,         5751],
[0,  2,  0,  0,      -2069,            0],
[2, -2, -1,  0,       2048,        -4950],
[2,  0,  1, -2,      -1773,         4130],
[2,  0,  0,  2,      -1595,            0],
[4, -1, -1,  0,       1215,        -3958],
[0,  0,  2,  2,      -1110,            0],
[3,  0, -1,  0,       -892,         3258],
[2,  1,  1,  0,       -810,         2616],
[4, -1, -2,  0,        759,        -1897],
[0,  2, -1,  0,       -713,        -2117],
[2,  2, -1,  0,       -700,         2354],
[2,  1, -2,  0,        691,            0],
[2, -1,  0, -2,        596,            0],
[4,  0,  1,  0,        549,        -1423],
[0,  0,  4,  0,        537,        -1117],
[4, -1,  0,  0,        520,        -1571],
[1,  0, -2,  0,       -487,        -1739],
[2,  1,  0, -2,       -399,            0],
[0,  0,  2, -2,       -381,        -4421],
[1,  1,  1,  0,        351,            0],
[3,  0, -2,  0,       -340,            0],
[4,  0, -3,  0,        330,            0],
[2, -1,  2,  0,        327,            0],
[0,  2,  1,  0,       -323,         1165],
[1,  1, -1,  0,        299,            0],
[2,  0,  3,  0,        294,            0],
[2,  0, -1, -2,          0,         8752]]


# Tab.  47b, 341
Lterms = [
[0,  0,  0,  1,      5128122],
[0,  0,  1,  1,       280602],
[0,  0,  1, -1,       277693],
[2,  0,  0, -1,       173237],
[2,  0, -1,  1,        55413],
[2,  0, -1, -1,        46271],
[2,  0,  0,  1,        32573],
[0,  0,  2,  1,        17198],
[2,  0,  1, -1,         9266],
[0,  0,  2, -1,         8822],
[2, -1,  0, -1,         8216],
[2,  0, -2, -1,         4324],
[2,  0,  1,  1,         4200],
[2,  1,  0, -1,        -3359],
[2, -1, -1,  1,         2463],
[2, -1,  0,  1,         2211],
[2, -1, -1, -1,         2065],
[0,  1, -1, -1,        -1870],
[4,  0, -1, -1,         1828],
[0,  1,  0,  1,        -1794],
[0,  0,  0,  3,        -1749],
[0,  1, -1,  1,        -1565],
[1,  0,  0,  1,        -1491],
[0,  1,  1,  1,        -1475],
[0,  1,  1, -1,        -1410],
[0,  1,  0, -1,        -1344],
[1,  0,  0, -1,        -1335],
[0,  0,  3,  1,         1107],
[4,  0,  0, -1,         1021],
[4,  0, -1,  1,          833],
[0,  0,  1, -3,          777],
[4,  0, -2,  1,          671],
[2,  0,  0, -3,          607],
[2,  0,  2, -1,          596],
[2, -1,  1, -1,          491],
[2,  0, -2,  1,         -451],
[0,  0,  3, -1,          439],
[2,  0,  2,  1,          422],
[2,  0, -3, -1,          421],
[2,  1, -1,  1,         -366],
[2,  1,  0,  1,         -351],
[4,  0,  0,  1,          331],
[2, -1,  1,  1,          315],
[2, -2,  0, -1,          302],
[0,  0,  1,  3,         -283],
[2,  1,  1, -1,         -229],
[1,  1,  0, -1,          223],
[1,  1,  0,  1,          223],
[0,  1, -2, -1,         -220],
[2,  1, -1, -1,         -220],
[1,  0,  1,  1,         -185],
[2, -1, -2, -1,          181],
[0,  1,  2,  1,         -177],
[4,  0, -2, -1,          176],
[4, -1, -1, -1,          166],
[1,  0,  1, -1,         -164],
[4,  0,  1, -1,          132],
[1,  0, -1, -1,         -119],
[4, -1,  0, -1,          115],
[2, -2,  0,  1,          107]]

def getJD(y, m, d) :
    if (m < 3):
        y = y - 1
        m = m + 12 
    a = (y // 100)
    b = 2 - a + (a // 4.0)
    jd = floor(365.25 * (y + 4716)) 
    jd = jd + int(30.6001 * (m + 1)) + d + b - 1524.5
    return jd

def getKalender(myjd):  # vgl.  Wikipedia
    tag = 0 # eingegliedert, für statische Methode
    tag = 0
    monat = 0
    jahr = 0
    stunde = 0
    minute = 0

    z = int((myjd + 0.5)) # Z = Int(JD + 0,5)
    f = (myjd + 0.5) - z # F = Frac(JD + 0,5)
    if (z < 2299161):
       a = z # wenn Z < 2299161 dann A = Z Ergebnis julianisch
    else:  # gregorianisch:
    	g = (int)((z - 1867216.25) / 36524.25) # g = Int((Z - 1867216,25) / 36524,25)
    	a = z + 1 + g - (int)(g / 4) # A = Z + 1 + g - Int(g/4)
    b = a + 1524 # B = A + 1524
    c = int((b - 122.1) / 365.25) # C = Int((B-122,1) / 365,25)
    d = int(365.25 * c) # D = Int(365,25 * C)
    e = int((b - d) / 30.6001) # E = Int((B-D) / 30,6001)
    tagDouble = b - d - int(30.6001 * e) + f # Tag = B - D - Int(30,6001*E) + F # Tag, inklusive Tagesbruchteil
    tag = int(tagDouble)
    if (e < 14) :
        monat = e - 1
    else :
        monat = e - 13 # wenn E<14 dann Monat = E - 1 sonst Monat = E - 13
    if (monat > 2):
       jahr = c - 4716
    else :
        jahr = c - 4715 # wenn Monat>2 dann Jahr = C - 4716 sonst Jahr = C - 4715
    bruchteil = tagDouble - tag
    stunde = int(24 * bruchteil)
    minute = int(60 * (24 * bruchteil - stunde))
    return jahr, monat, tag, stunde, minute



def mod360(x) :
    x = 360 * (x / 360 - int(x / 360))
    if (x < 0):
       x = x + 360
    return x

def moon(JD) :
    T = (JD - 2451545) / 36525
    T2 = T * T 
    T3 = T2 * T 
    T4 = T3 * T

    # mittl.  Länge samt Lichtzeit
    Ls = 218.3164477 + 481267.88123421 * T - 0.0015786 * T2 + T3 / 538841 - T4 / 65194000
    Ls = mod360(Ls)

    # mittlere Elongation
    D = 297.8501921 + 445267.1114034 * T - 0.0018819 * T2 + T3 / 545868 - T4 / 113065000
    D = mod360(D)

    # mittl.  Anomalie der Sonne
    M = 357.5291092 + 35999.0502909 * T - 0.0001536 * T2 + T3 / 24490000
    M = mod360(M)

    # mittl.  ANomalie des Mondes
    Ms = 134.9633964 + 477198.8675055 * T + 0.0087414 * T2 + T3 / 69699 - T4 / 14712000
    Ms = mod360(Ms)

    # Argument der Breite
    F = 93.2720950 + 483202.0175233 * T - 0.0036539 * T2 - T3 / 3526000 + T4 / 863310000
    F = mod360(F)

    A1 = 119.75 + 131.849 * T
    A1 = mod360(A1)
    A2 = 53.09 + 479264.290 * T
    A2 = mod360(A2)
    A3 = 313.45 + 481266.484 * T
    A3 = mod360(A3)

    E = 1 - 0.002516 * T - 0.0000074 * T2
    E2 = E * E

    # Summe der periodischen Terme für Länge und Abstand
    SL = 0
    Sr = 0
    for i in range(60):
        a = LDterms[i][0] * D + LDterms[i][1] * M + LDterms[i][2] * Ms + LDterms[i][3] * F
        a = a * math.pi / 180 # in Radiant
        if ((LDterms[i][1] == 2) or (LDterms[i][1] == -2)) :
            SL += E2 * LDterms[i][4] * sin(a)
            Sr += E2 * LDterms[i][5] * cos(a)
        elif ((LDterms[i][1] == 1) or (LDterms[i][1] == -1)) :
            SL += E * LDterms[i][4] * sin(a)
            Sr += E * LDterms[i][5] * cos(a)
        else :
            SL += LDterms[i][4] * sin(a)
            Sr += LDterms[i][5] * cos(a)

    # Summe der periodischen Terme für die Breite
    Sb = 0
    for i in range(60):
        a = Lterms[i][0] * D + Lterms[i][1] * M + Lterms[i][2] * Ms + Lterms[i][3] * F
        a = a * math.pi / 180 # in Radiant
        if ((Lterms[i][1] == 2) or (Lterms[i][1] == -2)) :
            Sb += E2 * Lterms[i][4] * sin(a)
        elif ((Lterms[i][1] == 1) or (Lterms[i][1] == -1)): 
            Sb += E * Lterms[i][4] * sin(a)
        else:
            Sb += Lterms[i][4] * sin(a)
    

    # Terme für Jupiter, Venus und Erdabplattung
    SL += 3958 * sin(toRad * A1) + 1962 * sin(toRad * (Ls - F)) + 318 * sin(toRad * A2)
    Sb += -2235 * sin(toRad * Ls) + 382 * sin(toRad * A3) + 175 * sin(toRad * (A1 - F)) + 175 * sin(toRad * (A1 + F)) + 127 * sin(toRad * (Ls - Ms)) - 115 * sin(toRad * (Ls + Ms))

    # die Koordinaten selbst
    delta = (385000.56 + Sr / 1000.0)  # / AU
    lambdaM = mod360(Ls + SL / 1000000.0)
    beta = Sb / 1000000.0
    pi = asin(6378.14 / delta) / toRad

    return lambdaM, beta, delta, pi, Ls, D, M, Ms, F, E

# Sonne: Kap.  25, 163
def sonne(JD) :

    T = (JD - 2451545) / 36525
    T2 = T * T 
    T3 = T2 * T 
    T4 = T3 * T 
    T5 = T4 * T
    Lo = 280.4664567 + 36000.76982779 * T + 0.003032028 * T2 + T3 / 49931 - T4 / 15299 - T5 / 1988000

    Lo = mod360(Lo)
    M = 357.52911 + 35999.050 * T - 0.0001537 * T2
    M = mod360(M)
    e = 0.016708617 - 0.000042037 * T - 0.0000001237 * T2
    C = (1.9146 - 0.004817 * T - 0.000014 * T2) * sin(toRad * M)
    C = C + (0.019993 - 0.000101 * T) * sin(2 * toRad * M)
    C = C + 0.00029 * sin(3 * toRad * M)
    C = mod360(C)

    trueLon = Lo + C # wahre Länge
    trueLon = mod360(trueLon)

    trueAnom = M + C
    R = (1.000001018 * (1.0 - e * e)) / (1.0 + e * cos(toRad * trueAnom))

    Omega = 125.0445479 - 1934.1362891 * T + 0.0020754 * T2 + T3 / 467441.0 - T4 / 60616000

    # scheinbare Länge
    lambdaS = trueLon - 0.00569 - 0.00478 * sin(toRad * Omega)
    lambdaS = mod360(lambdaS)
    
    return lambdaS, R, Omega, Lo

# Nutation und Schiefe der Ekliptik: Kap.  22 143
def nutation(JD, Omega, Lo, Ls) :
    T = (JD - 2451545) / 36525
    T2 = T * T 
    T3 = T2 * T
    deltaPsi = -17.20 * sin(toRad * Omega) - 1.32 * sin(2 * toRad * Lo) - 0.23 * sin(2 * toRad * Ls) + 0.21 * sin(2 * toRad * Omega)
    deltaPsi /= 3600

    deltaEps = 9.2 * cos(toRad * Omega) + 0.57 * cos(2 * toRad * Lo) + 0.10 * cos(2 * toRad * Ls) - 0.09 * cos(2 * toRad * Omega)
    deltaEps /= 3600

    epsilon = 23.43929111 + (-46.8150 * T - 0.00059 * T2 + 0.001813 * T3) / 3600 #22.2
    epsilon = 23.4392911 * 3600 - 46.8150 * T - 0.00059 * T2 + 0.001813 * T3
    epsilon = epsilon / 3600
    epsilon = epsilon + deltaEps
    return deltaPsi, deltaEps, epsilon

def colongitude(JD) :
 
    lambdaM, beta, delta, pi, Ls, D, M, Ms, F, E = moon(JD)
    lambdaS, R, Omega, Lo = sonne(JD)
    deltaPsi, deltaEps, epsilon = nutation(JD, Omega, Lo, Ls)

    T = (JD - 2451545) / 36525
    I = 1.54242 # Neigung des Mondäquators
    
    # die Libration 53.1
    W = lambdaM - Omega
    W = mod360(W)
    A = atan2((sin(toRad * W) * cos(toRad * beta) * cos(toRad * I) - sin(toRad * beta) * sin(toRad * I)), (cos(toRad * W) * cos(toRad * beta)))
    A = A / toRad
    A = mod360(A)
    ls = A - F
    sbs = -sin(toRad * W) * cos(toRad * beta) * sin(toRad * I) - sin(toRad * beta) * cos(toRad * I)
    bs = asin(sbs) / toRad


    K1 = 119.75 + 131.849 * T
    K2 = 72.56 + 20.186 * T

    rho = (-0.02752 * cos(toRad * Ms) - 0.02245 * sin(toRad * F) + 0.00684 * cos(toRad * Ms - 2 * toRad * F) + -0.00293 * cos(2 * toRad * F) - 0.00085 * cos(2 * toRad * F - 2 * toRad * D) - 0.00054 * cos(toRad * Ms - 2 * toRad * D) + -0.00020 * sin(toRad * Ms + toRad * F) - 0.00020 * cos(toRad * Ms + 2 * toRad * F) + -0.00020 * cos(toRad * Ms - toRad * F) + 0.00014 * cos(toRad * Ms + 2 * toRad * F - 2 * toRad * D))

    sigma = (-0.02816 * sin(toRad * Ms) + 0.02244 * cos(toRad * F) - 0.00682 * sin(toRad * Ms - 2 * toRad * F) + -0.00279 * sin(2 * toRad * F) - 0.00083 * sin(2 * toRad * F - 2 * toRad * D) + 0.00069 * sin(toRad * Ms - 2 * toRad * D) + 0.00040 * cos(toRad * Ms + toRad * F) - 0.00025 * sin(2 * toRad * Ms) + -0.00023 * sin(toRad * Ms + 2 * toRad * F) + 0.00020 * cos(toRad * Ms - toRad * F) + 0.00019 * sin(toRad * Ms - toRad * F) + 0.00013 * sin(toRad * Ms + 2 * toRad * F - 2 * toRad * D) + -0.00010 * cos(toRad * Ms - 3 * toRad * F))

    tau = (0.02520 * E * sin(toRad * M) + 0.00473 * sin(2 * toRad * Ms - 2 * toRad * F) - 0.00467 * sin(toRad * Ms) + 0.00396 * sin(toRad * K1) + 0.00276 * sin(2 * toRad * Ms - 2 * toRad * D) + 0.00196 * sin(toRad * Omega) + -0.00183 * cos(toRad * Ms - toRad * F) + 0.00115 * sin(toRad * Ms - 2 * toRad * D) + -0.00096 * sin(toRad * Ms - toRad * D) + 0.00046 * sin(2 * toRad * F - 2 * toRad * D) + -0.00039 * sin(toRad * Ms - toRad * F) - 0.00032 * sin(toRad * Ms - toRad * M - toRad * D) + 0.00027 * sin(2 * toRad * Ms - toRad * M - 2 * toRad * D) + 0.00023 * sin(toRad * K2) + -0.00014 * sin(2 * toRad * D) + 0.00014 * cos(2 * toRad * Ms - 2 * toRad * F) + -0.00012 * sin(toRad * Ms - 2 * toRad * F) - 0.00012 * sin(2 * toRad * Ms) + 0.00011 * sin(2 * toRad * Ms - 2 * toRad * M - 2 * toRad * D))

    lss = -tau + (rho * cos(toRad * A) + sigma * sin(toRad * A)) * tan(toRad * bs)
    bss = sigma * cos(toRad * A) - rho * sin(toRad * A)
    l = ls + lss
    b = bs + bss

    V = Omega + deltaPsi + sigma / sin(toRad * I)
    X = sin(toRad * I + toRad * rho) * sin(toRad * V)
    Y = sin(toRad * I + toRad * rho) * cos(toRad * V) * cos(toRad * epsilon) - cos(toRad * I + toRad * rho) * sin(toRad * epsilon)
    omega = atan2(X, Y) / toRad
    alpha = atan2(sin(toRad * lambdaM) * cos(toRad * epsilon) - tan(toRad * beta) * sin(toRad * epsilon), cos(toRad * lambdaM)) / toRad

    sP = sqrt(X * X + Y * Y) * cos(toRad * alpha - toRad * omega)
    sP = sP / cos(toRad * b)
    P = asin(sP) / toRad

    # die selenographische Sonnenposition
    AE = 149597870.7 # km

    lambdaH = lambdaS + 180 + delta / (AE * R) * 180 / math.pi * cos(toRad * beta) * sin(toRad * (lambdaS - lambdaM))
    betaH = delta / (AE * R) * beta
    # 53.1 neu
    Wo = lambdaH - deltaPsi - Omega
    Wo = mod360(Wo)
    Ao = atan2((sin(toRad * Wo) * cos(toRad * betaH) * cos(toRad * I) - sin(toRad * betaH) * sin(toRad * I)), (cos(toRad * Wo) * cos(toRad * betaH)))
    Ao = Ao / toRad
    Ao = mod360(Ao)

    lso = Ao - F
    sbso = -sin(toRad * Wo) * cos(toRad * betaH) * sin(toRad * I) - sin(toRad * betaH) * cos(toRad * I)
    bso = asin(sbso) / toRad

    lsso = -tau + (rho * cos(toRad * Ao) + sigma * sin(toRad * Ao)) * tan(toRad * bso)
    bsso = sigma * cos(toRad * Ao) - rho * sin(toRad * Ao)

    lo = lso + lsso
    bo = bso + bsso
    co = 90 - lo
    co = mod360(co)
    return co, b


# Hauptprogramm
j = 2020 # Startjahr
J = 1 # Anzahl Jahre
coZiel = 18.40 # was wird gesucht?  Hesiodus: coZiel=18,40°

vorigeCo = 7e77
JD = getJD(j, 1, 1.0)
deltaJD = 1 / 24.0 # 2,4h

for i in range(24 * 366 * J):
    co,b = colongitude(JD)
    if ((co > coZiel) and (vorigeCo < coZiel)):
        JDHesiodus=JD-1/24*(co-coZiel)/(co-vorigeCo);
        jahr, monat, tag, stunde, minute=getKalender(JDHesiodus)
        print(tag,monat,jahr," \t",stunde," ",minute," UT.  b=",b)
    JD = JD+deltaJD
    vorigeCo = co
Uwe Pilz, März 2020