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
28static double PI = 3.14159265358979323846;
29#endif
30static double MaxLat = 89.99;
31static double MaxLong = 179.99;
32
33using namespace KHolidays;
34using namespace SunRiseSet;
35
36static double degToRad( double degree )
37{
38 return ( degree * PI ) / 180.00;
39}
40
41static double radToDeg( double rad )
42{
43 return ( rad * 180.0 ) / PI;
44}
45
46//convert Julian Day to centuries since J2000.0.
47static 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)
53static 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)
61static 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)
70static 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)
83static 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)
90static 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)
103static 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)
112static 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)
121static 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)
132static 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)
140static 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)
163static 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
173QTime 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
190QTime 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