Module Sun
[hide private]
[frames] | no frames]

Source Code for Module Sun

  1  #!/usr/bin/env python  
  2  # -*- coding: iso-8859-1 -*- 
  3  """  
  4  SUNRISET.C computes Sun rise/set times, start/end of twilight, and 
  5  the length of the day at any date and latitude 
  6                            
  7  Written as DAYLEN.C, 1989-08-16 
  8   
  9  Modified to SUNRISET.C, 1992-12-01 
 10                            
 11  (c) Paul Schlyter, 1989, 1992 
 12                            
 13  Released to the public domain by 
 14  U{Paul Schlyter<http://stjarnhimlen.se/english.html>}, December 1992 
 15                            
 16  Direct conversion to Java  
 17  U{Sean Russell<ser@germane-software.com>} 
 18   
 19  Conversion to Python Class, 2002-03-21 
 20  U{Henrik Härkönen<radix@kortis.to>} 
 21   
 22   - Solar Altitude added by Miguel Tremblay 2005-01-16 
 23   - Solar flux, equation of time and import of python library 
 24     added by U{Miguel Tremblay<http://ptaff.ca/miguel/>} 2007-11-22 
 25   
 26  Source available at U{Henrik's Python page<http://kortis.to/radix/python/>} 
 27  """ 
 28   
 29  import math 
 30  from math import pi 
 31   
 32  import calendar 
 33   
