1 /*
2 This file is part of Caelum.
3 See http://www.ogre3d.org/wiki/index.php/Caelum
4 
5 Copyright (c) 2008 Caelum team. See Contributors.txt for details.
6 
7 Caelum is free software: you can redistribute it and/or modify
8 it under the terms of the GNU Lesser General Public License as published
9 by the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 Caelum is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 GNU Lesser General Public License for more details.
16 
17 You should have received a copy of the GNU Lesser General Public License
18 along with Caelum. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #include "CaelumPrecompiled.h"
22 #include "Astronomy.h"
23 
24 namespace Caelum
25 {
26     const LongReal Astronomy::PI = 3.1415926535897932384626433832795029L;
27 
28     const LongReal Astronomy::J2000 = 2451545.0;
29 
radToDeg(LongReal value)30     LongReal Astronomy::radToDeg (LongReal value)
31     {
32         return value * 180 / PI;
33     }
34 
degToRad(LongReal value)35     LongReal Astronomy::degToRad (LongReal value)
36     {
37         return value * PI / 180;
38     }
39 
sinDeg(LongReal x)40     LongReal Astronomy::sinDeg (LongReal x) {
41         return std::sin (degToRad (x));
42     }
43 
cosDeg(LongReal x)44     LongReal Astronomy::cosDeg (LongReal x) {
45         return std::cos (degToRad (x));
46     }
47 
atan2Deg(LongReal y,LongReal x)48     LongReal Astronomy::atan2Deg (LongReal y, LongReal x) {
49         return radToDeg(std::atan2 (y, x));
50     }
51 
normalizeDegrees(LongReal value)52     LongReal Astronomy::normalizeDegrees (LongReal value)
53     {
54         value = fmod (value, 360);
55         if (value < LongReal (0)) {
56             value += LongReal (360);
57         }
58         return value;
59     }
60 
convertEclipticToEquatorialRad(LongReal lon,LongReal lat,LongReal & rasc,LongReal & decl)61 	void Astronomy::convertEclipticToEquatorialRad (
62             LongReal lon, LongReal lat,
63             LongReal &rasc, LongReal &decl)
64 	{
65 		double ecl = Astronomy::degToRad(23.439281);
66 
67 		double x = cos(lon) * cos(lat);
68 		double y = cos(ecl) * sin(lon) * cos(lat) - sin(ecl) * sin(lat);
69 		double z = sin(ecl) * sin(lon) * cos(lat) + cos(ecl) * sin(lat);
70 
71         double r = sqrt(x * x + y * y);
72         rasc = atan2(y, x);
73         decl = atan2(z, r);
74 	}
75 
convertRectangularToSpherical(LongReal x,LongReal y,LongReal z,LongReal & rasc,LongReal & decl,LongReal & dist)76     void Astronomy::convertRectangularToSpherical (
77             LongReal x, LongReal y, LongReal z,
78             LongReal &rasc, LongReal &decl, LongReal &dist)
79     {
80         dist = sqrt (x * x + y * y + z * z);
81         rasc = atan2Deg (y, x);
82         decl = atan2Deg (z, sqrt (x * x + y * y));
83     }
84 
convertSphericalToRectangular(LongReal rasc,LongReal decl,LongReal dist,LongReal & x,LongReal & y,LongReal & z)85     void Astronomy::convertSphericalToRectangular (
86             LongReal rasc, LongReal decl, LongReal dist,
87             LongReal &x, LongReal &y, LongReal &z)
88     {
89         x = dist * cosDeg (rasc) * cosDeg (decl);
90         y = dist * sinDeg (rasc) * cosDeg (decl);
91         z = dist * sinDeg (decl);
92     }
93 
convertEquatorialToHorizontal(LongReal jday,LongReal longitude,LongReal latitude,LongReal rasc,LongReal decl,LongReal & azimuth,LongReal & altitude)94 	void Astronomy::convertEquatorialToHorizontal (
95             LongReal jday,
96             LongReal longitude,   LongReal latitude,
97             LongReal rasc,        LongReal decl,
98             LongReal &azimuth,    LongReal &altitude)
99     {
100         LongReal d = jday - 2451543.5;
101         LongReal w = LongReal (282.9404 + 4.70935E-5 * d);
102         LongReal M = LongReal (356.0470 + 0.9856002585 * d);
103         // Sun's mean longitude
104         LongReal L = w + M;
105         // Universal time of day in degrees.
106         LongReal UT = LongReal(fmod(d, 1) * 360);
107         LongReal hourAngle = longitude + L + LongReal (180) + UT - rasc;
108 
109         LongReal x = cosDeg (hourAngle) * cosDeg (decl);
110         LongReal y = sinDeg (hourAngle) * cosDeg (decl);
111         LongReal z = sinDeg (decl);
112 
113         LongReal xhor = x * sinDeg (latitude) - z * cosDeg (latitude);
114         LongReal yhor = y;
115         LongReal zhor = x * cosDeg (latitude) + z * sinDeg (latitude);
116 
117         azimuth = atan2Deg (yhor, xhor) + LongReal (180);
118         altitude = atan2Deg (zhor, sqrt (xhor * xhor + yhor * yhor));
119     }
120 
getHorizontalSunPosition(LongReal jday,LongReal longitude,LongReal latitude,LongReal & azimuth,LongReal & altitude)121     void Astronomy::getHorizontalSunPosition (
122             LongReal jday,
123             LongReal longitude, LongReal latitude,
124             LongReal &azimuth, LongReal &altitude)
125     {
126         // 2451544.5 == Astronomy::getJulianDayFromGregorianDateTime(2000, 1, 1, 0, 0, 0));
127         // 2451543.5 == Astronomy::getJulianDayFromGregorianDateTime(1999, 12, 31, 0, 0, 0));
128         LongReal d = jday - 2451543.5;
129 
130         // Sun's Orbital elements:
131         // argument of perihelion
132         LongReal w = LongReal (282.9404 + 4.70935E-5 * d);
133         // eccentricity (0=circle, 0-1=ellipse, 1=parabola)
134         LongReal e = 0.016709 - 1.151E-9 * d;
135         // mean anomaly (0 at perihelion; increases uniformly with time)
136         LongReal M = LongReal(356.0470 + 0.9856002585 * d);
137         // Obliquity of the ecliptic.
138         //LongReal oblecl = LongReal (23.4393 - 3.563E-7 * d);
139 
140         // Eccentric anomaly
141         LongReal E = M + radToDeg(e * sinDeg (M) * (1 + e * cosDeg (M)));
142 
143         // Sun's Distance(R) and true longitude(L)
144         LongReal xv = cosDeg (E) - e;
145         LongReal yv = sinDeg (E) * sqrt (1 - e * e);
146         //LongReal r = sqrt (xv * xv + yv * yv);
147         LongReal lon = atan2Deg (yv, xv) + w;
148         LongReal lat = 0;
149 
150 		LongReal lambda = degToRad(lon);
151 		LongReal beta = degToRad(lat);
152         LongReal rasc, decl;
153 		convertEclipticToEquatorialRad (lambda, beta, rasc, decl);
154 		rasc = radToDeg(rasc);
155 		decl = radToDeg(decl);
156 
157         // Horizontal spherical.
158         Astronomy::convertEquatorialToHorizontal (
159                 jday, longitude, latitude, rasc, decl, azimuth, altitude);
160     }
161 
getHorizontalSunPosition(LongReal jday,Ogre::Degree longitude,Ogre::Degree latitude,Ogre::Degree & azimuth,Ogre::Degree & altitude)162     void Astronomy::getHorizontalSunPosition (
163             LongReal jday,
164             Ogre::Degree longitude, Ogre::Degree latitude,
165             Ogre::Degree &azimuth, Ogre::Degree &altitude)
166     {
167         LongReal az, al;
168         getHorizontalSunPosition(jday, longitude.valueDegrees (), latitude.valueDegrees (), az, al);
169         azimuth = Ogre::Degree(az);
170         altitude = Ogre::Degree(al);
171     }
172 
getEclipticMoonPositionRad(LongReal jday,LongReal & lon,LongReal & lat)173 	void Astronomy::getEclipticMoonPositionRad (
174             LongReal jday,
175             LongReal &lon, LongReal &lat)
176 	{
177         // Julian centuries since January 1, 2000
178 		double T = (jday - 2451545.0L) / 36525.0L;
179 		double lprim = 3.8104L + 8399.7091L * T;
180 		double mprim = 2.3554L + 8328.6911L * T;
181 		double m = 6.2300L + 648.3019L * T;
182 		double d = 5.1985L + 7771.3772L * T;
183 		double f = 1.6280L + 8433.4663L * T;
184 		lon = lprim
185                 + 0.1098L * sin(mprim)
186                 + 0.0222L * sin(2.0L * d - mprim)
187                 + 0.0115L * sin(2.0L * d)
188 				+ 0.0037L * sin(2.0L * mprim)
189                 - 0.0032L * sin(m)
190                 - 0.0020L * sin(2.0L * f)
191                 + 0.0010L * sin(2.0L * d - 2.0L * mprim)
192 				+ 0.0010L * sin(2.0L * d - m - mprim)
193                 + 0.0009L * sin(2.0L * d + mprim)
194                 + 0.0008L * sin(2.0L * d - m)
195 				+ 0.0007L * sin(mprim - m)
196                 - 0.0006L * sin(d)
197                 - 0.0005L * sin(m + mprim);
198 		lat =
199                 + 0.0895L * sin(f)
200                 + 0.0049L * sin(mprim + f)
201                 + 0.0048L * sin(mprim - f)
202                 + 0.0030L * sin(2.0L * d - f)
203                 + 0.0010L * sin(2.0L * d + f - mprim)
204                 + 0.0008  * sin(2.0L * d - f - mprim)
205                 + 0.0006L * sin(2.0L * d + f);
206 	}
207 
getHorizontalMoonPosition(LongReal jday,LongReal longitude,LongReal latitude,LongReal & azimuth,LongReal & altitude)208     void Astronomy::getHorizontalMoonPosition (
209             LongReal jday,
210             LongReal longitude, LongReal latitude,
211             LongReal &azimuth, LongReal &altitude)
212     {
213         // Ecliptic spherical
214 		LongReal lonecl, latecl;
215 		Astronomy::getEclipticMoonPositionRad (jday, lonecl, latecl);
216 
217 		// Equatorial spherical
218         LongReal rasc, decl;
219 		Astronomy::convertEclipticToEquatorialRad (lonecl, latecl, rasc, decl);
220 
221 		// Radians to degrees (all angles are in radians up to this point)
222         rasc = radToDeg(rasc);
223 		decl = radToDeg(decl);
224 
225 		// Equatorial to horizontal
226         Astronomy::convertEquatorialToHorizontal (
227                 jday, longitude, latitude, rasc, decl, azimuth, altitude);
228     }
229 
getHorizontalMoonPosition(LongReal jday,Ogre::Degree longitude,Ogre::Degree latitude,Ogre::Degree & azimuth,Ogre::Degree & altitude)230     void Astronomy::getHorizontalMoonPosition (
231             LongReal jday,
232             Ogre::Degree longitude, Ogre::Degree latitude,
233             Ogre::Degree &azimuth, Ogre::Degree &altitude)
234     {
235         LongReal az, al;
236         getHorizontalMoonPosition(jday, longitude.valueDegrees (), latitude.valueDegrees (), az, al);
237         azimuth = Ogre::Degree(az);
238         altitude = Ogre::Degree(al);
239     }
240 
getHorizontalNorthEclipticPolePosition(LongReal jday,Ogre::Degree longitude,Ogre::Degree latitude,Ogre::Degree & azimuth,Ogre::Degree & altitude)241     void Astronomy::getHorizontalNorthEclipticPolePosition (
242 			LongReal jday,
243 			Ogre::Degree longitude, Ogre::Degree latitude,
244 			Ogre::Degree &azimuth, Ogre::Degree &altitude)
245     {
246 		// Get precise equatorial spherical coordinates of North Ecliptic Pole
247         LongReal rasc = 270.0;
248 		LongReal decl = 90.0 - 23.439281;
249 
250 		// Equatorial to horizontal
251         LongReal az, al;
252         Astronomy::convertEquatorialToHorizontal (
253 				jday, longitude.valueDegrees (), latitude.valueDegrees (), rasc, decl, az, al);
254         azimuth = Ogre::Degree(az);
255         altitude = Ogre::Degree(al);
256     }
257 
getJulianDayFromGregorianDate(int year,int month,int day)258     int Astronomy::getJulianDayFromGregorianDate(
259             int year, int month, int day)
260     {
261         // Formulas from http://en.wikipedia.org/wiki/Julian_day
262         // These are all integer divisions, but I'm not sure it works
263         // correctly for negative values.
264         int a = (14 - month) / 12;
265         int y = year + 4800 - a;
266         int m = month + 12 * a - 3;
267         return day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045;
268     }
269 
getJulianDayFromGregorianDateTime(int year,int month,int day,int hour,int minute,LongReal second)270     LongReal Astronomy::getJulianDayFromGregorianDateTime(
271             int year, int month, int day,
272             int hour, int minute, LongReal second)
273     {
274         ScopedHighPrecissionFloatSwitch precissionSwitch;
275 
276         int jdn = getJulianDayFromGregorianDate (year, month, day);
277         // These are NOT integer divisions.
278         LongReal jd = jdn + (hour - 12) / 24.0 + minute / 1440.0 + second / 86400.0;
279 
280         return jd;
281     }
282 
getJulianDayFromGregorianDateTime(int year,int month,int day,LongReal secondsFromMidnight)283     LongReal Astronomy::getJulianDayFromGregorianDateTime(
284             int year, int month, int day,
285             LongReal secondsFromMidnight)
286     {
287         int jdn = getJulianDayFromGregorianDate(year, month, day);
288         LongReal jd = jdn + secondsFromMidnight / 86400.0 - 0.5;
289         return jd;
290     }
291 
getGregorianDateFromJulianDay(int julianDay,int & year,int & month,int & day)292     void Astronomy::getGregorianDateFromJulianDay(
293             int julianDay, int &year, int &month, int &day)
294     {
295 
296         // From http://en.wikipedia.org/wiki/Julian_day
297         int J = julianDay;
298         int j = J + 32044;
299         int g = j / 146097;
300         int dg = j % 146097;
301         int c = (dg / 36524 + 1) * 3 / 4;
302         int dc = dg - c * 36524;
303         int b = dc / 1461;
304         int db = dc % 1461;
305         int a = (db / 365 + 1) * 3 / 4;
306         int da = db - a * 365;
307         int y = g * 400 + c * 100 + b * 4 + a;
308         int m = (da * 5 + 308) / 153 - 2;
309         int d = da - (m + 4) * 153 / 5 + 122;
310         year = y - 4800 + (m + 2) / 12;
311         month = (m + 2) % 12 + 1;
312         day = d + 1;
313     }
314 
getGregorianDateTimeFromJulianDay(LongReal julianDay,int & year,int & month,int & day,int & hour,int & minute,LongReal & second)315     void Astronomy::getGregorianDateTimeFromJulianDay(
316             LongReal julianDay, int &year, int &month, int &day,
317             int &hour, int &minute, LongReal &second)
318     {
319         // Integer julian days are at noon.
320         // static_cast<int)(floor( is more precise than Ogre::Math::IFloor.
321         // Yes, it does matter.
322         int fpmode = enterHighPrecissionFloatingPointMode();
323         julianDay += (LongReal)0.5;
324        int ijd = static_cast<int>(floor(julianDay));
325         getGregorianDateFromJulianDay(ijd, year, month, day);
326 
327 
328       LongReal s = (julianDay - (LongReal)ijd);
329       s *= 86400.0;
330         hour = static_cast<int>(floor(s / 3600));
331         s -= hour * 3600;
332         minute = static_cast<int>(floor(s / 60));
333         s -= minute * 60;
334         second = s;
335       restoreFloatingPointMode(fpmode);
336     }
337 
getGregorianDateFromJulianDay(LongReal julianDay,int & year,int & month,int & day)338     void Astronomy::getGregorianDateFromJulianDay(
339             LongReal julianDay, int &year, int &month, int &day)
340     {
341         int hour;
342         int minute;
343         LongReal second;
344         getGregorianDateTimeFromJulianDay(julianDay, year, month, day, hour, minute, second);
345     }
346 
347 #if (OGRE_PLATFORM == OGRE_PLATFORM_WIN32) && (OGRE_COMPILER == OGRE_COMPILER_MSVC)
enterHighPrecissionFloatingPointMode()348     int Astronomy::enterHighPrecissionFloatingPointMode ()
349     {
350         int oldMode = ::_controlfp (0, 0);
351         ::_controlfp (_PC_64, _MCW_PC);
352         return oldMode;
353     }
354 
restoreFloatingPointMode(int oldMode)355     void Astronomy::restoreFloatingPointMode (int oldMode)
356     {
357         ::_controlfp (oldMode, _MCW_PC);
358     }
359 #else
enterHighPrecissionFloatingPointMode()360     int Astronomy::enterHighPrecissionFloatingPointMode ()
361     {
362         // Meaningless
363         return 0xC0FFEE;
364     }
365 
restoreFloatingPointMode(int oldMode)366     void Astronomy::restoreFloatingPointMode (int oldMode)
367     {
368         // Useless check.
369         assert(oldMode == 0xC0FFEE);
370     }
371 #endif
372 }
373