VdS-Journal 99: Globaler Verlauf einer Sonnenfinsternis

Die Berechnung erfolgt dadurch, das engmaschig der lokale Verlauf analysiert wird. Im Verzeichnus muss die Datei besselworldmap2.py vorhanden sein.
# sofBessel3Worldmap : Sonnenfinsternisberechnung mit den Besselschen Elementen:
# Weltkarte

from math import *
from worldmap2 import *

# Kennwerte der Berechnung - hier ändern.

# Finsternis vom 10. Mai 1994, Maximum  17:11:27
x = [-0.173367, 0.4990629, 0.0000296, -0.00000563]  
y = [0.383484, 0.0869393, -0.0001183, -0.00000092]
d = [17.6861305, 0.0106420, -0.0000040, 0.0000000]
l1 = [0.566906, 0.0000318, 0.0000098, 0.0000000]
l2 = [0.020679, -0.0000317, -0.0000097, 0.0000000]
u = [75.909233, 15.001620, 0.000000, 0.000000]
tanf1 = 0.0046308
tanf2 = 0.0046077
yr=1994
mnt=5
day=10
T0=17 # h UT


# 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 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
    # ####### print("T",T)
    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
    # Omega=125.04-1934.136*T
    # Lambda=L-0.00569-0.00478*Sin(Omega)
    eps=schiefeDerEkliptik(T)
    # ##### print("L",L,eps)

    alpha=Atan2(Cos(eps)*Sin(L), Cos(L))
    delta=Asin(Sin(eps)*Sin(L))
    return(alpha, delta)

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, pSinPS, pCosPS, Lambda):
    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 # dt=0 für globale Betrachtung. - 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(t, U, V, L1S, L2S, D, H, phi):
    # 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=Atan(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=Atan(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
# ================

screen = Screen()
w=1200 # halbe Bildbreite. w=1200 und bruchteil=4 gibt schöne Bilder, rechnet 2-3 min
bruchteil=4 # Koordinaten feiner als 1 Grad: 2 =1/2 3=1/3 usw.  Nur Ganzzahlen
worldmap(1, "green", w, screen)
phi=90
for j in range(180*bruchteil): # alle Breitengrade
    Lambda=-180
    for k in range (360*bruchteil): # alle Längengrade
        t=0 
        # rechtwinklige geozentrische Koordinaten
        U = Atan(0.99664719*Tan(phi))
        pSinPS=0.99664719 * Sin(U)  # weltweit: h=0. +hoehe/6378140*Sin(phi)
        pCosPS=Cos(U)               # + hoehe/6378140 * Cos(phi)
        # Ergebnisse zum Maximums-Zeitpunkt per Iteration
        for i in range (5): # 5 Iterationen genügen
            a,b, U, V, n2, L1S, L2S, D, H = magnitudePoswinkel(t, pSinPS, pCosPS, Lambda)
            tm=-(U*a+V*b)/n2
            t=t+tm
        UTh, UTm, UTs, G, A, Zm, P= ergebnisberechnung(t,  U, V, L1S, L2S, D, H, phi)
        jd=JD(yr, mnt, day+(T0+t)/24.0)
        alpha,delta=sonnenkoordinaten(jd)
        theta0=sternzeit(jd) # für Greenwich
        h=shoehe(theta0, alpha, delta, Lambda, phi)
        if (h>0):
            if (G>0) and (h<1):
                worldpoint(Lambda/180.0*pi,phi/180.0*pi,3,"magenta", w, screen, 1)
            if (G>0) and (G<0.05):
                worldpoint(Lambda/180.0*pi,phi/180.0*pi,2,"blue", w, screen, 1)
            if (G>0.49) and (G<0.51): 
                worldpoint(Lambda/180.0*pi,phi/180.0*pi,2,"blue", w, screen, 1)
            if (G>0.96) and (G<=1):
                worldpoint(Lambda/180.0*pi,phi/180.0*pi,3,"red", w, screen, 1)
        Lambda=Lambda+1.0/bruchteil
    phi=phi-1.0/bruchteil
    # print (j, phi, Lambda, day+(T0+t)/24.0, h)
    update()
ende = input("Fertig.")
Uwe Pilz, Februar 2025