VdS-Journal 97: Zentrallinie einer Sonnenfinsternis in der Weltkarte

Die Berechnung erfolgt minütlich. Vor Benutzung müssen die sog. Bessel-Elemente der gewünschten Finsternis eingesetzt werden.
Das Modul worldmap2.py muss im Programmverzeichnise liegen.
# sofiBessel2Wordlmap : Sonnenfinsternisberechnung mit den Besselschen Elementen:
# Verlauf der Zentrallinie. In die Weltkarte geplottet
# für Heft 97

# 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 *
from worldmap2 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, screen, W):
    # 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)
    # Protokollausgabe nach Meeus, S. 13
    #print(X, Y, D, M, L2, XS, YS, w, p, b, c, y1, b1, b2, B)
    #print(H, phi, Lambda, L2S, a, n, dauer, h, K2, pfadbreite)
    #print(L1S, A,)
    # wirklich ausgeben (formatiert): 
    #   Zeit (Minuten seit Referenzstunde)
    #   geographische Koordinaten in Grad. West und Nord sind positiv
    #   Sonnenhöhe in Grad
    #   Breite des Totalitätspfades
    #   Größe der Mondscheibe in scheinbarem Sonnendurchmesser
    #   ringförmig oder total
    if Lambda>360:
        Lambda=Lambda-360
    if (Lambda>180):
        Lambda=Lambda-360
    # Tabellenausgabe ist auskommentiert.
    # print(t*60.0, Lambda, phi, dauer, h, pfadbreite, A, typ)
    # Druck in die Weltkarte
    worldpoint(Lambda/180.0*pi,phi/180.0*pi,4,"blue", W, screen, 1)
    

# Main
# ====
screen = Screen()
print("Weltkarte: 10 Sekunden warten")
w=1000 # Bildbreite
worldmap(1, "green", w, screen)
    
# für jede Minute 1h vor bis 2 h nach der Referenzstunde 
dt = 0  # Abweichung TD-UT. Ich rechne hier in TD (34)
# keine Überschrit in der Weltkarten-Variante.
# print("Zeit, Lambda/°(<0:Ost)   phi(pos:Nord)   Dauer/s          Hohe/°              Pfadbreite/km    rel.Monddurchm.  Typ")
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, screen, w)
# Keine Maximumsberechnung für die Kartenansicht
# print("Maximum:") 
# calcEclipse(21/60.0+8/3600.0, screen)
    

# formatiert drucken ohne Zeilenwechsel
# print("%5.3f " % (j * dt),end="")

# Prompt am Ende
update()
ende = input("Fertig.")

# ===================
Uwe Pilz, Februar 2025