日出/布景计算 [英] Sunrise / set calculations
问题描述
我正在尝试根据下面提供的链接使用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()
推荐答案
为什么所有对radians
和degrees
的调用?我以为输入数据已经是十进制数了.
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屋!