日出/日落计算
我正在尝试使用Python根据下面提供的链接来计算日出和日落的时间。
我通过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 个回答
1
我怀疑这和没有进行浮点数除法有关。在Python中，如果a和b都是整数，那么a / b的结果也是整数：
$ python
>>> 1 / 2
0
你可以选择将其中一个参数转换为浮点数（也就是说，不是用a/b，而是用float(a) / b），或者确保'/'在Python 3中正常工作：
$ python
>>> from __future__ import division
>>> 1 / 2
0.5
所以如果你在文件的开头加上那个导入语句，可能就能解决你的问题。这样一来，'/'总是会产生浮点数，如果你想要以前的行为，可以用'//'来代替。
10
你可以使用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
为什么要调用 radians 和 degrees 呢？我以为输入的数据已经是十进制度数了。
如果我：
- 去掉所有对 radians 和 degrees 的调用
- 把经纬度改成：
45.9和-66.6 - 把 time_in_utc 调整到 0 到 24 之间。
编辑：
正如 J. F. Sebastian 指出的，根据问题中链接的电子表格和使用 ephem 的 Observer 类得到的答案，这个地点的日出时间大约在 07:01 到 07:02 之间。
当我在评论中看到一个大致正确的结果（07:34）时，我就不再寻找 dassouki 实现的美国海军天文台算法中的错误了。
深入研究后发现，这个算法做了一些简化，而且关于什么算是“日出”也有不同的看法，部分内容可以在 这里 讨论。不过，按照我最近学到的，这些差异应该只会导致日出时间相差几分钟，而不是超过半小时。