VdS-Journal 95: Lokaler Verlauf einer Sonnenfinsternis, ausführlich

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 ausgegeben, also Ost = negativ. Das ist umgekehrt zur Google-Maps-Anzeige.
# sofiBessel3A : lokaler Verlauf einer Sonnenfinsternis
# mit 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 # dt=TD/ET-UT Abweichung der Weltzeit von der dynamischen Zeit
'''
1900  0 p.D
1950  29.1
1960  33.1
1970  40.18
1980  50.54
1990  56.86
2000  63.83
2005  64.69
2007  65.15
2016  68.109
2100  (240)
'''

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

def to360(w):
    w=w/360
    return 360*(w - int(w))

def To360(w):
    w=to360(w)
    if (w<0):w=w+360
    return w

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

def sternzeit(JD): # Meeus Kap. 12
    # einfache Methode, welche die Zahlengenauigkeit belastet
    T = (JD - 2451545.0) / 36525.0
    theta0 = 280.46061837 + 360.98564736629 * (JD - 2451545.0) +  0.000387933 * T * T - T * T * T / 39710000.0;
    theta0=To360(theta0)
    return theta0 # in Grad

def schiefeDerEkliptik(T): # ohne Nutation
    eps = 23 + 26 / 60.0 + 21.448 / 3600;
    eps = eps + (-46.8150 * T - 0.00059 * T * T + 0.001813 * T * T * T) / 3600;
    return eps

def sonnenkoordinaten(jd): # einfache Genauigkeit, Meeus Kap 25
    T = (jd - 2451545) / 36525    # Julianisches Jahrhundert
    L0 = To360(280.46646 + T * (36000.76983 + T * 0.0003032))
    M = To360(357.52911 + T * (35999.05029 - 0.0001537 * T))
    # e = 0.016708634 - T * (0.000042037 + 0.0000001267 * T)
    C=(1.914602-0.004817*T-0.000014*T*T)*Sin(M)
    C=C+(0.019993-0.000101*T)*Sin(M+M)
    C=C+0.000289*Sin(M+M+M)
    L=L0+C
    eps=schiefeDerEkliptik(T)
    alpha=Atan2(Cos(eps)*Sin(L), Cos(L))
    delta=Asin(Sin(eps)*Sin(L))
    return(alpha, delta)    

# Sonnenhöhe für einen Zeitpunkt und Ort.
# Erfordert die Sonnenkoordinaten und die Sternzeit
def shoehe(theta0, alpha, delta, long, lat):
    H=To360(theta0-long-alpha)
    # ### print("H", H)
    h=Asin(Sin(lat)*Sin(delta)+Cos(lat)*Cos(delta)*Cos(H))
    return h
           

# 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(jd):
    # 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)
    # Berechnung der Sonnenhöhe ü.H.
    alpha,delta=sonnenkoordinaten(jd)
    theta0=sternzeit(jd) # für Greenwich
    h=shoehe(theta0, alpha, delta, Lambda, phi)
    return UTh, UTm, UTs, G, A, Zm, P, h


# 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
    jd=JD(yr, mnt, day+(T0+t)/24.0)
    UTh, UTm, UTs, G, A, Zm, P, h = ergebnisberechnung(jd)
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")
print("Sonnenhöhe: ", int(h+0.5), "Grad")

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