Berechnung der aktuellen Sonnenposition

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Antworten
sonnenanbeter
User
Beiträge: 7
Registriert: Samstag 9. Januar 2016, 11:47

Hallo liebe Python-Kenner,
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))
Auf einigen Websites (z.B. http://www.sonnenverlauf.de) kann man die Sonnenbahn mitverfolgen, aber meine Ergebnisse weichen in beliebiger Richtung ab. Berechne ich das Julianische Datum falsch? Verwende ich die falschen Zeitfunktionen? Ich stehe völlig auf dem Schlauch und bin über jeden Hinweis dankbar.

Schöne Grüße
Ingo
Zuletzt geändert von Anonymous am Samstag 9. Januar 2016, 12:13, insgesamt 1-mal geändert.
Grund: Quelltext in Python-Codebox-Tags gesetzt.
sonnenanbeter
User
Beiträge: 7
Registriert: Samstag 9. Januar 2016, 11:47

Code: Alles auswählen

def Degr(x, y):
  if x < 0:
    x = x + y
  return((x - math.floor(x/y)*y))
Die Funktion hatte ich noch vergessen zu erwähnen...
BlackJack

@sonnenanbeter: Am besten gehst Du die Rechnung Schritt für Schritt durch und schaust ab wo die Werte nicht mehr stimmen. Notfalls kann man die Rechnung auch mit bekannten Ein- und Ausgabewerten von hinten nach vorne rechnen und mit den Werten im Programm vergleichen.

Die Funktion ist auch ein bisschen sehr lang, die könnte man bestimmt auf mehrere aufteilen, die man dann auch mal einzeln testen könnte. Beispielsweise die Berechnung des julianischen Datums.

Bei den Datumsgeschichten könnte man einiges etwas lesbarer schreiben wenn man nicht magische Indexzahlen in das Ergebnis von `gmtime()` verwenden würde, sondern die Attributnamen. ``lt_now.tm_min`` sagt beispielsweise ein wenig mehr aus als ``lt_now[4]``. ``int(time.strftime("%H", time.localtime()))`` ist etwas umständlich um ``time.localtime().tm_hour`` auszudrücken.

Wobei ich persönlich lieber das `datetime`-Modul verwende wenn es um Datumssachen geht.

Streng genommen darf man die aktuelle Zeit nur einmal ermitteln, also nicht UTC und lokale Zeit getrennt, denn zwischen den beiden Aufrufen vergeht ja Zeit die das Ergebnis verfälschen könnte. Du machst das glücklicherweise (zufällig?) in einer Reihenfolge bei der das nicht so schlimm ist.

Kann man `localtime()` überhaupt einfach so verwenden? Da steckt ja beispielsweise Winter-/Sommerzeit mit drin, diese Stunde Verschiebung verändert doch aber nicht die geographischen Gegebenheiten oder den Sonnenstand‽
sonnenanbeter
User
Beiträge: 7
Registriert: Samstag 9. Januar 2016, 11:47

@BlackJack: Danke für die Hinweise. Auf die Idee mit den Attributnamen bin ich nicht gekommen. Ich habe aber festgestellt, dass meine Berechnung int und float an einigen relevanten Stellen vermischt. Dass muss ich zunächst nochmal sauber trennen und die daraus resultierenden Rechenfehler prüfen/korrigieren. Vermutlich fange ich besser nochmal von vorne an! ;-)
BlackJack

@sonnenanbeter: An das Problem habe ich kurz gedacht, es dann aber verworfen weil Du `print` mit Klammern benutzt und dann ja eigentlich kein solches Problem besteht weil man entweder Python 3 verwendet oder `division` und `print_function` aus dem `__future__`-Modul importiert hat.

Edit: Vielleicht noch ein Hinweis auf den Style Guide for Python Code. Die Namensschreibweisen entsprechen dem nicht wirklich.
sonnenanbeter
User
Beiträge: 7
Registriert: Samstag 9. Januar 2016, 11:47

@BlackJack: Ich versuche in den nächsten Tagen, das Rechenproblem zu ergründen und dann kümmere ich mich um den Style Guide. Aufhübschen werde ich alles, wenn es funktioniert! Danke für die Hinweise und die Denkanstöße!
BlackJack

@sonnenanbeter: Ich würde das nicht als „aufhübschen” bezeichnen. Das ist ja keine Kosmetik sondern hat eine Funktion: Damit andere den Quelltext leichter lesen und verstehen können. Und Du kannst dann auch Quelltext von anderen leichter lesen und verstehen.
Benutzeravatar
pillmuncher
User
Beiträge: 1482
Registriert: Samstag 21. März 2009, 22:59
Wohnort: Pfaffenwinkel

