代码之家  ›  专栏  ›  技术社区  ›  dassouki

日出/日落计算

  •  6
  • dassouki  · 技术社区  · 15 年前

    我在用下面的链接计算日落时间。

    我通过excel和python得到的结果与实际值不匹配。你知道我做错了什么吗?

    我的Excel表格可以在。。 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()
    
    3 回复  |  直到 15 年前
        1
  •  3
  •   MattH    15 年前

    为什么所有的电话 radians degrees ? 我以为输入的数据已经是十进制的了。

    如果我:

    • 去掉所有的弧度和度数
    • 45.9 -66.6 分别地
    • 将utc中的时间\u更正为0到24之间。

    编辑:
    作为J。F塞巴斯蒂安指出,根据问题中链接的电子表格,这个位置的日出时间的答案,以及使用观察者类提供的答案 ephem 在07:01-07:02之间。

    一旦我在正确的范围内得到了一个数字(实现中的注释是07:34),我就不再寻找dassouki实现美国海军天文台算法的错误了。

    here . 然而,根据我最近对此事的了解,这些变化只会导致日出时间相差几分钟,而不是超过半小时。

        2
  •  10
  •   jfs    15 年前

    你可以用 ephem 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')
    

    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)
    

    输出:

    Beijing 06:02
    Berlin 19:45
    Moscow 06:53
    Fredericton 07:01
    
        3
  •  1
  •   rlotun    15 年前

    我怀疑这与没有实际执行浮点除法有关。在python中,如果a和b都是整数,那么a/b也是整数:

       $ python
       >>> 1 / 2
       0
    

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