Liebe Beobachter, im Heft 88 habe ich ein Programm besprochen, mit welchem man besondere Schattenwürfe berechnen kann. Hierzu wird der Sonnenort benötigt, und der Mondort. Beides habe ich aus dem "Meeus", 2. Auflage.
Anbei das Programm:
Code:
# hesiodus.py
# Wie in Meeus: AA werde alle Winkel in Grad geführt und erst in den
# Winkelfunktionen #umgewandelt.Dies erleichtert die Tests.
from math import *
import math
toRad = math.pi / 180
# Tab. 47.A, 339
LDterms = [
[0, 0, 1, 0, 6288774, -20905355],
[2, 0, -1, 0, 1274027, -3699111],
[2, 0, 0, 0, 658314, -2955968],
[0, 0, 2, 0, 213618, -569925],
[0, 1, 0, 0, -185116, 48888],
[0, 0, 0, 2, -114332, -3149],
[2, 0, -2, 0, 58793, 246158],
[2, -1, -1, 0, 57066, -152138],
[2, 0, 1, 0, 53322, -170733],
[2, -1, 0, 0, 45758, -204586],
[0, 1, -1, 0, -40923, -129620],
[1, 0, 0, 0, -34720, 108743],
[0, 1, 1, 0, -30383, 104755],
[2, 0, 0, -2, 15327, 10321],
[0, 0, 1, 2, -12528, 0],
[0, 0, 1, -2, 10980, 79661],
[4, 0, -1, 0, 10675, -34782],
[0, 0, 3, 0, 10034, -23210],
[4, 0, -2, 0, 8548, -21636],
[2, 1, -1, 0, -7888, 24208],
[2, 1, 0, 0, -6766, 30824],
[1, 0, -1, 0, -5163, -8379],
[1, 1, 0, 0, 4987, -16675],
[2, -1, 1, 0, 4036, -12831],
[2, 0, 2, 0, 3994, -10445],
[4, 0, 0, 0, 3861, -11650],
[2, 0, -3, 0, 3665, 14403],
[0, 1, -2, 0, -2689, -7003],
[2, 0, -1, 2, -2602, 0],
[2, -1, -2, 0, 2390, 10056],
[1, 0, 1, 0, -2348, 6322],
[2, -2, 0, 0, 2236, -9884],
[0, 1, 2, 0, -2120, 5751],
[0, 2, 0, 0, -2069, 0],
[2, -2, -1, 0, 2048, -4950],
[2, 0, 1, -2, -1773, 4130],
[2, 0, 0, 2, -1595, 0],
[4, -1, -1, 0, 1215, -3958],
[0, 0, 2, 2, -1110, 0],
[3, 0, -1, 0, -892, 3258],
[2, 1, 1, 0, -810, 2616],
[4, -1, -2, 0, 759, -1897],
[0, 2, -1, 0, -713, -2117],
[2, 2, -1, 0, -700, 2354],
[2, 1, -2, 0, 691, 0],
[2, -1, 0, -2, 596, 0],
[4, 0, 1, 0, 549, -1423],
[0, 0, 4, 0, 537, -1117],
[4, -1, 0, 0, 520, -1571],
[1, 0, -2, 0, -487, -1739],
[2, 1, 0, -2, -399, 0],
[0, 0, 2, -2, -381, -4421],
[1, 1, 1, 0, 351, 0],
[3, 0, -2, 0, -340, 0],
[4, 0, -3, 0, 330, 0],
[2, -1, 2, 0, 327, 0],
[0, 2, 1, 0, -323, 1165],
[1, 1, -1, 0, 299, 0],
[2, 0, 3, 0, 294, 0],
[2, 0, -1, -2, 0, 8752]]
# Tab. 47b, 341
Lterms = [
[0, 0, 0, 1, 5128122],
[0, 0, 1, 1, 280602],
[0, 0, 1, -1, 277693],
[2, 0, 0, -1, 173237],
[2, 0, -1, 1, 55413],
[2, 0, -1, -1, 46271],
[2, 0, 0, 1, 32573],
[0, 0, 2, 1, 17198],
[2, 0, 1, -1, 9266],
[0, 0, 2, -1, 8822],
[2, -1, 0, -1, 8216],
[2, 0, -2, -1, 4324],
[2, 0, 1, 1, 4200],
[2, 1, 0, -1, -3359],
[2, -1, -1, 1, 2463],
[2, -1, 0, 1, 2211],
[2, -1, -1, -1, 2065],
[0, 1, -1, -1, -1870],
[4, 0, -1, -1, 1828],
[0, 1, 0, 1, -1794],
[0, 0, 0, 3, -1749],
[0, 1, -1, 1, -1565],
[1, 0, 0, 1, -1491],
[0, 1, 1, 1, -1475],
[0, 1, 1, -1, -1410],
[0, 1, 0, -1, -1344],
[1, 0, 0, -1, -1335],
[0, 0, 3, 1, 1107],
[4, 0, 0, -1, 1021],
[4, 0, -1, 1, 833],
[0, 0, 1, -3, 777],
[4, 0, -2, 1, 671],
[2, 0, 0, -3, 607],
[2, 0, 2, -1, 596],
[2, -1, 1, -1, 491],
[2, 0, -2, 1, -451],
[0, 0, 3, -1, 439],
[2, 0, 2, 1, 422],
[2, 0, -3, -1, 421],
[2, 1, -1, 1, -366],
[2, 1, 0, 1, -351],
[4, 0, 0, 1, 331],
[2, -1, 1, 1, 315],
[2, -2, 0, -1, 302],
[0, 0, 1, 3, -283],
[2, 1, 1, -1, -229],
[1, 1, 0, -1, 223],
[1, 1, 0, 1, 223],
[0, 1, -2, -1, -220],
[2, 1, -1, -1, -220],
[1, 0, 1, 1, -185],
[2, -1, -2, -1, 181],
[0, 1, 2, 1, -177],
[4, 0, -2, -1, 176],
[4, -1, -1, -1, 166],
[1, 0, 1, -1, -164],
[4, 0, 1, -1, 132],
[1, 0, -1, -1, -119],
[4, -1, 0, -1, 115],
[2, -2, 0, 1, 107]]
def getJD(y, m, d) :
if (m < 3):
y = y - 1
m = m + 12
a = (y // 100)
b = 2 - a + (a // 4.0)
jd = floor(365.25 * (y + 4716))
jd = jd + int(30.6001 * (m + 1)) + d + b - 1524.5
return jd
def getKalender(myjd): # vgl. Wikipedia
tag = 0 # eingegliedert, für statische Methode
tag = 0
monat = 0
jahr = 0
stunde = 0
minute = 0
z = int((myjd + 0.5)) # Z = Int(JD + 0,5)
f = (myjd + 0.5) - z # F = Frac(JD + 0,5)
if (z < 2299161):
a = z # wenn Z < 2299161 dann A = Z Ergebnis julianisch
else: # gregorianisch:
g = (int)((z - 1867216.25) / 36524.25) # g = Int((Z - 1867216,25) / 36524,25)
a = z + 1 + g - (int)(g / 4) # A = Z + 1 + g - Int(g/4)
b = a + 1524 # B = A + 1524
c = int((b - 122.1) / 365.25) # C = Int((B-122,1) / 365,25)
d = int(365.25 * c) # D = Int(365,25 * C)
e = int((b - d) / 30.6001) # E = Int((B-D) / 30,6001)
tagDouble = b - d - int(30.6001 * e) + f # Tag = B - D - Int(30,6001*E) + F # Tag, inklusive Tagesbruchteil
tag = int(tagDouble)
if (e < 14) :
monat = e - 1
else :
monat = e - 13 # wenn E<14 dann Monat = E - 1 sonst Monat = E - 13
if (monat > 2):
jahr = c - 4716
else :
jahr = c - 4715 # wenn Monat>2 dann Jahr = C - 4716 sonst Jahr = C - 4715
bruchteil = tagDouble - tag
stunde = int(24 * bruchteil)
minute = int(60 * (24 * bruchteil - stunde))
return jahr, monat, tag, stunde, minute
def mod360(x) :
x = 360 * (x / 360 - int(x / 360))
if (x < 0):
x = x + 360
return x
def moon(JD) :
T = (JD - 2451545) / 36525
T2 = T * T
T3 = T2 * T
T4 = T3 * T
# mittl. Länge samt Lichtzeit
Ls = 218.3164477 + 481267.88123421 * T - 0.0015786 * T2 + T3 / 538841 - T4 / 65194000
Ls = mod360(Ls)
# mittlere Elongation
D = 297.8501921 + 445267.1114034 * T - 0.0018819 * T2 + T3 / 545868 - T4 / 113065000
D = mod360(D)
# mittl. Anomalie der Sonne
M = 357.5291092 + 35999.0502909 * T - 0.0001536 * T2 + T3 / 24490000
M = mod360(M)
# mittl. ANomalie des Mondes
Ms = 134.9633964 + 477198.8675055 * T + 0.0087414 * T2 + T3 / 69699 - T4 / 14712000
Ms = mod360(Ms)
# Argument der Breite
F = 93.2720950 + 483202.0175233 * T - 0.0036539 * T2 - T3 / 3526000 + T4 / 863310000
F = mod360(F)
A1 = 119.75 + 131.849 * T
A1 = mod360(A1)
A2 = 53.09 + 479264.290 * T
A2 = mod360(A2)
A3 = 313.45 + 481266.484 * T
A3 = mod360(A3)
E = 1 - 0.002516 * T - 0.0000074 * T2
E2 = E * E
# Summe der periodischen Terme für Länge und Abstand
SL = 0
Sr = 0
for i in range(60):
a = LDterms[i][0] * D + LDterms[i][1] * M + LDterms[i][2] * Ms + LDterms[i][3] * F
a = a * math.pi / 180 # in Radiant
if ((LDterms[i][1] == 2) or (LDterms[i][1] == -2)) :
SL += E2 * LDterms[i][4] * sin(a)
Sr += E2 * LDterms[i][5] * cos(a)
elif ((LDterms[i][1] == 1) or (LDterms[i][1] == -1)) :
SL += E * LDterms[i][4] * sin(a)
Sr += E * LDterms[i][5] * cos(a)
else :
SL += LDterms[i][4] * sin(a)
Sr += LDterms[i][5] * cos(a)
# Summe der periodischen Terme für die Breite
Sb = 0
for i in range(60):
a = Lterms[i][0] * D + Lterms[i][1] * M + Lterms[i][2] * Ms + Lterms[i][3] * F
a = a * math.pi / 180 # in Radiant
if ((Lterms[i][1] == 2) or (Lterms[i][1] == -2)) :
Sb += E2 * Lterms[i][4] * sin(a)
elif ((Lterms[i][1] == 1) or (Lterms[i][1] == -1)):
Sb += E * Lterms[i][4] * sin(a)
else:
Sb += Lterms[i][4] * sin(a)
# Terme für Jupiter, Venus und Erdabplattung
SL += 3958 * sin(toRad * A1) + 1962 * sin(toRad * (Ls - F)) + 318 * sin(toRad * A2)
Sb += -2235 * sin(toRad * Ls) + 382 * sin(toRad * A3) + 175 * sin(toRad * (A1 - F)) + 175 * sin(toRad * (A1 + F)) + 127 * sin(toRad * (Ls - Ms)) - 115 * sin(toRad * (Ls + Ms))
# die Koordinaten selbst
delta = (385000.56 + Sr / 1000.0) # / AU
lambdaM = mod360(Ls + SL / 1000000.0)
beta = Sb / 1000000.0
pi = asin(6378.14 / delta) / toRad
return lambdaM, beta, delta, pi, Ls, D, M, Ms, F, E
# Sonne: Kap. 25, 163
def sonne(JD) :
T = (JD - 2451545) / 36525
T2 = T * T
T3 = T2 * T
T4 = T3 * T
T5 = T4 * T
Lo = 280.4664567 + 36000.76982779 * T + 0.003032028 * T2 + T3 / 49931 - T4 / 15299 - T5 / 1988000
Lo = mod360(Lo)
M = 357.52911 + 35999.050 * T - 0.0001537 * T2
M = mod360(M)
e = 0.016708617 - 0.000042037 * T - 0.0000001237 * T2
C = (1.9146 - 0.004817 * T - 0.000014 * T2) * sin(toRad * M)
C = C + (0.019993 - 0.000101 * T) * sin(2 * toRad * M)
C = C + 0.00029 * sin(3 * toRad * M)
C = mod360(C)
trueLon = Lo + C # wahre Länge
trueLon = mod360(trueLon)
trueAnom = M + C
R = (1.000001018 * (1.0 - e * e)) / (1.0 + e * cos(toRad * trueAnom))
Omega = 125.0445479 - 1934.1362891 * T + 0.0020754 * T2 + T3 / 467441.0 - T4 / 60616000
# scheinbare Länge
lambdaS = trueLon - 0.00569 - 0.00478 * sin(toRad * Omega)
lambdaS = mod360(lambdaS)
return lambdaS, R, Omega, Lo
# Nutation und Schiefe der Ekliptik: Kap. 22 143
def nutation(JD, Omega, Lo, Ls) :
T = (JD - 2451545) / 36525
T2 = T * T
T3 = T2 * T
deltaPsi = -17.20 * sin(toRad * Omega) - 1.32 * sin(2 * toRad * Lo) - 0.23 * sin(2 * toRad * Ls) + 0.21 * sin(2 * toRad * Omega)
deltaPsi /= 3600
deltaEps = 9.2 * cos(toRad * Omega) + 0.57 * cos(2 * toRad * Lo) + 0.10 * cos(2 * toRad * Ls) - 0.09 * cos(2 * toRad * Omega)
deltaEps /= 3600
epsilon = 23.43929111 + (-46.8150 * T - 0.00059 * T2 + 0.001813 * T3) / 3600 #22.2
epsilon = 23.4392911 * 3600 - 46.8150 * T - 0.00059 * T2 + 0.001813 * T3
epsilon = epsilon / 3600
epsilon = epsilon + deltaEps
return deltaPsi, deltaEps, epsilon
def colongitude(JD) :
lambdaM, beta, delta, pi, Ls, D, M, Ms, F, E = moon(JD)
lambdaS, R, Omega, Lo = sonne(JD)
deltaPsi, deltaEps, epsilon = nutation(JD, Omega, Lo, Ls)
T = (JD - 2451545) / 36525
I = 1.54242 # Neigung des Mondäquators
# die Libration 53.1
W = lambdaM - Omega
W = mod360(W)
A = atan2((sin(toRad * W) * cos(toRad * beta) * cos(toRad * I) - sin(toRad * beta) * sin(toRad * I)), (cos(toRad * W) * cos(toRad * beta)))
A = A / toRad
A = mod360(A)
ls = A - F
sbs = -sin(toRad * W) * cos(toRad * beta) * sin(toRad * I) - sin(toRad * beta) * cos(toRad * I)
bs = asin(sbs) / toRad
K1 = 119.75 + 131.849 * T
K2 = 72.56 + 20.186 * T
rho = (-0.02752 * cos(toRad * Ms) - 0.02245 * sin(toRad * F) + 0.00684 * cos(toRad * Ms - 2 * toRad * F) + -0.00293 * cos(2 * toRad * F) - 0.00085 * cos(2 * toRad * F - 2 * toRad * D) - 0.00054 * cos(toRad * Ms - 2 * toRad * D) + -0.00020 * sin(toRad * Ms + toRad * F) - 0.00020 * cos(toRad * Ms + 2 * toRad * F) + -0.00020 * cos(toRad * Ms - toRad * F) + 0.00014 * cos(toRad * Ms + 2 * toRad * F - 2 * toRad * D))
sigma = (-0.02816 * sin(toRad * Ms) + 0.02244 * cos(toRad * F) - 0.00682 * sin(toRad * Ms - 2 * toRad * F) + -0.00279 * sin(2 * toRad * F) - 0.00083 * sin(2 * toRad * F - 2 * toRad * D) + 0.00069 * sin(toRad * Ms - 2 * toRad * D) + 0.00040 * cos(toRad * Ms + toRad * F) - 0.00025 * sin(2 * toRad * Ms) + -0.00023 * sin(toRad * Ms + 2 * toRad * F) + 0.00020 * cos(toRad * Ms - toRad * F) + 0.00019 * sin(toRad * Ms - toRad * F) + 0.00013 * sin(toRad * Ms + 2 * toRad * F - 2 * toRad * D) + -0.00010 * cos(toRad * Ms - 3 * toRad * F))
tau = (0.02520 * E * sin(toRad * M) + 0.00473 * sin(2 * toRad * Ms - 2 * toRad * F) - 0.00467 * sin(toRad * Ms) + 0.00396 * sin(toRad * K1) + 0.00276 * sin(2 * toRad * Ms - 2 * toRad * D) + 0.00196 * sin(toRad * Omega) + -0.00183 * cos(toRad * Ms - toRad * F) + 0.00115 * sin(toRad * Ms - 2 * toRad * D) + -0.00096 * sin(toRad * Ms - toRad * D) + 0.00046 * sin(2 * toRad * F - 2 * toRad * D) + -0.00039 * sin(toRad * Ms - toRad * F) - 0.00032 * sin(toRad * Ms - toRad * M - toRad * D) + 0.00027 * sin(2 * toRad * Ms - toRad * M - 2 * toRad * D) + 0.00023 * sin(toRad * K2) + -0.00014 * sin(2 * toRad * D) + 0.00014 * cos(2 * toRad * Ms - 2 * toRad * F) + -0.00012 * sin(toRad * Ms - 2 * toRad * F) - 0.00012 * sin(2 * toRad * Ms) + 0.00011 * sin(2 * toRad * Ms - 2 * toRad * M - 2 * toRad * D))
lss = -tau + (rho * cos(toRad * A) + sigma * sin(toRad * A)) * tan(toRad * bs)
bss = sigma * cos(toRad * A) - rho * sin(toRad * A)
l = ls + lss
b = bs + bss
V = Omega + deltaPsi + sigma / sin(toRad * I)
X = sin(toRad * I + toRad * rho) * sin(toRad * V)
Y = sin(toRad * I + toRad * rho) * cos(toRad * V) * cos(toRad * epsilon) - cos(toRad * I + toRad * rho) * sin(toRad * epsilon)
omega = atan2(X, Y) / toRad
alpha = atan2(sin(toRad * lambdaM) * cos(toRad * epsilon) - tan(toRad * beta) * sin(toRad * epsilon), cos(toRad * lambdaM)) / toRad
sP = sqrt(X * X + Y * Y) * cos(toRad * alpha - toRad * omega)
sP = sP / cos(toRad * b)
P = asin(sP) / toRad
# die selenographische Sonnenposition
AE = 149597870.7 # km
lambdaH = lambdaS + 180 + delta / (AE * R) * 180 / math.pi * cos(toRad * beta) * sin(toRad * (lambdaS - lambdaM))
betaH = delta / (AE * R) * beta
# 53.1 neu
Wo = lambdaH - deltaPsi - Omega
Wo = mod360(Wo)
Ao = atan2((sin(toRad * Wo) * cos(toRad * betaH) * cos(toRad * I) - sin(toRad * betaH) * sin(toRad * I)), (cos(toRad * Wo) * cos(toRad * betaH)))
Ao = Ao / toRad
Ao = mod360(Ao)
lso = Ao - F
sbso = -sin(toRad * Wo) * cos(toRad * betaH) * sin(toRad * I) - sin(toRad * betaH) * cos(toRad * I)
bso = asin(sbso) / toRad
lsso = -tau + (rho * cos(toRad * Ao) + sigma * sin(toRad * Ao)) * tan(toRad * bso)
bsso = sigma * cos(toRad * Ao) - rho * sin(toRad * Ao)
lo = lso + lsso
bo = bso + bsso
co = 90 - lo
co = mod360(co)
return co, b
# Hauptprogramm
j = 2020 # Startjahr
J = 1 # Anzahl Jahre
coZiel = 18.40 # was wird gesucht? Hesiodus: coZiel=18,40°
vorigeCo = 7e77
JD = getJD(j, 1, 1.0)
deltaJD = 1 / 24.0 # 1 h
for i in range(24 * 366 * J):
co,b = colongitude(JD)
if ((co > coZiel) and (vorigeCo < coZiel)):
JDHesiodus=JD-1/24*(co-coZiel)/(co-vorigeCo);
jahr, monat, tag, stunde, minute=getKalender(JDHesiodus)
print(tag,monat,jahr," \t",stunde," ",minute," UT. b=",b)
JD = JD+deltaJD
vorigeCo = co
Registrieren
Anmelden
FAQ
Suche
