I'm trying to calculate the sunset / rise times using python based on the link provided below.
My results done through excel and python do not match the real values. Any ideas on what I could be doing wrong?
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()
Why all the calls to radians and degrees? I thought the input data was already in decimal degrees.
I get a result of 7:37am if I:
45.9 and -66.6 respectivelyEdit:
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.
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.