VdS-Journal 95: Lokaler Verlauf einer Finsternis

Die Berechnung erfolgt minütlich. Vor Benutzung müssen die sog. Bessel-Elemente der gewünschten Finsternis eingesetzt werden. Die Datumsangabe muss aktiviert ("einkommentiert") werden.
Geographische Koordinaten werden astronomisch richtig ausgegebe, also Ost = negativ. Das ist umgekehrt zur Google-Maps-Anzeige.
# sofiBessel3 : lokaler Verlauf einer Sonnenfinsternis
# ohne Berücksichtigung der lokalen Sonnenhöhe

from math import *

# Kennwerte der Berechnung - hier ändern.

# geographische Koordinaten, Ost =negativ (umgekehrt wie Google Maps)
# kann aus Maps kopiert werden, dann nur noch das Vz ändern.
P=39.57627485850586
L=-2.6532853318252094

# 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


dt=61 # Abweichung der Weltzeit von der dynamischen Zeit
'''
1970      40.18
1980      50.54
1990      56.86
2000      63.83
2005      64.69
2007      65.15
2016      68.109
'''

# um die Schreibarbeit der kopierten Koordinaten zu verringern dort nur 1 Buchstabe
Lambda= L
phi=P
hoehe=0 # wird nicht berücksichtigt. Das Programm kann das aber.

# 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)


# Ergebnisse der Fundamentalebene für einen Zeitpunkt
def magnitudePoswinkel(t):
    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

    # Koordinaten bezüglich der Fundamentalebene
    H = M - Lambda - 0.00417807 * dt
    xi = pCosPS * Sin(H)
    eta = pSinPS * Cos(D) - pCosPS * Cos(H) * Sin(D)
    zeta = pSinPS * Sin(D) + pCosPS * Cos(H) * Cos(D)
    # stündliche Änderungen
    xiS = 0.01745329 * u[1] * pCosPS * Cos(H)
    etaS = 0.01745329 * (u[1] * xi * Sin(D) - zeta * d[1])

    U = X - xi
    V = Y - eta
    a = XS - xiS
    b = YS - etaS
    L1S = L1 - zeta * tanf1
    L2S = L2 - zeta * tanf2
    n2 = a * a + b * b
    return a, b, U, V, n2, L1S, L2S, D, H


# Ergebnisse für die Ausgabe
def ergebnisberechnung():
    # Zeit des lokalen Maximums
    TD = T0 + t
    # Magnitude G
    m = sqrt(U * U + V * V)
    G = (L1S - m) / (L1S + L2S)
    # Durchmesserverhältnis
    A = (L1S - L2S) / (L1S + L2S)
    # Positionswinklel zur maximalen Verfinsterung, vom Nordpol der Sonne
    Pm = Atan2(U , V)
    if Pm < 0:
        Pm = Pm + 360.0
    # Positionswinkel Zm bezogen auf die Zenit-Richtung
    sinH = Sin(D) * Sin(phi) + Cos(D) * Cos(phi) * Cos(H)
    h = Asin(sinH)
    sinq = (Cos(phi) * Sin(H)) / Cos(h)
    q = Asin(sinq)
    Zm = Pm - q
    # Positionswinkel
    P = Atan2(U, V)
    if P < 0:
        P = P + 360
    UT = TD - dt / 3600
    UTh = int(UT)
    UTm0 = 60 * (UT - UTh)
    UTm = int(UTm0)
    UTs = int(60 * (UTm0 - UTm) + 0.5)
    return UTh, UTm, UTs, G, A, Zm, P


# Hauptprogramm
t = 0

# rechtwinklige geozentrische Koordinaten
U = Atan(0.99664719 * Tan(phi))
pSinPS = 0.99664719 * Sin(U) + hoehe / 6378140 * Sin(phi)
pCosPS = Cos(U) + hoehe / 6378140 * Cos(phi)

# 1) Ergebnisse zum Maximums-Zeitpunkt
for i in range(5):  # 5 Iterationen genügen
    a, b, U, V, n2, L1S, L2S, D, H = magnitudePoswinkel(t)
    # Zeitpunkt der maximalen Finsternis per Iteration
    tm = -(U * a + V * b) / n2
    t = t + tm
    UTh, UTm, UTs, G, A, Zm, P = ergebnisberechnung()
print("Uhrzeit des Maximums:", UTh, "h", UTm, "m", UTs, "s UT")  # dynamische Zeit
print("Magnitude: ", int(1000 * G) / 10, "%")
print("Verhältnis Durchmesser Mond/Sonne: ", int(1000 * A) / 1000)
print("Positionswinkel, bezogen auf Nord: ", int(P+0.5), "Grad")
print("Positionswinkel, bezogen auf Zenit: ", int(Zm+0.5), "Grad")

# 2) Ergebnisse aller 10 Minuten
tmax = t
t = tmax - 3 
print("UT             Magnitude Nord-  Zenit- Posw.")
for i in range(480):  # 3 h aller Minuten
    a, b, U, V, n2, L1S, L2S, D, H = magnitudePoswinkel(t)
    UTh, UTm, UTs, G, A, Zm, P = ergebnisberechnung()
    if G >= 0:
        print("%2.0f h %02.0f m %02.0f s  %5.1f%%     %3.0f° %3.0f°" % (UTh, UTm, UTs, 100*G, P+0.5, Zm+0.5))
        # print(UTh,"h", UTm,"m", UTs,"s  ", int(1000 * G) / 10, "%  ", int(P+0.5),"°  ", int(Zm+0.5), "°")
    t = t + 1 / 60
Uwe Pilz, Februar 2025