Sonnenaufgang/Set Berechnungen

Question

Ich versuche, die Sonnenuntergang/Anstieg Zeiten mit Python basierend auf dem unten angegebenen Link zu berechnen.

Meine Ergebnisse über Excel und Python stimmen nicht mit den tatsächlichen Werten überein. Irgendwelche Ideen, was ich falsch machen könnte?

Mein Excel-Tabelle finden Sie unter .. http://transpotools.com/sun_time.xls

# Created on 2010-03-28 

# @author: dassouki 
# @source: [http://williams.best.vwh.net/sunrise_sunset_algorithm.htm][2] 
# @summary: this is based on the Nautical Almanac Office, United States Naval 
# Observatory. 

import math, sys 

class TimeOfDay(object): 

    def calculate_time(self, in_day, in_month, in_year, 
        lat, long, is_rise, utc_time_zone): 

    # is_rise is a bool when it's true it indicates rise, 
    # and if it's false it indicates setting time 

    #set Zenith 
    zenith = 96 
    # offical  = 90 degrees 50' 
    # civil  = 96 degrees 
    # nautical  = 102 degrees 
    # astronomical = 108 degrees 


    #1- calculate the day of year 
    n1 = math.floor(275 * in_month/9) 
    n2 = math.floor((in_month + 9)/12) 
    n3 = (1 + math.floor(in_year - 4 * math.floor(in_year/4) + 2)/3) 

    new_day = n1 - (n2 * n3) + in_day - 30 

    print "new_day ", new_day 

    #2- calculate rising/setting time 
    if is_rise: 
     rise_or_set_time = new_day + ((6 - (long/15))/24) 
    else: 
     rise_or_set_time = new_day + ((18 - (long/ 15))/24) 

    print "rise/set", rise_or_set_time 

    #3- calculate sun mean anamoly 
    sun_mean_anomaly = (0.9856 * rise_or_set_time) - 3.289 
    print "sun mean anomaly", sun_mean_anomaly 

    #4 calculate true longitude 
    true_long = (sun_mean_anomaly + 
        (1.916 * math.sin(math.radians(sun_mean_anomaly))) + 
        (0.020 * math.sin( 2 * math.radians(sun_mean_anomaly))) + 
        282.634) 
    print "true long ", true_long 

    # make sure true_long is within 0, 360 
    if true_long < 0: 
     true_long = true_long + 360 
    elif true_long > 360: 
     true_long = true_long - 360 
    else: 
     true_long 

    print "true long (360 if) ", true_long 

    #5 calculate s_r_a (sun_right_ascenstion) 
    s_r_a = math.degrees(math.atan(0.91764 * math.tan(math.radians(true_long)))) 
    print "s_r_a is ", s_r_a 

    #make sure it's between 0 and 360 
    if s_r_a < 0: 
     s_r_a = s_r_a + 360 
    elif true_long > 360: 
     s_r_a = s_r_a - 360 
    else: 
     s_r_a 

    print "s_r_a (modified) is ", s_r_a 

    # s_r_a has to be in the same Quadrant as true_long 
    true_long_quad = (math.floor(true_long/90)) * 90 
    s_r_a_quad = (math.floor(s_r_a/90)) * 90 
    s_r_a = s_r_a + (true_long_quad - s_r_a_quad) 

    print "s_r_a (quadrant) is ", s_r_a 

    # convert s_r_a to hours 
    s_r_a = s_r_a/15 

    print "s_r_a (to hours) is ", s_r_a 


    #6- calculate sun diclanation in terms of cos and sin 
    sin_declanation = 0.39782 * math.sin(math.radians (true_long)) 
    cos_declanation = math.cos(math.asin(sin_declanation)) 

    print " sin/cos declanations ", sin_declanation, ", ", cos_declanation 

    # sun local hour 
    cos_hour = (math.cos(math.radians(zenith)) - 
       (sin_declanation * math.sin(math.radians (lat)))/
       (cos_declanation * math.cos(math.radians (lat)))) 

    print "cos_hour ", cos_hour 

    # extreme north/south 
    if cos_hour > 1: 
     print "Sun Never Rises at this location on this date, exiting" 
     # sys.exit() 
    elif cos_hour < -1: 
     print "Sun Never Sets at this location on this date, exiting" 
     # sys.exit() 

    print "cos_hour (2)", cos_hour 


    #7- sun/set local time calculations 
    if is_rise: 
     sun_local_hour = (360 - math.degrees(math.acos(cos_hour)))/15 
    else: 
     sun_local_hour = math.degrees(math.acos(cos_hour))/15 


    print "sun local hour ", sun_local_hour 

    sun_event_time = sun_local_hour + s_r_a - (0.06571 * 
               rise_or_set_time) - 6.622 

    print "sun event time ", sun_event_time 

    #final result 
    time_in_utc = sun_event_time - (long/15) + utc_time_zone 

    return time_in_utc 



