Mathematik zur Berechnung der terrestrischen Länge direkt unter der Sonne mit der Zeit

Ich versuche, den Längengrad auf der Erde zu berechnen, wo es irgendwann Mittag ist. (Das heißt, der Längengrad, der koplanar mit der Ebene ist, die durch die Sonne und die Erdachse definiert ist.)

Hier ist mein Python-Code:

from math import sin, cos, tan, atan, pi
def sun_longitude(when):
    """Given time in ms since 1/1/1970, return the longitude the sun is over at that moment"""
    # https://en.wikipedia.org/wiki/Position_of_the_Sun
    jdn = 2440587.5 + when / (1000.0 * 3600 * 24)
    n = jdn - 2451545.0 # 1/1/2000
    L = (280.460 + 0.9856474 * n) % 360.0
    g = (357.528 + 0.9856004 * n) % 360.0
    degtorad = 2.0 * pi / 360.0
    lambda_ = L + 1.915 * sin(g * degtorad) + 0.020 * sin(2 * g * degtorad)

    # https://en.wikipedia.org/wiki/Axial_tilt#Short_term
    T = n / (365 * 100.0) # julian centuries from 2000
    # ε = 23° 26′ 21.406″ − 46.836769″ T − 0.0001831″ T2
    epsilon = 23.4392794 - 0.780612817 * T - 5.0861e-8 * T * T;
    alpha = atan(cos(epsilon * degtorad) * tan(lambda_ * degtorad)) / degtorad

    # https://en.wikipedia.org/wiki/Sidereal_time
    # https://en.wikipedia.org/wiki/Hour_angle
    GMST = (18.697374558 + 24.06570982441908 * n) % 24.0
    print "jdn = {jdn}, n = {n}, L = {L}, g = {g}, lambda_ = {lambda_}, T = {T}, epsilon = {epsilon}, alpha = {alpha}, GMST = {GMST}".format(**locals())
    LHA = GMST * 360/24 + alpha;

    return LHA

from datetime import datetime
def sun_long_from_str(whenstr):
    when = datetime.strptime(whenstr, "%Y-%m-%d %H:%M")
    secs = (when - datetime(1970, 1, 1)).total_seconds()
    return sun_longitude(secs * 1000.0)

Mir ist bewusst, dass ich bei der Berechnung von Alpha nicht den richtigen Quadranten korrigiere, aber ich habe an anderer Stelle Fehler; Zum Beispiel erhalte ich für Mittag GMT am 1. Januar 2000 279 ° und nicht nahe 0, was ich erwarten würde.

Dies gibt mir keine richtigen Antworten, und ich bin ein bisschen ratlos, wie ich es debuggen soll. Kann jemand meine Fehler finden oder mich auf einen vernünftigen Beispielcode oder ein durchgearbeitetes Beispiel dafür hinweisen?

Denn sobald ich den Algorithmus richtig verstehe, werde ich für ein eingebettetes Gerät neu implementieren, ich kann nicht einfach ein Implementierungspaket verwenden, und ich habe keine Bibliothek gefunden, die einfach genug ist, um zu verstehen, wie ich dies implementieren kann konkrete Frage.

Danke für ein paar Fehlerkorrekturen. Wenn ich es jetzt ausführe, erhalte ich für diese Beispiele Folgendes:

2000-01-01 12:00:
  jdn = 2451545.0, n = 0.0, L = 280.46, g = 357.528, lambda_ = 280.375680197, T = 0.0, epsilon = 23.4392794, alpha = -78.7141369122, GMST = 18.697374558
  result: 201.74648145775257

2000-03-20 07:35:
  jdn = 2451623.81597, n = 78.815972222, L = 358.144758099, g = 75.2090537484, lambda_ = 360.006175505, T = 0.00215934170471, epsilon = 23.4375937902, alpha = 0.00566598789588, GMST = 19.4596915825
  291.9010397253072

Wie Sie sehen können, hat die erste Abfrage (mittags am 1. Januar 2000) jetzt eine korrekte julianische Tageszahl. Um Mittag GMT würde ich erwarten, dass die Sonne über einem Längengrad nahe 0 steht, also ist 201 falsch.

Die zweite Abfrage ist die Zeit der Frühlings-Tagundnachtgleiche im Jahr 2000, von der ich erwartet hatte, dass sie dazu führen würde, dass der siderische Teil der Berechnung auf Null gesetzt wird, aber das ist nicht der Fall.

Ich habe am 1.1.2000 mittags GMT versucht, die NASA HORIZONS-Webschnittstelle nach Ephemeriden für die Sonne zu konsultieren, sowohl für "geozentrische" als auch für Greenwich-Standorte, aber ich weiß nicht, wie ich das Ergebnis auswerten und mit meiner Arbeit vergleichen soll Oben.

Geozentrisch:

 Date__(UT)__HR:MN     R.A._(ICRF/J2000.0)_DEC  APmag  S-brt            delta      deldot    S-O-T /r    S-T-O
