3 /* Copyright (C) 2002 Brad Jorsch <anomie@users.sourceforge.net>
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20 /* Algorithms from http://www.srrb.noaa.gov/highlights/sunrise/azel.html */
30 /* Purpose: calculate the Geometric Mean Longitude of the Sun (degrees) */
31 double calcGeomMeanLongSun(double t
) {
32 double L0
= 280.46646 + t
* (36000.76983 + 0.0003032 * t
);
43 /* Purpose: calculate the Geometric Mean Anomaly of the Sun (degrees) */
44 double calcGeomMeanAnomalySun(double t
) {
45 return 357.52911 + t
* (35999.05029 - 0.0001537 * t
);
49 /* Purpose: calculate the eccentricity of earth's orbit */
50 double calcEccentricityEarthOrbit(double t
) {
51 return 0.016708634 - t
* (0.000042037 + 0.0000001267 * t
);
55 /* Purpose: calculate the equation of center for the sun (degrees) */
56 double calcSunEqOfCenter(double t
) {
57 double m
= deg2rad(calcGeomMeanAnomalySun(t
));
59 return sin(m
) * (1.914602 - t
* (0.004817 + 0.000014 * t
)) + sin(m
+m
) * (0.019993 - 0.000101 * t
) + sin(m
+m
+m
) * 0.000289;
63 /* Purpose: calculate the true longitude of the sun (degrees) */
64 double calcSunTrueLong(double t
) {
65 return calcGeomMeanLongSun(t
) + calcSunEqOfCenter(t
);
69 /* Purpose: calculate the apparent longitude of the sun (degrees) */
70 double calcSunApparentLong(double t
) {
71 return calcSunTrueLong(t
) - 0.00569 - 0.00478 * sin(deg2rad(125.04-1934.136*t
));
75 /* Purpose: calculate the mean obliquity of the ecliptic (degrees) */
76 double calcMeanObliquityOfEcliptic(double t
) {
77 return 23.0 + (26.0 + ((21.448 - t
*(46.8150 + t
*(0.00059 - t
*(0.001813))))/60.0))/60.0;
81 /* Purpose: calculate the corrected obliquity of the ecliptic (degrees) */
82 double calcObliquityCorrection(double t
) {
83 return calcMeanObliquityOfEcliptic(t
) + 0.00256*cos(deg2rad(125.04-1934.136*t
));
87 /* Purpose: calculate the declination of the sun (degrees) */
88 double calcSunDeclination(double t
) {
89 return rad2deg(asin(sin(deg2rad(calcObliquityCorrection(t
))) *
90 sin(deg2rad(calcSunApparentLong(t
)))));
94 /* Purpose: calculate the difference between true solar time and mean
95 * solar time (minutes)
97 double calcEquationOfTime(double t
) {
98 double l0
= deg2rad(calcGeomMeanLongSun(t
));
99 double e
= calcEccentricityEarthOrbit(t
);
100 double m
= deg2rad(calcGeomMeanAnomalySun(t
));
101 double y
= tan(deg2rad(calcObliquityCorrection(t
))/2.0);
102 double sinm
= sin(m
);
106 return rad2deg(y
*sin(2.0*l0
) - 2.0*e
*sinm
+ 4.0*e
*y
*sinm
*cos(2.0*l0
)
107 - 0.5*y
*y
*sin(4.0*l0
) - 1.25*e
*e
*sin(2.0*m
))*4.0;
111 double calcSolarZenith(double latitude
, double longitude
, int year
, int month
, int day
, int timeUTC
){
112 double T
, trueSolarTime
, hourAngle
, solarDec
, csz
, zenith
, exoatmElevation
, te
, refractionCorrection
;
114 T
=jd2jcent(mdy2jd(year
, month
, day
) + timeUTC
/1440.0);
115 trueSolarTime
= timeUTC
+ calcEquationOfTime(T
) - 4.0 * longitude
;
116 hourAngle
= trueSolarTime
/ 4.0 - 180.0;
117 solarDec
= calcSunDeclination(T
);
118 csz
= sin(deg2rad(latitude
)) * sin(deg2rad(solarDec
)) +
119 cos(deg2rad(latitude
)) * cos(deg2rad(solarDec
)) *
120 cos(deg2rad(hourAngle
));
121 zenith
=rad2deg(acos(csz
));
122 exoatmElevation
= 90.0 - zenith
;
123 if (exoatmElevation
> 85.0) {
124 refractionCorrection
= 0.0;
126 te
= tan(deg2rad(exoatmElevation
));
127 if (exoatmElevation
> 5.0) {
128 refractionCorrection
= 58.1/te
- 0.07/(te
*te
*te
) +
129 0.000086/(te
*te
*te
*te
*te
);
130 } else if (exoatmElevation
> -0.575) {
131 refractionCorrection
= 1735.0 + exoatmElevation
*(-518.2 + exoatmElevation
*(103.4 + exoatmElevation
*(-12.79 + exoatmElevation
*0.711)));
133 refractionCorrection
= -20.774 / te
;
135 refractionCorrection
= refractionCorrection
/ 3600.0;
137 return zenith
- refractionCorrection
;