title: Sunrise / set calculations
body: <p>I'm trying to calculate the sunset / rise times using python based on the link provided below.</p>
  
  <p>My results done through excel and python do not match the real values. Any ideas on what I could be doing wrong?</p>
  
  <p>My Excel sheet can be found under .. <a href="http://transpotools.com/sun_time.xls" rel="noreferrer">http://transpotools.com/sun_time.xls</a></p>
  
  <pre><code># 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 &lt; 0:
        true_long = true_long + 360
      elif true_long &gt; 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 &lt; 0:
        s_r_a = s_r_a + 360
      elif true_long &gt; 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 &gt; 1:
        print "Sun Never Rises at this location on this date, exiting"
        # sys.exit()
      elif cos_hour &lt; -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()
  </code></pre>