leider haben mich Lesen, Ausprobieren, Tutorials durcharbeiten und Suchfunktion sowie Google nicht zu einer Lösung geführt.
Folgendes Problem möchte ich lösen: Anhand der üblicherweise genutzten Formeln in der Astronomie möchte ich den aktuellen Sonnenstand (Höhenwinkel und Azimut) berechnen. Dazu muss aus dem aktuellen Datum (Zeit und Tag) ein Julianisches Datum berechnet werden. Mit diesem kann dann die aktuelle Sonnenposition erechnet werden. Grundsätzlich funktioniert meine Funktion auch, aber leider sind die Ergebnisse unbrauchbar. Und ich finde meinen Gedankenfehler nicht.
Code: Alles auswählen
def calculateSunPosition(lonGrad1, latGrad1):
lonGrad = float(lonGrad1)
latGrad = float(latGrad1)
#Aus dem Kalenderdatum ein Julianisches Datum erzeugen
lt_now = time.gmtime()
#Julianisches Datum JD und JD0
Jahr = lt_now[0]
Monat = lt_now[1]
if Monat <= 2:
Jahr = Jahr - 1
Monat = Monat + 12
Stunden = int(time.strftime("%H",time.localtime())) / 24 + lt_now[4] / 1440 + lt_now[5] / 86400
Jahrhunderte = lt_now[0] / 100
A = 2 - Jahrhunderte + Jahrhunderte / 4
B = 365.25 * (Jahr + 4716)
C = (30.6001 * (Monat + 1))
JD0 = (A + B + C + lt_now[2] - 1524.5)
JD = (A + B + C + lt_now[2] + Stunden -1524.5)
T0 = (JD0 - 2451545.0) / 36525.0
T = Stunden * 24
#JD0 = JD
#Berechnung von n
n = JD - 2451545.0
#Berechnung der ekliptikalen Länge L der Sonne
L = Degr(280.460*2*math.pi/360 + 0.9856474*2*math.pi/360*n, 2*math.pi)
#Berechnung der ekliptikalen Länge gGamma inkl. Bahnelliptizität, mittlerer Anomalie g und Exzentrizität e=0,0167
g = Degr(357.528*2*math.pi/360 + 0.9856003*2*math.pi/360*n, 2*math.pi)
gGamma = L + 1.915*math.sin(g) + 0.01997*math.sin(2*g)
#Berechnung der Rektaszension alpha mit der Schiefe der Ekliptik epsilon
epsilon = 23.439*2*math.pi/360 - 0.0000004*2*math.pi/360*n
alpha = math.atan(math.cos(epsilon)*math.sin(gGamma)/math.cos(gGamma))
if math.cos(gGamma) < 0:
alpha = alpha + math.pi #180.0
else:
alpha = alpha
#Berechnung der Deklination delta
delta = math.asin((math.sin(epsilon)*math.sin(gGamma)))
#Berechnung von Theta Greenwich
thetaG = Degr(6.697376+2400.05134*T0+1.002738*JD*24, 24)
theta = thetaG * 0.2617993878 + lonGrad*2*math.pi/360
tau = theta - alpha
#Berechnung des Azimuts a
phi = latGrad*2*math.pi/360
#Nenner Azimut na
na = math.cos(tau)*math.sin(phi)-math.tan(delta)*math.cos(lonGrad*2*math.pi/360)
rz = math.atan(math.sin(tau) / na) + math.pi
if na < 0:
rz = rz + math.pi
a = Degr(rz, 2*math.pi)
#Berechnung des Höhenwinkels h
h = math.asin(math.cos(delta)*math.cos(tau)*math.cos(phi)+math.sin(delta)*math.sin(phi))
#Korrekturwert für Refraktion
kr = 0.0002967059728 / math.tan(h + 0.179768913 / (h + 0.08918632478))
h = h + kr
print(a*360/(2*math.pi))
print(h*360/(2*math.pi))
Schöne Grüße
Ingo