VdS-Journal 95: Zentrallinie einer Sonnenfinsternis

Die Berechnung erfolgt minütlich. Vor Benutzung müssen die sog. Bessel-Elemente der gewünschten Finsternis eingesetzt werden.
Geographische Koordinaten werden astronomisch richtig ausgegebe, also Ost = negativ. Das ist umgekehrt zur Google-Maps-Anzeige.
# sofiBessel2 : Zentrallinie einer Sonnenfinsternis
# Meeus, Elements of solar eclipses §13

# Besselsche Elemente
# yr=2026; mnt=8; day=12
T0 = 18
x =  [ 0.4755140, 0.5189249, -0.0000773, -0.0000080]
y =  [ 0.7711830, -0.2301680, -0.0001246, 0.0000038]
d =  [ 14.7966700, -0.0120650, -0.0000030, 0.0]
l1 = [ 0.5379550, 0.0000939, -0.0000121, 0.0]
l2 = [ -0.0081420, 0.0000935, -0.0000121, 0.0]
u =  [ 88.7477900, 15.0030900, 0.0000000, 0.0]
tanf1 = 0.0046141
tanf2 = 0.0045911

from math import *


# Winkelfunktionen in Gradmaß, wie im Meeus
def Sin(w):    return sin(pi * w / 180)
def Cos(w):    return cos(pi * w / 180)
def Tan(w):    return tan(pi * w / 180)
def Atan2(y, x):    return 180 * atan2(y, x) / pi
def Asin(x):    return 180 * asin(x) / pi
def Acos(x):    return 180 * acos(x) / pi
def Atan(x):    return 180 * atan(x) / pi
def sq(x):    return x*x


def solveQuadrant(sinX,cosX):
    if (sinX>=0 and cosX>=0): return Asin(sinX);
    if (sinX<0 and cosX>=0): return Asin(sinX);
    if (sinX<0 and cosX<0): return -Acos(cosX);
    if (sinX>=0 and cosX<0): return Acos(cosX);

def calcEclipse(t):
    # Koordinaten zum Zeitpunkt
    X = 0
    Y = 0
    D = 0
    L1 = 0
    L2 = 0
    M = 0
    for j in range(3, -1, -1):
        X = X * t + x[j]        # die Koordinaten
        Y = Y * t + y[j] 
        D = D * t + d[j]        # Deklination ...
        M = M * t + u[j]        # ... und Stundenwinkel im orstfesten Äquatorsystem
        L1 = L1 * t + l1[j]     # die Radien der Schattenkegel, Halbschatten ...
        L2 = L2 * t + l2[j]     # ... und Kernschatten
    # die Ableitungen / Veränderungen des Ortes i.f. Fundamental-E.
    XS = x[1] + 2 * x[2] * t + 3 * x[3] * t * t
    YS = y[1] + 2 * y[2] * t + 3 * y[3] * t * t
    cd = Cos(D)
    w = 1.0 / (sqrt(1 - 0.006694385 * cd * cd))
    p = u[1] / 57.2957795
    b = YS - p * X * Sin(D)
    c = XS + p * Y * Sin(D)
    y1 = w * Y
    b1 = w * Sin(D)
    b2 = 0.99664719 * w * Cos(D)
    if 1 - X * X - y1 * y1 < 0:
        return
    B = sqrt(1 - X * X - y1 * y1)
    phi1 = Asin(B * b1 + y1 * b2)
    # von Greg Miller
    sinH=X / Cos(phi1)
    cosH=(B * b2 - y1 * b1) / Cos(phi1)
    H=solveQuadrant(sinH, cosH)
   
    
    # Berechnung der topozentrischen Koordinaten der Schattenachse = Zentral-Linie
    phi = Atan(1.00336409 * Tan(phi1))
    Lambda = M - H - 0.00417807 * dt  # dt in Sekunden
    #  Berechnung des Totalitätsstreifens und der Dauer der Finsternis im Zentrum
    L2S = L2 - B * tanf2 # l2S negativ: Ringförmig
    typ="r"
    if L2S<0:
        typ="t"
    a = c - p * B * Cos(D)
    n2 = a * a + b * b
    if n2 < 0:
        n = -777  # kein n
        dauer = -777  # keine Dauer
    else:
        n=sqrt(n2)
        dauer = abs(7200 * L2S / n)
    sinh = Sin(D) * Sin(phi) + Cos(D) * Cos(phi) * Cos(H)
    h = Asin(sinh)
    K2 = B * B + (X * a + Y * b) * (X * a + Y * b) / (n * n)
    K = sqrt(K2)
    pfadbreite = 12756 * abs(L2S) / K
    # Berechnung der Breite der scheinbaren Größe des Mondes (in Sonnendurchmessern)
    L1S=l1[0]+l1[1]*t+l2[2]*t*t-B*tanf1
    A = (L1S - L2S) / (L1S + L2S)
    if Lambda>360:
        Lambda=Lambda-360
    if (Lambda>180):
        Lambda=Lambda-360
    print("%3.0f  %8.4f  %8.4f  %3.0f %4.1f    %5.1f   %5.3f       %s" %
          (t*60.0, Lambda, phi, dauer, h, pfadbreite, A, typ))

# Main
# ====
    
# für jede Minute 1h vor bis 2 h nach der Referenzstunde 10.00 h UT
dt = 0  # Abweichung TD-UT. Ich rechne hier in TD (34)
print("Zeit  Lambda    phi      Dauer  Hohe  Pfad   Monddurchm. Typ")
print("      °,<0:O    °,+:N     s     °    (km)    ")
print("---------------------------------------------------------------")
for i in range(-120, 120, 1): # von 2h vor bis 2h nach der Referenzzeit in 10-min-Schritten
    t = i/60.0  # Zeit in Stunden
    calcEclipse(t)
print("Maximum:")
calcEclipse(21/60.0+8/3600.0)
# Prompt am Ende
ende = input("Fertig.")
Uwe Pilz, Februar 2025