**************************************************************************************************************
 2000-Jan-01 12:00     18 45 09.36 -23 01 59.7 -26.78 -10.59 0.98332762653520  -0.0127281   0.0000 /?   0.0000

Greenwich:

 Date__(UT)__HR:MN     R.A._(ICRF/J2000.0)_DEC  APmag  S-brt            delta      deldot    S-O-T /r    S-T-O
**************************************************************************************************************
 2000-Jan-01 12:00 *m  18 45 09.36 -23 02 08.3 -26.78 -10.59 0.98331613178086  -0.0166655   0.0000 /?   0.0000

Nochmals vielen Dank für jede Hilfe, die Sie leisten können.

[Bearbeitet, um die Berechnung des Julianischen Tages und die Verwendung von Grad vs. Bogenmaß für Alpha zu korrigieren und Beispiele hinzuzufügen]

Allgemeiner Debugging-Vorschlag: Drucken Sie die Variablen aus, nachdem Sie sie berechnet haben, um zu sehen, wo Sie falsch liegen? Ich denke auch, dass Sie dies tun können, indem Sie einfach die Zeitgleichung und die Tatsache verwenden, dass sich die Erde alle 24 Stunden dreht (in Bezug auf die Sonne).
Könntest du uns auch die Ausgabe zeigen?
Als erstes ist mir aufgefallen, dass die Unix-Zeit in Sekunden ab Mitternacht gemessen wird, während Julianische Daten ab Mittag gemessen werden. Versuchenjdn = 2440587.5 + when / (1000.0 * 3600 * 24)
@barrycarter es hängt davon ab, wie genau Sie sein wollen. Da sich die Erde in einer elliptischen Umlaufbahn (ändernde Geschwindigkeit) befindet, variiert der wörtliche Sonnentag leicht von 24 Stunden, je nachdem, wo Sie sich in der Umlaufbahn befinden. Der Sonnentag dauert im Durchschnitt genau 24 Stunden.
Das nächste, was seltsam erscheint, ist, dass "Alpha" in Radiant berechnet wird, anstatt in den üblichen Stunden Min Sec. und Lambda ist in Grad.
@Aabaakawad Ja, deshalb habe ich die Zeitgleichung erwähnt (als Korrektur des 24-Stunden-Rotationszyklus mit / gegen die Sonne)
Danke für alle Kommentare. @barrycarter -- Es druckt die Werte kurz vor dem Ende.
@james-kifiger: danke; Ich habe diese Korrektur für die Berechnung der Sternzeit vorgenommen (siehe n - 0,5), aber aus irgendeinem Grund hatte ich den Eindruck, dass die für die anderen Berechnungen verwendeten julianischen Tage um Mitternacht begannen. Und natürlich hast du Recht mit Alpha!

Antworten (1)

Es gibt zwei Dinge, die ich in diesem Code sehen kann:

Erstens wird das julianische Datum ab Mittag gemessen, während die Unix-Epoche ab Mitternacht gemessen wird. jdn = 2440587.5 + when / (1000.0 * 3600 * 24)sollte der richtige Ausdruck sein.

Zweitens, alpha = atan(cos(epsilon * degtorad) * tan(lambda_ * degtorad))berechnet eine Rektaszension im Bogenmaß, sollten Sie diese für den letzten Teil der Berechnungen in Grad umrechnen.

Bearbeiten Sie müssen den Quadranten korrigieren, zum Beispiel wenn n == 0 ist, sollten Sie eine rechte Aszension von 281,285 erhalten. Der letzte Fehler liegt jedoch in der Berechnung des Längengrades.

Die relevante Gleichung ist L H EIN = G M S T − λ − a . Wenn die Sonne auf dem Meridian steht, ist LHA = 0 und (Umordnung) λ = G M S T − a . So LHA = GMST * 360/24 + alphasollte die Zeile lauten longitude = GMST*360/24 - alpha. (und return longitude) Wenn Sie dies mit n=0 tun, erhalten Sie einen Längengrad von -1,5 Grad. Die Sonne steht am Mittag des 1. Januar 2000 auf dem Meridian, wenn Sie sich 1,5 Grad westlich von Greenwich befinden.

Sie scheinen eine sehr genaue Ephemeride für die Sonne zu versuchen. Sie können Ihre Genauigkeit anhand von PyEphem oder den Ephemeriden von Nasa Horizon überprüfen .

Danke, ich habe diese Korrekturen vorgenommen, sehr hilfreich. Die Ergebnisse scheinen jedoch immer noch falsch zu sein.
Ich denke, der Fehler liegt in der Berechnung von LHA, die LHA = GMST - Alpha sein sollte. Auch der Wert von Alpha sollte im Bereich von 0 bis 360 liegen, während Atan im Bereich von -90 bis 90 zurückgibt, Sie müssen über die Korrektur des Quadranten nachdenken
Danke, die werde ich mir anschauen. Sollte ich nicht erwarten, dass die GMST zur Frühlings-Tagundnachtgleiche Null ist?
Mathematik zur Berechnung der terrestrischen Länge direkt unter der Sonne mit der Zeit
AntwortenHierDatenschutz-Bestimmungen1y, 0v