1 | /* |
2 | This file is part of the kholidays library. |
3 | |
4 | Adapted from the Javascript found at http://www.esrl.noaa.gov/gmd/grad/solcalc |
5 | which is in the public domain. |
6 | |
7 | Copyright (c) 2012 Allen Winter <winter@kde.org> |
8 | |
9 | This library is free software; you can redistribute it and/or |
10 | modify it under the terms of the GNU Library General Public |
11 | License as published by the Free Software Foundation; either |
12 | version 2 of the License, or (at your option) any later version. |
13 | |
14 | This library is distributed in the hope that it will be useful, |
15 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
17 | Library General Public License for more details. |
18 | |
19 | You should have received a copy of the GNU Library General Public License |
20 | along with this library; see the file COPYING.LIB. If not, write to |
21 | the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, |
22 | Boston, MA 02110-1301, USA. |
23 | */ |
24 | |
25 | #include "sunriseset.h" |
26 | #include <cmath> |
27 | #ifndef KDE_USE_FINAL |
28 | static double PI = 3.14159265358979323846; |
29 | #endif |
30 | static double MaxLat = 89.99; |
31 | static double MaxLong = 179.99; |
32 | |
33 | using namespace KHolidays; |
34 | using namespace SunRiseSet; |
35 | |
36 | static double degToRad( double degree ) |
37 | { |
38 | return ( degree * PI ) / 180.00; |
39 | } |
40 | |
41 | static double radToDeg( double rad ) |
42 | { |
43 | return ( rad * 180.0 ) / PI; |
44 | } |
45 | |
46 | //convert Julian Day to centuries since J2000.0. |
47 | static double calcTimeJulianCent( int jday ) |
48 | { |
49 | return ( double( jday ) - 2451545.0 ) / 36525.0; |
50 | } |
51 | |
52 | //calculate the mean obliquity of the ecliptic (in degrees) |
53 | static double calcMeanObliquityOfEcliptic( double t ) |
54 | { |
55 | double seconds = 21.448 - t * ( 46.8150 + t * ( 0.00059 - t * 0.001813 ) ); |
56 | double e0 = 23.0 + ( 26.0 + ( seconds / 60.0 ) ) / 60.0; |
57 | return e0; // in degrees |
58 | } |
59 | |
60 | //calculate the corrected obliquity of the ecliptic (in degrees) |
61 | static double calcObliquityCorrection( double t ) |
62 | { |
63 | double e0 = calcMeanObliquityOfEcliptic( t ); |
64 | double omega = 125.04 - 1934.136 * t; |
65 | double e = e0 + 0.00256 * cos( degToRad( omega ) ); |
66 | return e; // in degrees |
67 | } |
68 | |
69 | //calculate the Geometric Mean Longitude of the Sun (in degrees) |
70 | static double calcGeomMeanLongSun( double t ) |
71 | { |
72 | double L0 = 280.46646 + t * ( 36000.76983 + 0.0003032 * t ); |
73 | while ( L0 > 360.0 ) { |
74 | L0 -= 360.0; |
75 | } |
76 | while ( L0 < 0.0 ) { |
77 | L0 += 360.0; |
78 | } |
79 | return L0; // in degrees |
80 | } |
81 | |
82 | //calculate the Geometric Mean Anomaly of the Sun (in degrees) |
83 | static double calcGeomMeanAnomalySun( double t ) |
84 | { |
85 | double M = 357.52911 + t * ( 35999.05029 - 0.0001537 * t ); |
86 | return M; // in degrees |
87 | } |
88 | |
89 | //calculate the equation of center for the sun (in degrees) |
90 | static double calcSunEqOfCenter( double t ) |
91 | { |
92 | double m = calcGeomMeanAnomalySun( t ); |
93 | double mrad = degToRad( m ); |
94 | double sinm = sin( mrad ); |
95 | double sin2m = sin( mrad + mrad ); |
96 | double sin3m = sin( mrad + mrad + mrad ); |
97 | double C = sinm * ( 1.914602 - t * ( 0.004817 + 0.000014 * t ) ) + |
98 | sin2m * ( 0.019993 - 0.000101 * t ) + sin3m * 0.000289; |
99 | return C; // in degrees |
100 | } |
101 | |
102 | //calculate the true longitude of the sun (in degrees) |
103 | static double calcSunTrueLong( double t ) |
104 | { |
105 | double l0 = calcGeomMeanLongSun( t ); |
106 | double c = calcSunEqOfCenter( t ); |
107 | double O = l0 + c; |
108 | return O; // in degrees |
109 | } |
110 | |
111 | //calculate the apparent longitude of the sun (in degrees) |
112 | static double calcSunApparentLong( double t ) |
113 | { |
114 | double o = calcSunTrueLong( t ); |
115 | double omega = 125.04 - 1934.136 * t; |
116 | double lambda = o - 0.00569 - 0.00478 * sin( degToRad( omega ) ); |
117 | return lambda; // in degrees |
118 | } |
119 | |
120 | //calculate the declination of the sun (in degrees) |
121 | static double calcSunDeclination( double t ) |
122 | { |
123 | double e = calcObliquityCorrection( t ); |
124 | double lambda = calcSunApparentLong( t ); |
125 | |
126 | double sint = sin( degToRad( e ) ) * sin( degToRad( lambda ) ); |
127 | double theta = radToDeg( asin( sint ) ); |
128 | return theta; // in degrees |
129 | } |
130 | |
131 | //calculate the eccentricity of earth's orbit (unitless) |
132 | static double calcEccentricityEarthOrbit( double t ) |
133 | { |
134 | double e = 0.016708634 - t * ( 0.000042037 + 0.0000001267 * t ); |
135 | return e; // unitless |
136 | } |
137 | |
138 | //calculate the difference between true solar time and mean solar time |
139 | //(output: equation of time in minutes of time) |
140 | static double calcEquationOfTime( double t ) |
141 | { |
142 | double epsilon = calcObliquityCorrection( t ); |
143 | double l0 = calcGeomMeanLongSun( t ); |
144 | double e = calcEccentricityEarthOrbit( t ); |
145 | double m = calcGeomMeanAnomalySun( t ); |
146 | |
147 | double y = tan( degToRad( epsilon ) / 2.0 ); |
148 | y *= y; |
149 | |
150 | double sin2l0 = sin( 2.0 * degToRad( l0 ) ); |
151 | double sinm = sin( degToRad( m ) ); |
152 | double cos2l0 = cos( 2.0 * degToRad( l0 ) ); |
153 | double sin4l0 = sin( 4.0 * degToRad( l0 ) ); |
154 | double sin2m = sin( 2.0 * degToRad( m ) ); |
155 | |
156 | double Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * |
157 | y * sinm * cos2l0 - 0.5 * |
158 | y * y * sin4l0 - 1.25 * e * e * sin2m; |
159 | return radToDeg( Etime ) * 4.0; // in minutes of time |
160 | } |
161 | |
162 | //calculate the hour angle of the sun at sunrise for the latitude (in radians) |
163 | static double calcHourAngleSunrise( double latitude, double solarDec ) |
164 | { |
165 | double latRad = degToRad( latitude ); |
166 | double sdRad = degToRad( solarDec ); |
167 | double HAarg = ( cos( degToRad( 90.833 ) ) / |
168 | ( cos( latRad ) * cos( sdRad ) ) - tan( latRad ) * tan( sdRad ) ); |
169 | double HA = acos( HAarg ); |
170 | return HA; // in radians (for sunset, use -HA) |
171 | } |
172 | |
173 | QTime KHolidays::SunRiseSet::utcSunrise( const QDate &date, double latitude, double longitude ) |
174 | { |
175 | latitude = qMax( qMin( MaxLat, latitude ), -MaxLat ); |
176 | longitude = qMax( qMin( MaxLong, longitude ), -MaxLong ); |
177 | |
178 | double t = calcTimeJulianCent( date.toJulianDay() ); |
179 | double eqTime = calcEquationOfTime( t ); |
180 | double solarDec = calcSunDeclination( t ); |
181 | double hourAngle = calcHourAngleSunrise( latitude, solarDec ); |
182 | double delta = longitude + radToDeg( hourAngle ); |
183 | QTime timeUTC; |
184 | timeUTC = timeUTC.addSecs( ( 720 - ( 4.0 * delta ) - eqTime ) * 60 ); |
185 | return QTime( timeUTC.hour(), |
186 | timeUTC.second() > 29 ? timeUTC.minute() + 1 : timeUTC.minute(), |
187 | 0 ); |
188 | } |
189 | |
190 | QTime KHolidays::SunRiseSet::utcSunset( const QDate &date, double latitude, double longitude ) |
191 | { |
192 | latitude = qMax( qMin( MaxLat, latitude ), -MaxLat ); |
193 | longitude = qMax( qMin( MaxLong, longitude ), -MaxLong ); |
194 | |
195 | double t = calcTimeJulianCent( date.toJulianDay() ); |
196 | double eqTime = calcEquationOfTime( t ); |
197 | double solarDec = calcSunDeclination( t ); |
198 | double hourAngle = -calcHourAngleSunrise( latitude, solarDec ); |
199 | double delta = longitude + radToDeg( hourAngle ); |
200 | QTime timeUTC; |
201 | timeUTC = timeUTC.addSecs( ( 720 - ( 4.0 * delta ) - eqTime ) * 60 ); |
202 | return QTime( timeUTC.hour(), |
203 | timeUTC.second() > 29 ? timeUTC.minute() + 1 : timeUTC.minute(), |
204 | 0 ); |
205 | } |
206 | |