34 -class Sun:
35
36 - def __init__(self):
37 """""" 38 39 # Some conversion factors between radians and degrees 40 self.RADEG = 180.0 / pi 41 self.DEGRAD = pi / 180.0 42 self.INV360 = 1.0 / 360.0
43 44
45 - def daysSince2000Jan0(self, y, m, d):
46 """A macro to compute the number of days elapsed since 2000 Jan 0.0 47 (which is equal to 1999 Dec 31, 0h UT)""" 48 return (367*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530)
49 50 51 # The trigonometric functions in degrees
52 - def sind(self, x):
53 """Returns the sin in degrees""" 54 return math.sin(x * self.DEGRAD)
55
56 - def cosd(self, x):
57 """Returns the cos in degrees""" 58 return math.cos(x * self.DEGRAD)
59
60 - def tand(self, x):
61 """Returns the tan in degrees""" 62 return math.tan(x * self.DEGRAD)
63
64 - def atand(self, x):
65 """Returns the arc tan in degrees""" 66 return math.atan(x) * self.RADEG
67
68 - def asind(self, x):
69 """Returns the arc sin in degrees""" 70 return math.asin(x) * self.RADEG
71
72 - def acosd(self, x):
73 """Returns the arc cos in degrees""" 74 return math.acos(x) * self.RADEG
75
76 - def atan2d(self, y, x):
77 """Returns the atan2 in degrees""" 78 return math.atan2(y, x) * self.RADEG
79 80 # Following are some macros around the "workhorse" function __daylen__ 81 # They mainly fill in the desired values for the reference altitude 82 # below the horizon, and also selects whether this altitude should 83 # refer to the Sun's center or its upper limb. 84
85 - def dayLength(self, year, month, day, lon, lat):
86 """ 87 This macro computes the length of the day, from sunrise to sunset. 88 Sunrise/set is considered to occur when the Sun's upper limb is 89 35 arc minutes below the horizon (this accounts for the refraction 90 of the Earth's atmosphere). 91 """ 92 return self.__daylen__(year, month, day, lon, lat, -35.0/60.0, 1)
93
94 - def dayCivilTwilightLength(self, year, month, day, lon, lat):
95 """ 96 This macro computes the length of the day, including civil twilight. 97 Civil twilight starts/ends when the Sun's center is 6 degrees below 98 the horizon. 99 """ 100 return self.__daylen__(year, month, day, lon, lat, -6.0, 0)
101
102 - def dayNauticalTwilightLength(self, year, month, day, lon, lat):
103 """ 104 This macro computes the length of the day, incl. nautical twilight. 105 Nautical twilight starts/ends when the Sun's center is 12 degrees 106 below the horizon. 107 """ 108 return self.__daylen__(year, month, day, lon, lat, -12.0, 0)
109
110 - def dayAstronomicalTwilightLength(self, year, month, day, lon, lat):
111 """ 112 This macro computes the length of the day, incl. astronomical twilight. 113 Astronomical twilight starts/ends when the Sun's center is 18 degrees 114 below the horizon. 115 """ 116 return self.__daylen__(year, month, day, lon, lat, -18.0, 0)
117
118 - def sunRiseSet(self, year, month, day, lon, lat):
119 """ 120 This macro computes times for sunrise/sunset. 121 Sunrise/set is considered to occur when the Sun's upper limb is 122 35 arc minutes below the horizon (this accounts for the refraction 123 of the Earth's atmosphere). 124 """ 125 return self.__sunriset__(year, month, day, lon, lat, -35.0/60.0, 1)
126 127
128 - def civilTwilight(self, year, month, day, lon, lat):
129 """ 130 This macro computes the start and end times of civil twilight. 131 Civil twilight starts/ends when the Sun's center is 6 degrees below 132 the horizon. 133 """ 134 return self.__sunriset__(year, month, day, lon, lat, -6.0, 0)
135
136 - def nauticalTwilight(self, year, month, day, lon, lat):
137 """ 138 This macro computes the start and end times of nautical twilight. 139 Nautical twilight starts/ends when the Sun's center is 12 degrees 140 below the horizon. 141 """ 142 return self.__sunriset__(year, month, day, lon, lat, -12.0, 0)
143
144 - def astronomicalTwilight(self, year, month, day, lon, lat):
145 """ 146 This macro computes the start and end times of astronomical twilight. 147 Astronomical twilight starts/ends when the Sun's center is 18 degrees 148 below the horizon. 149 """ 150 return self.__sunriset__(year, month, day, lon, lat, -18.0, 0)
151 152 # The "workhorse" function for sun rise/set times
153 - def __sunriset__(self, year, month, day, lon, lat, altit, upper_limb):
154 """ 155 Note: 156 - year,month,date = calendar date, 1801-2099 only. 157 Eastern longitude positive, Western longitude negative 158 Northern latitude positive, Southern latitude negative 159 The longitude value IS critical in this function! 160 - altit = the altitude which the Sun should cross 161 Set to -35/60 degrees for rise/set, -6 degrees 162 for civil, -12 degrees for nautical and -18 163 degrees for astronomical twilight. 164 - upper_limb: non-zero -> upper limb, zero -> center 165 Set to non-zero (e.g. 1) when computing rise/set 166 times, and to zero when computing start/end of 167 twilight. 168 - *rise = where to store the rise time 169 - *set = where to store the set time 170 Both times are relative to the specified altitude, 171 and thus this function can be used to compute 172 various twilight times, as well as rise/set times 173 - Return value: 174 - 0 = sun rises/sets this day, times stored at 175 *trise and *tset. 176 - +1 = sun above the specified 'horizon' 24 hours. 177 *trise set to time when the sun is at south, 178 minus 12 hours while *tset is set to the south 179 time plus 12 hours. 'Day' length = 24 hours 180 - -1 = sun is below the specified 'horizon' 24 hours 181 'Day' length = 0 hours, *trise and *tset are 182 both set to the time when the sun is at south. 183 """ 184 # Compute d of 12h local mean solar time 185 d = self.daysSince2000Jan0(year,month,day) + 0.5 - (lon/360.0) 186 187 # Compute local sidereal time of this moment 188 sidtime = self.revolution(self.GMST0(d) + 180.0 + lon) 189 190 # Compute Sun's RA + Decl at this moment 191 res = self.sunRADec(d) 192 sRA = res[0] 193 sdec = res[1] 194 sr = res[2] 195 196 # Compute time when Sun is at south - in hours UT 197 tsouth = 12.0 - self.rev180(sidtime - sRA)/15.0; 198 199 # Compute the Sun's apparent radius, degrees 200 sradius = 0.2666 / sr; 201 202 # Do correction to upper limb, if necessary 203 if upper_limb: 204 altit = altit - sradius 205 206 # Compute the diurnal arc that the Sun traverses to reach 207 # the specified altitude altit: 208 209 cost = (self.sind(altit) - self.sind(lat) * self.sind(sdec))/\ 210 (self.cosd(lat) * self.cosd(sdec)) 211 212 if cost >= 1.0: 213 rc = -1 214 t = 0.0 # Sun always below altit 215 216 elif cost <= -1.0: 217 rc = +1 218 t = 12.0; # Sun always above altit 219 220 else: 221 t = self.acosd(cost)/15.0 # The diurnal arc, hours 222 223 224 # Store rise and set times - in hours UT 225 return (tsouth-t, tsouth+t)
226 227
228 - def __daylen__(self, year, month, day, lon, lat, altit, upper_limb):
229 """ 230 Note: 231 - year,month,date = calendar date, 1801-2099 only. 232 Eastern longitude positive, Western longitude negative 233 Northern latitude positive, Southern latitude negative 234 The longitude value is not critical. Set it to the correct 235 longitude if you're picky, otherwise set to, say, 0.0 236 - The latitude however IS critical - be sure to get it correct 237 - altit = the altitude which the Sun should cross 238 Set to -35/60 degrees for rise/set, -6 degrees 239 for civil, -12 degrees for nautical and -18 240 degrees for astronomical twilight. 241 - upper_limb: non-zero -> upper limb, zero -> center 242 Set to non-zero (e.g. 1) when computing day length 243 and to zero when computing day+twilight length. 244 245 """ 246 247 # Compute d of 12h local mean solar time 248 d = self.daysSince2000Jan0(year,month,day) + 0.5 - (lon/360.0) 249 250 # Compute obliquity of ecliptic (inclination of Earth's axis) 251 obl_ecl = 23.4393 - 3.563E-7 * d 252 253 # Compute Sun's position 254 res = self.sunpos(d) 255 slon = res[0] 256 sr = res[1] 257 258 # Compute sine and cosine of Sun's declination 259 sin_sdecl = self.sind(obl_ecl) * self.sind(slon) 260 cos_sdecl = math.sqrt(1.0 - sin_sdecl * sin_sdecl) 261 262 # Compute the Sun's apparent radius, degrees 263 sradius = 0.2666 / sr 264 265 # Do correction to upper limb, if necessary 266 if upper_limb: 267 altit = altit - sradius 268 269 270 cost = (self.sind(altit) - self.sind(lat) * sin_sdecl) / \ 271 (self.cosd(lat) * cos_sdecl) 272 if cost >= 1.0: 273 t = 0.0 # Sun always below altit 274 275 elif cost <= -1.0: 276 t = 24.0 # Sun always above altit 277 278 else: 279 t = (2.0/15.0) * self.acosd(cost); # The diurnal arc, hours 280 281 return t
282 283
284 - def sunpos(self, d):
285 """ 286 Computes the Sun's ecliptic longitude and distance 287 at an instant given in d, number of days since 288 2000 Jan 0.0. The Sun's ecliptic latitude is not 289 computed, since it's always very near 0. 290 """ 291 292 # Compute mean elements 293 M = self.revolution(356.0470 + 0.9856002585 * d) 294 w = 282.9404 + 4.70935E-5 * d 295 e = 0.016709 - 1.151E-9 * d 296 297 # Compute true longitude and radius vector 298 E = M + e * self.RADEG * self.sind(M) * (1.0 + e * self.cosd(M)) 299 x = self.cosd(E) - e 300 y = math.sqrt(1.0 - e*e) * self.sind(E) 301 r = math.sqrt(x*x + y*y) #Solar distance 302 v = self.atan2d(y, x) # True anomaly 303 lon = v + w # True solar longitude 304 if lon >= 360.0: 305 lon = lon - 360.0 # Make it 0..360 degrees 306 307 return (lon,r)
308 309
310 - def sunRADec(self, d):
311 """ 312 Returns the angle of the Sun (RA) 313 the declination (dec) and the distance of the Sun (r) 314 for a given day d. 315 """ 316 317 # Compute Sun's ecliptical coordinates 318 res = self.sunpos(d) 319 lon = res[0] # True solar longitude 320 r = res[1] # Solar distance 321 322 # Compute ecliptic rectangular coordinates (z=0) 323 x = r * self.cosd(lon) 324 y = r * self.sind(lon) 325 326 # Compute obliquity of ecliptic (inclination of Earth's axis) 327 obl_ecl = 23.4393 - 3.563E-7 * d 328 329 # Convert to equatorial rectangular coordinates - x is unchanged 330 z = y * self.sind(obl_ecl) 331 y = y * self.cosd(obl_ecl) 332 333 # Convert to spherical coordinates 334 RA = self.atan2d(y, x) 335 dec = self.atan2d(z, math.sqrt(x*x + y*y)) 336 337 return (RA, dec, r)
338 339
340 - def revolution(self, x):
341 """ 342 This function reduces any angle to within the first revolution 343 by subtracting or adding even multiples of 360.0 until the 344 result is >= 0.0 and < 360.0 345 346 Reduce angle to within 0..360 degrees 347 """ 348 return (x - 360.0 * math.floor(x * self.INV360))
349
350 - def rev180(self, x):
351 """ 352 Reduce angle to within +180..+180 degrees 353 """ 354 return (x - 360.0 * math.floor(x * self.INV360 + 0.5))
355
356 - def GMST0(self, d):
357 """ 358 This function computes GMST0, the Greenwich Mean Sidereal Time 359 at 0h UT (i.e. the sidereal time at the Greenwhich meridian at 360 0h UT). GMST is then the sidereal time at Greenwich at any 361 time of the day. I've generalized GMST0 as well, and define it 362 as: GMST0 = GMST - UT -- this allows GMST0 to be computed at 363 other times than 0h UT as well. While this sounds somewhat 364 contradictory, it is very practical: instead of computing 365 GMST like: 366 367 GMST = (GMST0) + UT * (366.2422/365.2422) 368 369 where (GMST0) is the GMST last time UT was 0 hours, one simply 370 computes: 371 372 GMST = GMST0 + UT 373 374 where GMST0 is the GMST "at 0h UT" but at the current moment! 375 Defined in this way, GMST0 will increase with about 4 min a 376 day. It also happens that GMST0 (in degrees, 1 hr = 15 degr) 377 is equal to the Sun's mean longitude plus/minus 180 degrees! 378 (if we neglect aberration, which amounts to 20 seconds of arc 379 or 1.33 seconds of time) 380 """ 381 # Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr 382 # L = M + w, as defined in sunpos(). Since I'm too lazy to 383 # add these numbers, I'll let the C compiler do it for me. 384 # Any decent C compiler will add the constants at compile 385 # time, imposing no runtime or code overhead. 386 387 sidtim0 = self.revolution((180.0 + 356.0470 + 282.9404) + 388 (0.9856002585 + 4.70935E-5) * d) 389 return sidtim0;
390
391 - def solar_altitude(self, latitude, year, month, day):
392 """ 393 Compute the altitude of the sun. No atmospherical refraction taken 394 in account. 395 Altitude of the southern hemisphere are given relative to 396 true north. 397 Altitude of the northern hemisphere are given relative to 398 true south. 399 Declination is between 23.5° North and 23.5° South depending 400 on the period of the year. 401 Source of formula for altitude is PhysicalGeography.net 402 http://www.physicalgeography.net/fundamentals/6h.html 403 """ 404 # Compute declination 405 N = self.daysSince2000Jan0(year, month, day) 406 res = self.sunRADec(N) 407 declination = res[1] 408 sr = res[2] 409 410 # Compute the altitude 411 altitude = 90.0 - latitude + declination 412 413 # In the tropical and in extreme latitude, values over 90 may occurs. 414 if altitude > 90: 415 altitude = 90 - (altitude-90) 416 417 if altitude < 0: 418 altitude = 0 419 420 return altitude
421
422 - def get_max_solar_flux(self, latitude, year, month, day):
423 """ 424 Compute the maximal solar flux to reach the ground for this date and 425 latitude. 426 Originaly comes from Environment Canada weather forecast model. 427 Information was of the public domain before release by Environment Canada 428 Output is in W/M^2. 429 """ 430 431 (fEot, fR0r, tDeclsc) = self.equation_of_time(year, month, day, latitude) 432 fSF = (tDeclsc[0]+tDeclsc[1])*fR0r 433 434 # In the case of a negative declinaison, solar flux is null 435 if fSF < 0: 436 fCoeff = 0 437 else: 438 fCoeff = -1.56e-12*fSF**4 + 5.972e-9*fSF**3 -\ 439 8.364e-6*fSF**2 + 5.183e-3*fSF - 0.435 440 441 fSFT = fSF * fCoeff 442 443 if fSFT < 0: 444 fSFT=0 445 446 return fSFT
447
448 - def equation_of_time(self, year, month, day, latitude):
449 """ 450 Subroutine computing the part of the equation of time 451 needed in the computing of the theoritical solar flux 452 Correction originating of the CMC GEM model. 453 454 @type year: int 455 @param year: Year 456 @type month: int 457 @param month: Month 458 @type day: int 459 @param day: Day 460 @type latitude: float 461 @param latitude: latitude to comput the equation of time 462 463 464 @rtype: tuple 465 @return: (double fEot, double fR0r, tuple tDeclsc) 466 - dEot: Correction for the equation of time 467 - dR0r: Corrected solar constant for the equation of time 468 - tDeclsc: Declinaison 469 """ 470 # Julian date 471 nJulianDate = self.Julian(year, month, day) 472 # Check if it is a leap year 473 if(calendar.isleap(year)): 474 fDivide = 366.0 475 else: 476 fDivide = 365.0 477 # Correction for "equation of time" 478 fA = nJulianDate/fDivide*2*pi 479 fR0r = self.__Solcons(fA)*0.1367e4 480 fRdecl = 0.412*math.cos((nJulianDate+10.0)*2.0*pi/fDivide-pi) 481 fDeclsc1 = self.sind(latitude)*math.sin(fRdecl) 482 fDeclsc2 = self.cosd(latitude)*math.cos(fRdecl) 483 tDeclsc = (fDeclsc1, fDeclsc2) 484 # in minutes 485 fEot = 0.002733 -7.343*math.sin(fA)+ .5519*math.cos(fA) \ 486 - 9.47*math.sin(2.0*fA) - 3.02*math.cos(2.0*fA) \ 487 - 0.3289*math.sin(3.*fA) -0.07581*math.cos(3.0*fA) \ 488 -0.1935*math.sin(4.0*fA) -0.1245*math.cos(4.0*fA) 489 # Express in fraction of hour 490 fEot = fEot/60.0 491 # Express in radians 492 fEot = fEot*15*pi/180.0 493 494 return (fEot, fR0r, tDeclsc)
495
496 - def __Solcons(self, dAlf):
497 """ 498 Description: Statement function that calculates the variation of the 499 solar constant as a function of the julian day. (dAlf, in radians) 500 501 @type dAlf: double 502 @param dAlf: Solar constant to correct the excentricity 503 504 @rtype: double 505 @return: Variation of the solar constant 506 """ 507 508 dVar = 1.0/(1.0-9.464e-4*math.sin(dAlf)-0.01671*math.cos(dAlf)- \ 509 + 1.489e-4*math.cos(2.0*dAlf)-2.917e-5*math.sin(3.0*dAlf)- \ 510 + 3.438e-4*math.cos(4.0*dAlf))**2 511 return dVar
512 513
514 - def Julian(self, year, month, day):
515 """ 516 Return julian day. 517 """ 518 if calendar.isleap(year): # Bissextil year, 366 days 519 lMonth = [0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366] 520 else: # Normal year, 365 days 521 lMonth = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365] 522 523 nJulian = lMonth[month-1] + day 524 return nJulian
525 526 527 if __name__ == "__main__": 528 529 k = Sun() 530 print k.get_max_solar_flux(46.2, 2004, 01, 30) 531 # print k.sunRiseSet(2002, 3, 22, 25.42, 62.15) 532