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