Sunrise / Set of calculations

Question

I am trying to calculate the sunset / rise time using python based on the link below.

My results done with excel and python do not match the actual values. Any ideas on what I might be doing wrong?

My Excel sheet can be found at 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 true it indicates rise, # and if it 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 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() 

Answer 1

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

I get a result of 7:37 if I:

• delete all calls to radians and degrees
• fix lat / long for: 45.9 and -66.6 respectively
• The correct time_in_utc value must be between 0 and 24.

---

Edit:
As JF Sebastian points out, the answer to the sunrise time in this place according to the spreadsheet related to the question and the answer obtained using the Observer ephem class are in the area 07: 01-07: 02.

I stopped looking for errors in dassouki using the US Observatory algorithm when I got the number in the right step (07:34 in the implementation comments).

Looking at this, this algorithm makes some simplifications, and there is a variation as to what constitutes a “sunrise”, some of which are discussed here . However, in my opinion, from what I recently learned on this issue, these variations should lead to a difference of several minutes during sunrise, and not half an hour later.

Answer 2

You can use the ephem python module:

 #!/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') 

Example:

 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) 

Output:

 Beijing 06:02 Berlin 19:45 Moscow 06:53 Fredericton 07:01 

Answer 3

I suspect this is due to the fact that floating point separation is not actually performed. In python, if a and b are integers, a / b is also an integer:

  $ python >>> 1 / 2 0 

Your options are to force one of your arguments (i.e., instead of a / b, make float (a) / b) or make sure that "/" behaves correctly in the Python 3K path:

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

So, if you stick with this import statement at the top of the file, it may solve your problem. Now / will always create a float, and to get the old behavior you can use // instead.