Cálculos de amanecer/conjunto

Question

Estoy tratando de calcular los tiempos de puesta de sol/subida usando Python basado en el enlace proporcionado a continuación.

Mis resultados realizados a través de Excel y Python no coinciden con los valores reales. ¿Alguna idea sobre lo que podría estar haciendo mal?

hoja de Excel Mi se encuentra bajo .. 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

¿Por qué todas las llamadas a radians y degrees? Pensé que los datos de entrada ya estaban en grados decimales.

puedo obtener un resultado de 07:37 am si:

• tira a cabo todas las llamadas a radianes y grados
• corregir la latitud/longitud a: 45.9 y -66.6 respectivamente
• time_in_utc correcta a caer dentro de 0 y 24.

---

Edit:
Como señala JF Sebastian, la respuesta para la hora del amanecer en esta ubicación de acuerdo con la hoja de cálculo vinculada en la pregunta y la respuesta proporcionada mediante el uso de la clase Observer del ephem están en la región de 07: 01- 07:02.

Dejé de buscar errores en la implementación de dassouki del algoritmo del Observatorio Naval de los EE. UU. Una vez que obtuve una cifra en el estadio correcto (07:34 en los comentarios en la implementación).

Al observarlo, este algoritmo realiza algunas simplificaciones y hay una variación acerca de lo que constituye 'amanecer', parte de esto se discute here. Sin embargo, en mi opinión, por lo que aprendí recientemente sobre este tema, estas variaciones solo deberían conducir a una diferencia de unos pocos minutos en el tiempo del amanecer, en lugar de una media hora.

Answer 2

Sospecho que esto tiene algo que ver con no realizar realmente división de coma flotante. En Python si a y b son dos números enteros, a/b es también un entero:

$ python 
    >>> 1/2 
    0 

Sus opciones son o coaccionar a flotar uno de sus argumentos (es decir, en lugar de a/b hacer un flotador (a)/b) o para asegurarse de que la '/' se comporta adecuadamente en un pitón 3K manera:

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

Así que si usted se pega que la declaración de importación en la parte superior de su archivo, puede arreglar el problema. Ahora/siempre producirá un flotador, y para obtener el comportamiento anterior puede usar // en su lugar.

Answer 3

Usted podría utilizar ephem módulo de 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') 

Ejemplo:

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) 

Salida:

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