Calcoli alba/set

Question

Sto provando a calcolare i tempi di tramonto/aumento utilizzando python in base al collegamento fornito di seguito.

I miei risultati eseguiti tramite Excel e Python non corrispondono ai valori reali. Qualche idea su cosa potrei fare male?

foglio My Excel si trova sotto .. 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

Perché tutte le chiamate a radians e degrees? Pensavo che i dati di input fossero già in gradi decimali.

ottengo un risultato di 7:37 se:

• striscia fuori tutte le chiamate a radianti e gradi
• correggere il lat/lungo per: 45.9 e -66.6 rispettivamente
• corretta time_in_utc rientrare 0 e 24.

---

Modifica:
Come fa notare JF Sebastian, la risposta per il tempo di alba in questa posizione in base al foglio di calcolo collegato nella domanda e la risposta fornita utilizzando la classe Observer dello ephem sono nell'area 07: 01- 07:02.

Ho smesso di cercare errori nell'implementazione dell'algoritmo del US Naval Observer da parte di Dassouki, una volta che ho ottenuto una figura nel campo di destra (07:34 nei commenti nell'implementazione).

Guardando in esso, questo algoritmo rende alcune semplificazioni e c'è variazione su ciò che costituisce 'alba', parte di questo è discusso here. Tuttavia, a mio avviso, da quanto ho appreso di recente su questo argomento, queste variazioni dovrebbero solo portare ad una differenza di pochi minuti in ora di alba, piuttosto che oltre mezz'ora.

Answer 2

ho il sospetto che questo ha qualcosa a che fare con realtà non eseguire la divisione in virgola mobile. In Python se a e b sono entrambi interi, a/b è anche un numero intero:

$ python 
    >>> 1/2 
    0 

Le opzioni sono o di costringere a stare a galla uno dei tuoi argomenti (cioè, invece di a/b fare un galleggiante (a)/b) o per assicurarsi che il '/' si comporta correttamente in un Python 3K modo:

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

Quindi, se si tiene fede che la dichiarazione di importazione nella parte superiore del file, si può risolvere il problema. Ora/produrrò sempre un float e, per ottenere il vecchio comportamento, puoi usare // invece.

Answer 3

Si potrebbe utilizzare ephem modulo python:

#!/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') 

Esempio:

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) 

uscita:

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