Mặt trời mọc/đặt tính toán

Question

Tôi đang cố gắng tính thời gian hoàng hôn/tăng lên bằng cách sử dụng python dựa trên liên kết được cung cấp bên dưới.

Kết quả của tôi được thực hiện thông qua excel và python không khớp với giá trị thực. Bất kỳ ý tưởng về những gì tôi có thể làm sai?

My Excel tờ có thể được tìm thấy dưới .. 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

Tại sao tất cả các cuộc gọi đến radians và degrees? Tôi nghĩ dữ liệu đầu vào đã ở độ thập phân.

tôi nhận được một kết quả của 07:37 nếu tôi:

• dải ra tất cả các cuộc gọi đến radian và độ
• sửa lat/dài để: 45.9 và -66.6 tương ứng
• time_in_utc đúng rơi trong 0 và 24.

---

Chỉnh sửa:
Như JF Sebastian đã chỉ ra, câu trả lời cho thời gian mặt trời mọc tại vị trí này theo bảng tính được liên kết trong câu hỏi và câu trả lời được cung cấp bằng cách sử dụng lớp Observer của ephem nằm trong vùng 07: 01- 07:02.

Tôi đã ngừng tìm lỗi trong việc triển khai thuật toán của Đài quan sát Hải quân Hoa Kỳ của dassouki khi tôi đã có một hình bóng ở sân chơi bóng bên phải (07:34 trong các ý kiến ​​trong quá trình thực hiện).

Nhìn vào nó, thuật toán này làm cho một số đơn giản hóa và có sự thay đổi về những gì tạo thành 'mặt trời mọc', một số điều này được thảo luận here. Tuy nhiên, theo ý kiến ​​của tôi từ những gì tôi đã học gần đây về vấn đề này, các biến thể này sẽ chỉ dẫn đến một sự khác biệt của một vài phút trong thời gian mặt trời mọc, thay vì hơn nửa giờ.

Answer 2

Tôi nghi ngờ điều này có cái gì để làm với không thực sự thực hiện phân dấu chấm động. Trong python nếu a và b đều số nguyên, a/b cũng là một số nguyên:

$ python 
    >>> 1/2 
    0 

lựa chọn của bạn là một trong hai để ép buộc nổi một trong những lập luận của bạn (có nghĩa là, thay vì a/b làm một phao (a)/b) hoặc để đảm bảo '/' cư xử đúng cách trong một Python 3K cách:

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

Vì vậy, nếu bạn dính rằng tuyên bố nhập khẩu ở phía trên cùng của tập tin của bạn, nó có thể khắc phục vấn đề của bạn. Bây giờ/sẽ luôn tạo ra một phao, và để có được hành vi cũ, bạn có thể sử dụng // thay vào đó.

Answer 3

Bạn có thể sử dụng ephem python module:

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

Ví dụ:

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) 

Output:

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