1
2
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
35
37 """"""
38
39
40 self.RADEG = 180.0 / pi
41 self.DEGRAD = pi / 180.0
42 self.INV360 = 1.0 / 360.0
43
44
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
53 """Returns the sin in degrees"""
54 return math.sin(x * self.DEGRAD)
55
57 """Returns the cos in degrees"""
58 return math.cos(x * self.DEGRAD)
59
61 """Returns the tan in degrees"""
62 return math.tan(x * self.DEGRAD)
63
65 """Returns the arc tan in degrees"""
66 return math.atan(x) * self.RADEG
67
69 """Returns the arc sin in degrees"""
70 return math.asin(x) * self.RADEG
71
73 """Returns the arc cos in degrees"""
74 return math.acos(x) * self.RADEG
75
77 """Returns the atan2 in degrees"""
78 return math.atan2(y, x) * self.RADEG
79
80
81
82
83
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
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
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
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
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
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
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
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
185 d = self.daysSince2000Jan0(year,month,day) + 0.5 - (lon/360.0)
186
187
188 sidtime = self.revolution(self.GMST0(d) + 180.0 + lon)
189
190
191 res = self.sunRADec(d)
192 sRA = res[0]
193 sdec = res[1]
194 sr = res[2]
195
196
197 tsouth = 12.0 - self.rev180(sidtime - sRA)/15.0;
198
199
200 sradius = 0.2666 / sr;
201
202
203 if upper_limb:
204 altit = altit - sradius
205
206
207
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
215
216 elif cost <= -1.0:
217 rc = +1
218 t = 12.0;
219
220 else:
221 t = self.acosd(cost)/15.0
222
223
224
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
248 d = self.daysSince2000Jan0(year,month,day) + 0.5 - (lon/360.0)
249
250
251 obl_ecl = 23.4393 - 3.563E-7 * d
252
253
254 res = self.sunpos(d)
255 slon = res[0]
256 sr = res[1]
257
258
259 sin_sdecl = self.sind(obl_ecl) * self.sind(slon)
260 cos_sdecl = math.sqrt(1.0 - sin_sdecl * sin_sdecl)
261
262
263 sradius = 0.2666 / sr
264
265
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
274
275 elif cost <= -1.0:
276 t = 24.0
277
278 else:
279 t = (2.0/15.0) * self.acosd(cost);
280
281 return t
282
283
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
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
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)
302 v = self.atan2d(y, x)
303 lon = v + w
304 if lon >= 360.0:
305 lon = lon - 360.0
306
307 return (lon,r)
308
309
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
318 res = self.sunpos(d)
319 lon = res[0]
320 r = res[1]
321
322
323 x = r * self.cosd(lon)
324 y = r * self.sind(lon)
325
326
327 obl_ecl = 23.4393 - 3.563E-7 * d
328
329
330 z = y * self.sind(obl_ecl)
331 y = y * self.cosd(obl_ecl)
332
333
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
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
351 """
352 Reduce angle to within +180..+180 degrees
353 """
354 return (x - 360.0 * math.floor(x * self.INV360 + 0.5))
355
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
382
383
384
385
386
387 sidtim0 = self.revolution((180.0 + 356.0470 + 282.9404) +
388 (0.9856002585 + 4.70935E-5) * d)
389 return sidtim0;
390
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
405 N = self.daysSince2000Jan0(year, month, day)
406 res = self.sunRADec(N)
407 declination = res[1]
408 sr = res[2]
409
410
411 altitude = 90.0 - latitude + declination
412
413
414 if altitude > 90:
415 altitude = 90 - (altitude-90)
416
417 if altitude < 0:
418 altitude = 0
419
420 return altitude
421
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
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
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
471 nJulianDate = self.Julian(year, month, day)
472
473 if(calendar.isleap(year)):
474 fDivide = 366.0
475 else:
476 fDivide = 365.0
477
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
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
490 fEot = fEot/60.0
491
492 fEot = fEot*15*pi/180.0
493
494 return (fEot, fR0r, tDeclsc)
495
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):
519 lMonth = [0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366]
520 else:
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
532