#test through main 
def main(): 
    print "Time of day App " 
    # test: fredericton, NB 
    # answer: 7:34 am 
    long = 66.6 
    lat = -45.9 
    utc_time = -4 
    d = 3 
    m = 3 
    y = 2010 
    is_rise = True 

    tod = TimeOfDay() 
    print "TOD is ", tod.calculate_time(d, m, y, lat, long, is_rise, utc_time) 


if __name__ == "__main__": 
    main() 

Answer 1

Warum alle Anrufe zu radians und degrees? Ich dachte, die Eingabedaten wären bereits in Dezimalgraden.

bekomme ich ein Ergebnis von 7.37, wenn ich: 45.9 und -66.6 bzw.

korrekt time_in_utc:

• Streifen aus alle Anrufe in Radiant und Grad
• die lat/long zu korrigieren fallen zwischen 0 und 24.

---

Edit:
Wie JF Sebastian darauf hinweist, liegen die Antworten für die Sonnenaufgangszeit an diesem Standort gemäß der in der Frage verlinkten Tabelle und der Antwort der Observer Klasse ephem in der Region 07: 01- 07:02.

Ich habe aufgehört, nach Fehlern in Dassoukis Implementierung des Algorithmus des US Naval Observatory zu suchen, sobald ich eine Zahl im rechten Baseballstadion bekommen habe (07:34 in den Kommentaren in der Implementierung).

Ein Blick darauf, dieser Algorithmus macht einige Vereinfachungen und es gibt Variationen darüber, was "Sonnenaufgang" ausmacht, einige davon werden diskutiert here. Meiner Meinung nach, was ich kürzlich in dieser Sache gelernt habe, sollten diese Variationen jedoch nur zu einem Unterschied von ein paar Minuten in der Sonnenaufgangszeit und nicht über eine halbe Stunde führen.

Answer 2

Ich vermute, das etwas mit nicht wirklich durchführen Gleitkomma-Division zu tun hat. In Python, wenn a und b beide ganze Zahlen sind, a/b ist auch eine ganze Zahl:

$ python 
    >>> 1/2 
    0 

Ihre Optionen sind entweder zu zwingen, eines Ihrer Argumente zu schweben (das heißt, anstelle von a/b machen einen Schwimmer (a)/b) oder um sicherzustellen, dass die ‚/‘ richtig in einem Python-3K verhält:

$ python 
    >>> from __future__ import division 
    >>> 1/2 
    0.5 

Also, wenn Sie diese Import-Anweisung am Anfang der Datei bleiben, ist es Ihr Problem beheben kann. Jetzt wird/immer eine Gleitkommazahl erzeugen, und um das alte Verhalten zu erhalten, kann man stattdessen // verwenden.

Answer 3

könnten Sie ephem Python-Modul verwenden:

#!/usr/bin/env python 
import datetime 
import ephem # to install, type$ pip install pyephem 

def calculate_time(d, m, y, lat, long, is_rise, utc_time): 
    o = ephem.Observer() 
    o.lat, o.long, o.date = lat, long, datetime.date(y, m, d) 
    sun = ephem.Sun(o) 
    next_event = o.next_rising if is_rise else o.next_setting 
    return ephem.Date(next_event(sun, start=o.date) + utc_time*ephem.hour 
        ).datetime().strftime('%H:%M') 

Beispiel:

for town, kwarg in { "Fredericton": dict(d=3, m=3, y=2010, 
             lat='45.959045', long='-66.640509', 
             is_rise=True, 
             utc_time=20), 

        "Beijing": dict(d=29, m=3, y=2010, 
            lat='39:55', long='116:23', 
            is_rise=True, 
            utc_time=+8), 

        "Berlin": dict(d=4, m=4, y=2010, 
            lat='52:30:2', long='13:23:56', 
            is_rise=False, 
            utc_time=+2) , 

        "Moscow": dict(d=4, m=4, y=2010, 
            lat='55.753975', long='37.625427', 
            is_rise=True, 
            utc_time=4) }.items(): 
    print town, calculate_time(**kwarg) 

Ausgang:

Beijing 06:02 
Berlin 19:45 
Moscow 06:53 
Fredericton 07:01