BlackJack hat geschrieben:@sonnenanbeter: ... Damit andere den Quelltext leichter lesen und verstehen können. ...
Einer der von BlackJack genannten anderen bist übrigens du selbst, nächste Woche.
In specifications, Murphy's Law supersedes Ohm's.
sonnenanbeter
User
Beiträge: 7
Registriert: Samstag 9. Januar 2016, 11:47

Ich bin euch noch eine Antwort schuldig. Einmal Wegwerfen und neu anfangen hat geholfen, den Überblick zu gewinnen. Es waren tatsächlich 2 grundlegende Probleme: In den Formel für das Julianische Datum mussten zwei ganzzahlige Werte verwendet werden. Und in der Berechnung des Höhenwinkels müssen die Werte der geografischen Breite nach Osten mit einem Minus vorweg übergeben werden. Das ist aber noch nicht vollständig implementiert. Die Style Guides habe ich versucht zu berücksichtigen, aber bei den 79 Spalten pro Zeile habe ich aufgegeben... Vielen Dank für die Anregungen!!!
Falls das Ergebnis noch jemanden interessiert:

Code: Alles auswählen

def calc_jul_Date():
    jahr = int(time.strftime("%Y", time.gmtime()))
    monat = int(time.strftime("%m", time.gmtime()))
    tag = int(time.strftime("%d", time.gmtime()))
    stunden = int(time.strftime("%H", time.gmtime()))
    minuten = int(time.strftime("%M", time.gmtime()))
    sekunden = int(time.strftime("%S", time.gmtime()))

    if (monat<=2):
        jahr = jahr - 1
        monat = monat + 12

    stundenX = stunden / 24 + minuten / 1440 + sekunden / 86400
    jahrhunderte = int(jahr / 100)
    a = (2 - jahrhunderte + int(jahrhunderte / 4))
    b = int((365.25 * (jahr + 4716)))
    c = int(((30.6001 * (monat + 1))))
    jd0 = (a + b + c + tag - 1524.5)
    jd = (a + b + c + tag + stundenX - 1524.5)
    t0 = (jd0 - 2451545.0) / 36525.0
    t = stundenX * 24
    return jd0, jd, t0, t

def calc_Sun_Position(longGrad, latGrad, julDate0, julDate, time0, time):
    #Definition der Variablen
    laenge = longGrad
    breite = latGrad
    jd0 = julDate0
    jd = julDate
    t0 = time0
    t = time
    
    """Koordinaten Positionsbestimmung
     n,L,g,A
     Deklination, Ekliptik, Nenner, Nenner2, Rektaszension, Refraktion 
     StundenwinkelFP, GStundenwinkelFP, OStundenwinkelFP, OStundenwinkel
     Koordinaten Position"""

    #Definition der Konstanten M_PI
    M_PI = math.pi

    #Berechnung der einzelnen Formelelemente für die Bestimmung der Sonnenposition auf Nord = 0° bezogen
    n = jd - 2451545.0
    l = (280.46 + 0.9856474 * n) % 360
    g = (357.528 + 0.9856003 * n) % 360
    a = l + (math.e * 360 / M_PI) * math.sin(g * M_PI / 180) + (225 * math.e * math.e / M_PI) * math.sin(2 * g * M_PI / 180)
    ekliptik = 23.439 - 0.0000004 * n 
    nenner = math.cos(a * M_PI / 180)
    rektaszension = math.atan((math.cos(ekliptik * M_PI / 180) * math.sin(a * M_PI / 180)) / nenner) * 180 / M_PI
    
    if (nenner<0): 
        rektaszension = rektaszension + 180
 
    deklination = math.asin(math.sin(ekliptik * M_PI / 180) * math.sin(a * M_PI / 180)) * 180 / M_PI
    stundenwinkelFP = (6.697376 + (2400.05134 * t0) + (1.002738 * t)) % 24 
    gStundenwinkelFP = stundenwinkelFP * 15
    oStundenwinkelFP = gStundenwinkelFP + breite
    oStundenwinkel = oStundenwinkelFP - rektaszension
    nenner2 = math.cos(oStundenwinkel * M_PI / 180) * math.sin(laenge * M_PI / 180) - math.tan(deklination * M_PI / 180) * math.cos(laenge * M_PI / 180)
    azimut = math.atan(math.sin(oStundenwinkel * M_PI / 180) / nenner2) * 180 / M_PI
    
    if (nenner2<0):
        azimut = azimut + 180 # korrektur

    azimut = azimut + 180
    azimut = azimut % 360 # korrektur
    
    hoehenwinkel = math.asin(math.cos(deklination * M_PI / 180) * math.cos(oStundenwinkel * M_PI / 180) * math.cos(laenge * M_PI / 180) + math.sin(deklination * M_PI / 180) * math.sin(laenge * M_PI / 180)) * 180 / M_PI
    refraktion = 1.02 / math.tan((hoehenwinkel + (10.3 / (hoehenwinkel + 5.11))) * M_PI / 180) 
    hoehenwinkel= hoehenwinkel + refraktion / 60 
    return azimut, hoehenwinkel
