attachment:sunrisesunset.py of SunriseSunset

   1 #!/usr/bin/env python
   2 #
   3 # sunrisesunset.py
   4 #
   5 # This code is valid for dates from 1901 to 2099, and will not calculate
   6 # sunrise/sunset times for latitudes above/below 63/-63 degrees.
   7 #
   8 # No external packages are used when using the class SunriseSunset. If you
   9 # run the tests you will need to install pytz as shown below or use your
  10 # package installer if it's available.
  11 #
  12 # $ sudo easy_install --upgrade pytz
  13 #
  14 # CVS/SVN Info
  15 #----------------------------------
  16 # $Author: cnobile $
  17 # $Date: 2009-08-03 00:47:13 $
  18 # $Revision: 1.9 $
  19 #----------------------------------
  20 #
  21 # Contributions by:
  22 #    Darryl Smith -- Noticed a bug with atan fixed with atan2.
  23 #
  24 ##########################################################################
  25 # Copyright (c) 2009 Carl J. Nobile.
  26 # All rights reserved. This program and the accompanying materials
  27 # are made available under the terms of the Eclipse Public License v1.0
  28 # which accompanies this distribution, and is available at
  29 # http://www.eclipse.org/legal/epl-v10.html
  30 #
  31 # This Copyright covers only the implimentation of the algorithm used in
  32 # this code and does not cover the algorithm itself which is in the
  33 # public domain.
  34 #
  35 # Contributors:
  36 #    Carl J. Nobile - initial API and implementation
  37 ##########################################################################
  38 
  39 import datetime
  40 from math import degrees, radians, atan2, cos, sin, pi, sqrt, fabs
  41 
  42 
  43 class SunriseSunset(object):
  44     """
  45     This class determines the sunrise and sunset for zenith standards for the
  46     given day. It can also tell you if the given time is during the nigh or
  47     day.
  48     """
  49     __ZENITH = {'official': -0.833,
  50                 'civil': -6.0,
  51                 'nautical': -12.0,
  52                 'amateur': -15.0,
  53                 'astronomical': -18.0}
  54 
  55     def __init__(self, date, lat, lon, zenith='official'):
  56         """
  57         Set the values for the sunrise and sunset calculation.
  58 
  59         @param date: A localized datetime object that is timezone aware.
  60         @param lat: The latitude.
  61         @param lon: The longitude.
  62         @keyword zenith: The zenith name.
  63         """
  64         if not isinstance(date, datetime.datetime) or date.tzinfo is None:
  65             msg = "The date must be a datetime object with timezone info."
  66             raise ValueError(msg)
  67 
  68         if zenith not in self.__ZENITH:
  69             msg = "Invalid zenith name [%s] must be one of: %s"
  70             raise ValueError(msg % (zenith, self.__ZENITH.keys()))
  71 
  72         if abs(lat) > 63:
  73             raise ValueError('Invalid latitude: %s' % lat)
  74 
  75         self.__dateLocal = date
  76         self.__lat = lat
  77         self.__lon = lon
  78         self.__zenith = zenith
  79         localTuple = date.timetuple()
  80         utcTuple = date.utctimetuple()
  81         self.__offsetUTC = (localTuple[3] - utcTuple[3]) + \
  82                            (localTuple[4] - utcTuple[4]) / 60.0
  83         self.__sunrise = None
  84         self.__sunset = None
  85         self.__determineRiseSet()
  86 
  87     def isNight(self, collar=0):
  88         """
  89         Check if it is day or night. If the 'collar' keyword argument is
  90         changed it will skew the results to either before or after the real
  91         sunrise and sunset. This is useful if lead and lag timea are needed
  92         around the actual sunrise and sunset.
  93 
  94         Note::
  95             If collar == 30 then this method will say it is daytime 30
  96             minutes before the actual sunrise and likewise 30 minutes after
  97             sunset it would indicate it is night.
  98 
  99         @keyword collar: The minutes before or after sunrise and sunset.
 100         @return: True if it is night else False if day.
 101         """
 102         result = False
 103         delta = datetime.timedelta(minutes=collar)
 104 
 105         if (self.__sunrise - delta) > self.__dateLocal or \
 106                self.__dateLocal > (self.__sunset + delta):
 107             result = True
 108 
 109         return result
 110 
 111     def getSunRiseSet(self):
 112         """
 113         Get the sunrise and sunset.
 114 
 115         @return: A C{datetime} object in a tuple (sunrise, sunset).
 116         """
 117         return self.__sunrise, self.__sunset
 118 
 119     def __determineRiseSet(self):
 120         """
 121         Determine both the sunrise and sunset.
 122         """
 123         year = self.__dateLocal.year
 124         month = self.__dateLocal.month
 125         day = self.__dateLocal.day
 126         # Ephemeris
 127         ephem2000Day = 367 * year - (7 * (year + (month + 9) / 12) / 4) + \
 128                        (275 * month / 9) + day - 730531.5
 129         self.__sunrise = self.__determineRiseOrSet(ephem2000Day, 1)
 130         self.__sunset = self.__determineRiseOrSet(ephem2000Day, -1)
 131 
 132     def __determineRiseOrSet(self, ephem2000Day, rs):
 133         """
 134         Determine either the sunrise or the sunset.
 135 
 136         @param ephem2000Day: The Ephemeris from the beginning of the
 137                              21st century.
 138         @param rs: The factor that determines either sunrise or sunset where
 139                    1 equals sunrise and -1 sunset.
 140         @return: Either the sunrise or sunset as a C{datetime} object.
 141         """
 142         utold = pi
 143         utnew = 0
 144         altitude = self.__ZENITH[self.__zenith]
 145         sinAlt = sin(radians(altitude))       # solar altitude
 146         sinPhi = sin(radians(self.__lat))     # viewer's latitude
 147         cosPhi = cos(radians(self.__lat))     #
 148         lon = radians(self.__lon)             # viewer's longitude
 149         ct = 0
 150         #print rs, ephem2000Day, sinAlt, sinPhi, cosPhi, lon
 151 
 152         while fabs(utold - utnew) > 0.001 and ct < 35:
 153             ct += 1
 154             utold = utnew
 155             days = ephem2000Day + utold / (2 * pi)
 156             t = days / 36525
 157             # The magic numbers are orbital elements of the sun.
 158             l = self.__getRange(4.8949504201433 + 628.331969753199 * t)
 159             g = self.__getRange(6.2400408 + 628.3019501 * t)
 160             ec = 0.033423 * sin(g) + 0.00034907 * sin(2 * g)
 161             lam = l + ec
 162             e = -1 * ec + 0.0430398 * sin(2 * lam) - 0.00092502 * sin(4 * lam)
 163             obl = 0.409093 - 0.0002269 * t
 164             delta = sin(obl) * sin(lam)
 165             delta = atan2(delta, sqrt(1 - delta * delta))
 166             gha = utold - pi + e
 167             cosc = (sinAlt - sinPhi * sin(delta)) / (cosPhi * cos(delta))
 168 
 169             if cosc > 1:
 170                 correction = 0
 171             elif cosc < -1:
 172                 correction = pi
 173             else:
 174                 correction = atan2((sqrt(1 - cosc * cosc)), cosc)
 175 
 176             #print cosc, correction, utold, utnew
 177             utnew = self.__getRange(utold - (gha + lon + rs * correction))
 178 
 179         decimalTime = degrees(utnew) / 15
 180         #print utnew, decimalTime
 181         return self.__get24HourLocalTime(rs, decimalTime)
 182 
 183     def __getRange(self, value):
 184         """
 185         Get the range of the value.
 186 
 187         @param value: The domain.
 188         @return: The resultant range.
 189         """
 190         tmp1 = value / (2.0 * pi)
 191         tmp2 = (2.0 * pi) * (tmp1 - int(tmp1))
 192         if tmp2 < 0.0: tmp2 += (2.0 * pi)
 193         return tmp2
 194 
 195     def __get24HourLocalTime(self, rs, decimalTime):
 196         """
 197         Convert the decimal time into a local time (C{datetime} object)
 198         and correct for a 24 hour clock.
 199 
 200         @param rs: The factor that determines either sunrise or sunset where
 201                    1 equals sunrise and -1 sunset.
 202         @param decimalTime: The decimal time.
 203         @return: The C{datetime} objects set to either sunrise or sunset.
 204         """
 205         decimalTime += self.__offsetUTC
 206         #print decimalTime
 207 
 208         if decimalTime < 0.0:
 209             decimalTime += 24.0
 210         elif decimalTime > 24.0:
 211             decimalTime -= 24.0
 212 
 213         #print decimalTime
 214         hour = int(decimalTime)
 215         tmp = (decimalTime - hour) * 60
 216         minute = int(tmp)
 217         tmp = (tmp - minute) * 60
 218         second = int(tmp)
 219         micro = int(round((tmp - second) * 1000000))
 220         localDT = self.__dateLocal.replace(hour=hour, minute=minute,
 221                                            second=second, microsecond=micro)
 222         return localDT
 223 
 224 
 225 def __getRiseSet(date, lat=35.9513, lon=-83.9142, zenith='official'):
 226     """
 227     The default lat and lon are for Knoxville, TN. The default zenith is
 228     'official'.
 229     """
 230     rs = SunriseSunset(date, lat, lon, zenith=zenith)
 231     riseTime, setTime = rs.getSunRiseSet()
 232     print "Using zenith: %s" % zenith
 233     print "Date/Time now: %s" % date
 234     print "Sunrise: %s" % riseTime
 235     print " Sunset: %s" % setTime
 236     print "Is night: %s\n" % rs.isNight()
 237 
 238 
 239 if __name__ == '__main__':
 240     import sys, pytz
 241     zone = pytz.timezone("US/Eastern")
 242 
 243     # Get sunrise and sunset for now and all the zenith types.
 244     now = datetime.datetime.now(zone)
 245     # Get the zenith types.
 246     zenithKeys = SunriseSunset._SunriseSunset__ZENITH.keys()
 247     print "Test zenith"
 248 
 249     for zenith in zenithKeys:
 250         __getRiseSet(now, zenith=zenith)
 251 
 252     # Get sunrise sunset for every hour of the day using the default zenith.
 253     print "\nTest 24 hours"
 254 
 255     for hour in range(24):
 256         for minute in range(0, 60, 10):
 257             date = now.replace(hour=hour, minute=minute, second=0,
 258                                microsecond=0)
 259             __getRiseSet(date)
 260 
 261     sys.exit(0)