日出/布景计算 [英] Sunrise / set calculations

查看:60
本文介绍了日出/布景计算的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我正在尝试根据下面提供的链接使用python计算日落/上升时间.

I'm trying to calculate the sunset / rise times using python based on the link provided below.

我通过excel和python完成的结果与实际值不匹配.对我可能做错的事情有任何想法吗?

My results done through excel and python do not match the real values. Any ideas on what I could be doing wrong?

我的Excel工作表位于.. http://transpotools.com/sun_time.xls

My Excel sheet can be found under .. 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()

推荐答案

为什么所有对radiansdegrees的调用?我以为输入数据已经是十进制数了.

Why all the calls to radians and degrees? I thought the input data was already in decimal degrees.

如果我:7:37 am,我得到的结果是

I get a result of 7:37am if I:

  • 消除所有对弧度和度的呼叫
  • 将纬度/经度分别校正为:45.9-66.6
  • 将time_in_utc正确设置为0到24之间.
  • strip out all of the calls to radians and degrees
  • correct the lat / long to: 45.9 and -66.6 respectively
  • correct time_in_utc to fall within 0 and 24.

修改:
正如JF Sebastian指出的那样,根据问题中链接的电子表格,使用ephem的Observer类提供的答案,此位置处的日出时间的答案在07:01-07:02范围内.


As J. F. Sebastian points out, the answer for the sunrise time at this location according to the spreadsheet linked in the question and the answer provided by using the Observer class of the ephem are in the region of 07:01-07:02.

一旦我在正确的地点得到了一个数字(在实施注释中的07:34),我就停止寻找dassouki实施美国海军天文台算法的错误.

I stopped looking for errors in dassouki's implementation of the US Naval Observatory's algorithm once I got a figure in the right ballpark (07:34 in the comments in the implementation).

对此算法进行了简化,对构成日出"的内容也有所不同,其中的一些内容在此处.但是,从我最近在此问题上获得的经验来看,这些变化只会导致日出时间差几分钟,而不是半个小时.

Looking into it, this algorithm makes some simplifications and there is variation about what constitutes 'sunrise', some of this is discussed here. However, in my opinion from what I've recently learnt on this matter, these variations should only lead to a difference of a few minutes in sunrise time, rather than over half an hour.

这篇关于日出/布景计算的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