Sirius3
User
Beiträge: 17711
Registriert: Sonntag 21. Oktober 2012, 17:20

@sonnenanbeter: das mehrfache Aufrufen von gmtime ist ein Fehler, wie BlackJack schon bemerkt hatte. Zumal das ermitteln der ganzen Werte ohne den Umweg über Strings viel einfacher geht:

Code: Alles auswählen

jahr, monat, tag, stunden, minuten, sekunden = time.gmtime()[:6]
Warum werden in calc_Sun_Position alle Parameter nochmal umkopiert? Der Schritt ist unnötig, weg damit.
Das Umrechnen von Grad in Radians macht man mit math.radians. Modulo 360 ist überflüssig, weil sin und cos periodische Funktionen sind. Für mich sind im Code zu viele magische Zahlen, die ich als Konstanten schreiben würde.
sonnenanbeter
User
Beiträge: 7
Registriert: Samstag 9. Januar 2016, 11:47

Vielen Dank nochmal für die letzten Hinweise. Das gmtime-Konstrukt hatte ich zuerst nicht verstanden. Die Konstanten habe ich auch noch übernommen. Und die Variablen hatte ich nochmal kopiert, weil IDLE3 zunächst ohne Kopieren rummeckerte. Jetzt geht's! :K
Vielen Dank nochmal an alle. Nach fast 20 Jahren Programmierpause (damals mit Java noch im Studium) muss man sich doch wieder ganz von vorne rankämpfen...
Lobi_Solar
User
Beiträge: 2
Registriert: Dienstag 17. Mai 2022, 10:52

Hallo, habe eine Anwendung für Dein Programm.
Habe es getestet, aber leider kommen keine vernümpftigen Wert heraus.
Wo finde ich die Grundlage Deiner Berechnung?
Kannst Du mir zwischenwerte schicken, damit ich Schritt für Schritt mich durch die Berechnung wühle?

Danke
Gruß
Lobi
__deets__
User
Beiträge: 14493
Registriert: Mittwoch 14. Oktober 2015, 14:29

Ich bezweifele ja, dass der OP sich da noch melden wird. Das ist 6 Jahre her. Aber wenn du zeigst, was du probiert hast (Code-Tags nicht vergessen!), dann kann man dir ggf. helfen.
Lobi_Solar
User
Beiträge: 2
Registriert: Dienstag 17. Mai 2022, 10:52

Hallo zusammen,

ich habe mir vorgenommen, eine drehbare und schwenkbare Solaranlage zu bauen.
Bin 52 und Maschinenbauer.
Mit Python programmieren habe ich mir mit Büchern und dem Internet selbst bei gebracht.

Der Turm mit Drehpunkt ist fertig.
Foto wollte ich anhängen, geht aber nicht wie erwartet.
Für das Schwenken liegt ein Linearantrieb im Keller.
Den Schrittmotor habe ich Gersten auch endlich in den Griff bekommen.
Also alle Stellglieder sind fertig.
Bei der Planung bin ich auf Euer Forum gestoßen.
Im obigen Verlaufsteht ein Programm, das ich in Teilen in mein Programm einfügen wollte.
Berechnung der beiden Winkel, die ich zum Nachführen benötige.
Also kopiert und mit Anfang und Ende versehen.
Es kamen auch Werte heraus, aber wenn ich die mit der Seite: Sonnenverlauf.de vergleiche, passt es einfach nicht.
Sonnenanbeter hatte geschrieben, dass er zwei Fehler gefunden hatte und die Berechnung jetzt funktioniert würde.
Ob es an mir liegt oder an den Programmzeilen????

Jahr, monat........Sekunden konnte ich prüfen, das schein zu passen.
Aber bei den anderen Werten benötige ich die Berechnungsgrundlage die angewendet wurde.

Wühle mich gerade bei Wiki durch die Seite Sonnenstand.
Ähnlich, aber nicht gleich.

Gruß
Lobi
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

Vielleicht ist auch ein Blick auf `astropy` sinnvoll, das ist ein relativ großes und weitverbreitetes Paket zu astronomischen Berechnungen. Für die Sonnenposition haben die sogar schon ein fertiges Beispiel (ganz unten).
sonnenanbeter
User
Beiträge: 7
Registriert: Samstag 9. Januar 2016, 11:47

Hallo zusammen,
entschuldigt die verspätete Antwort, aber Python ist gerade nicht im Fokus. Da ich zwar die Steuerung und einen kleinen Testaufbau, aber keine große Gartenanlage realilsieren konnte, habe ich das Projekt vollständig eingestellt. Es gibt mittlerweile viele andere und einfacher gelöste vergleichbare Projekte.
Viel Spaß und viel ERfolg!
Der Sonnenanbeter
Antworten