1 // astro.cpp
2 //
3 // Copyright (C) 2001-2009, the Celestia Development Team
4 //
5 // This program is free software; you can redistribute it and/or
6 // modify it under the terms of the GNU General Public License
7 // as published by the Free Software Foundation; either version 2
8 // of the License, or (at your option) any later version.
9
10 #include <cstring>
11 #include <cmath>
12 #include <iomanip>
13 #include <cstdio>
14 #include <time.h>
15 #include <celutil/basictypes.h>
16 #include <celmath/mathlib.h>
17 #include "celestia.h"
18 #include "astro.h"
19 #include <celutil/util.h>
20
21 using namespace std;
22
23 const double astro::speedOfLight = 299792.458; // km/s
24
25 // epoch J2000: 12 UT on 1 Jan 2000
26 const double astro::J2000 = 2451545.0;
27
28 const double astro::G = 6.672e-11; // N m^2 / kg^2
29
30 const double astro::SolarMass = 1.989e30;
31 const double astro::EarthMass = 5.976e24;
32 const double astro::LunarMass = 7.354e22;
33
34 const double astro::SOLAR_IRRADIANCE = 1367.6; // Watts / m^2
35 const double astro::SOLAR_POWER = 3.8462e26; // Watts
36
37 // Angle between J2000 mean equator and the ecliptic plane.
38 // 23 deg 26' 21".448 (Seidelmann, _Explanatory Supplement to the
39 // Astronomical Almanac_ (1992), eqn 3.222-1.
40 const double astro::J2000Obliquity = degToRad(23.4392911);
41
42 static const Quatd ECLIPTIC_TO_EQUATORIAL_ROTATION =
43 Quatd::xrotation(-astro::J2000Obliquity);
44 static const Mat3d ECLIPTIC_TO_EQUATORIAL_MATRIX = ECLIPTIC_TO_EQUATORIAL_ROTATION.toMatrix3();
45
46 // Equatorial to galactic coordinate transformation
47 // North galactic pole at:
48 // RA 12h 51m 26.282s (192.85958 deg)
49 // Dec 27 d 07' 42.01" (27.1283361 deg)
50 // Zero longitude at position angle 122.932
51 // (J2000 coordinates)
52 static const double GALACTIC_NODE = 282.85958;
53 static const double GALACTIC_INCLINATION = 90.0 - 27.1283361;
54 static const double GALACTIC_LONGITUDE_AT_NODE = 32.932;
55
56 static const Quatd EQUATORIAL_TO_GALACTIC_ROTATION =
57 Quatd::zrotation(degToRad(GALACTIC_NODE)) *
58 Quatd::xrotation(degToRad(GALACTIC_INCLINATION)) *
59 Quatd::zrotation(degToRad(-GALACTIC_LONGITUDE_AT_NODE));
60 static const Mat3d EQUATORIAL_TO_GALACTIC_MATRIX = EQUATORIAL_TO_GALACTIC_ROTATION.toMatrix3();
61
62 // epoch B1950: 22:09 UT on 21 Dec 1949
63 #define B1950 2433282.423
64
65 // Difference in seconds between Terrestrial Time and International
66 // Atomic Time
67 static const double dTA = 32.184;
68
69 struct LeapSecondRecord
70 {
71 int seconds;
72 double t;
73 };
74
75
76 // Table of leap second insertions. The leap second always
77 // appears as the last second of the day immediately prior
78 // to the date in the table.
79 static const LeapSecondRecord LeapSeconds[] =
80 {
81 { 10, 2441317.5 }, // 1 Jan 1972
82 { 11, 2441499.5 }, // 1 Jul 1972
83 { 12, 2441683.5 }, // 1 Jan 1973
84 { 13, 2442048.5 }, // 1 Jan 1974
85 { 14, 2442413.5 }, // 1 Jan 1975
86 { 15, 2442778.5 }, // 1 Jan 1976
87 { 16, 2443144.5 }, // 1 Jan 1977
88 { 17, 2443509.5 }, // 1 Jan 1978
89 { 18, 2443874.5 }, // 1 Jan 1979
90 { 19, 2444239.5 }, // 1 Jan 1980
91 { 20, 2444786.5 }, // 1 Jul 1981
92 { 21, 2445151.5 }, // 1 Jul 1982
93 { 22, 2445516.5 }, // 1 Jul 1983
94 { 23, 2446247.5 }, // 1 Jul 1985
95 { 24, 2447161.5 }, // 1 Jan 1988
96 { 25, 2447892.5 }, // 1 Jan 1990
97 { 26, 2448257.5 }, // 1 Jan 1991
98 { 27, 2448804.5 }, // 1 Jul 1992
99 { 28, 2449169.5 }, // 1 Jul 1993
100 { 29, 2449534.5 }, // 1 Jul 1994
101 { 30, 2450083.5 }, // 1 Jan 1996
102 { 31, 2450630.5 }, // 1 Jul 1997
103 { 32, 2451179.5 }, // 1 Jan 1999
104 { 33, 2453736.5 }, // 1 Jan 2006
105 { 34, 2454832.5 }, // 1 Jan 2009
106 };
107
108
109 #if !defined(__GNUC__) || defined(_WIN32)
110 static const char* MonthAbbrList[12] =
111 { "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
112 #endif
113
114
115 static Mat3d equatorialToCelestiald = Mat3d::xrotation(astro::J2000Obliquity);
116 static Mat3f equatorialToCelestial = Mat3f::xrotation((float) astro::J2000Obliquity);
117
118
lumToAbsMag(float lum)119 float astro::lumToAbsMag(float lum)
120 {
121 return (float) (SOLAR_ABSMAG - log(lum) * LN_MAG);
122 }
123
124 // Return the apparent magnitude of a star with lum times solar
125 // luminosity viewed at lyrs light years
lumToAppMag(float lum,float lyrs)126 float astro::lumToAppMag(float lum, float lyrs)
127 {
128 return absToAppMag(lumToAbsMag(lum), lyrs);
129 }
130
absMagToLum(float mag)131 float astro::absMagToLum(float mag)
132 {
133 return (float) exp((SOLAR_ABSMAG - mag) / LN_MAG);
134 }
135
appMagToLum(float mag,float lyrs)136 float astro::appMagToLum(float mag, float lyrs)
137 {
138 return absMagToLum(appToAbsMag(mag, lyrs));
139 }
140
lightYearsToParsecs(float ly)141 float astro::lightYearsToParsecs(float ly)
142 {
143 return ly / (float) LY_PER_PARSEC;
144 }
145
lightYearsToParsecs(double ly)146 double astro::lightYearsToParsecs(double ly)
147 {
148 return ly / (double) LY_PER_PARSEC;
149 }
150
parsecsToLightYears(float pc)151 float astro::parsecsToLightYears(float pc)
152 {
153 return pc * (float) LY_PER_PARSEC;
154 }
155
parsecsToLightYears(double pc)156 double astro::parsecsToLightYears(double pc)
157 {
158 return pc * (double) LY_PER_PARSEC;
159 }
160
lightYearsToKilometers(float ly)161 float astro::lightYearsToKilometers(float ly)
162 {
163 return ly * (float) KM_PER_LY;
164 }
165
lightYearsToKilometers(double ly)166 double astro::lightYearsToKilometers(double ly)
167 {
168 return ly * KM_PER_LY;
169 }
170
kilometersToLightYears(float km)171 float astro::kilometersToLightYears(float km)
172 {
173 return km / (float) KM_PER_LY;
174 }
175
kilometersToLightYears(double km)176 double astro::kilometersToLightYears(double km)
177 {
178 return km / KM_PER_LY;
179 }
180
lightYearsToAU(float ly)181 float astro::lightYearsToAU(float ly)
182 {
183 return ly * (float) AU_PER_LY;
184 }
185
lightYearsToAU(double ly)186 double astro::lightYearsToAU(double ly)
187 {
188 return ly * AU_PER_LY;
189 }
190
AUtoKilometers(float au)191 float astro::AUtoKilometers(float au)
192 {
193 return au * (float) KM_PER_AU;
194 }
195
AUtoKilometers(double au)196 double astro::AUtoKilometers(double au)
197 {
198 return au * (double) KM_PER_AU;
199 }
200
kilometersToAU(float km)201 float astro::kilometersToAU(float km)
202 {
203 return km / (float) KM_PER_AU;
204 }
205
kilometersToAU(double km)206 double astro::kilometersToAU(double km)
207 {
208 return km / KM_PER_AU;
209 }
210
secondsToJulianDate(double sec)211 double astro::secondsToJulianDate(double sec)
212 {
213 return sec / 86400.0;
214 }
215
julianDateToSeconds(double jd)216 double astro::julianDateToSeconds(double jd)
217 {
218 return jd * 86400.0;
219 }
220
decimalToDegMinSec(double angle,int & degrees,int & minutes,double & seconds)221 void astro::decimalToDegMinSec(double angle, int& degrees, int& minutes, double& seconds)
222 {
223 double A, B, C;
224
225 degrees = (int) angle;
226
227 A = angle - (double) degrees;
228 B = A * 60.0;
229 minutes = (int) B;
230 C = B - (double) minutes;
231 seconds = C * 60.0;
232 }
233
degMinSecToDecimal(int degrees,int minutes,double seconds)234 double astro::degMinSecToDecimal(int degrees, int minutes, double seconds)
235 {
236 return (double)degrees + (seconds/60.0 + (double)minutes)/60.0;
237 }
238
decimalToHourMinSec(double angle,int & hours,int & minutes,double & seconds)239 void astro::decimalToHourMinSec(double angle, int& hours, int& minutes, double& seconds)
240 {
241 double A, B;
242
243 A = angle / 15.0;
244 hours = (int) A;
245 B = (A - (double) hours) * 60.0;
246 minutes = (int) (B);
247 seconds = (B - (double) minutes) * 60.0;
248 }
249
250 // Compute the fraction of a sphere which is illuminated and visible
251 // to a viewer. The source of illumination is assumed to be at (0, 0, 0)
sphereIlluminationFraction(Point3d,Point3d)252 float astro::sphereIlluminationFraction(Point3d,
253 Point3d)
254 {
255 return 1.0f;
256 }
257
microLightYearsToKilometers(float ly)258 float astro::microLightYearsToKilometers(float ly)
259 {
260 return ly * ((float) KM_PER_LY * 1e-6f);
261 }
262
microLightYearsToKilometers(double ly)263 double astro::microLightYearsToKilometers(double ly)
264 {
265 return ly * (KM_PER_LY * 1e-6);
266 }
267
kilometersToMicroLightYears(float km)268 float astro::kilometersToMicroLightYears(float km)
269 {
270 return km / ((float) KM_PER_LY * 1e-6f);
271 }
272
kilometersToMicroLightYears(double km)273 double astro::kilometersToMicroLightYears(double km)
274 {
275 return km / (KM_PER_LY * 1e-6);
276 }
277
microLightYearsToAU(float ly)278 float astro::microLightYearsToAU(float ly)
279 {
280 return ly * (float) AU_PER_LY * 1e-6f;
281 }
282
microLightYearsToAU(double ly)283 double astro::microLightYearsToAU(double ly)
284 {
285 return ly * AU_PER_LY * 1e-6;
286 }
287
AUtoMicroLightYears(float au)288 float astro::AUtoMicroLightYears(float au)
289 {
290 return au / ((float) AU_PER_LY * 1e-6f);
291 }
292
AUtoMicroLightYears(double au)293 double astro::AUtoMicroLightYears(double au)
294 {
295 return au / (AU_PER_LY * 1e-6);
296 }
297
298
299 // Convert the position in univeral coordinates to a star-centric
300 // coordinates in units of kilometers. Note that there are three different
301 // precisions used here: star coordinates are stored as floats in units of
302 // light years, position within a solar system are doubles in units of
303 // kilometers, and p is highest-precision in units of light years.
heliocentricPosition(const UniversalCoord & universal,const Point3f & starPosition)304 Point3d astro::heliocentricPosition(const UniversalCoord& universal,
305 const Point3f& starPosition)
306 {
307 // Get the offset vector
308 Vec3d v = universal - Point3d(starPosition.x * 1e6,
309 starPosition.y * 1e6,
310 starPosition.z * 1e6);
311
312 // . . . and convert it to kilometers
313 return Point3d(microLightYearsToKilometers(v.x),
314 microLightYearsToKilometers(v.y),
315 microLightYearsToKilometers(v.z));
316 }
317
318 // universalPosition is the inverse operation of heliocentricPosition
universalPosition(const Point3d & heliocentric,const Point3f & starPosition)319 UniversalCoord astro::universalPosition(const Point3d& heliocentric,
320 const Point3f& starPosition)
321 {
322 return UniversalCoord(Point3d(starPosition.x * 1e6,
323 starPosition.y * 1e6,
324 starPosition.z * 1e6)) +
325 Vec3d(kilometersToMicroLightYears(heliocentric.x),
326 kilometersToMicroLightYears(heliocentric.y),
327 kilometersToMicroLightYears(heliocentric.z));
328 }
329
330 // universalPosition is the inverse operation of heliocentricPosition
universalPosition(const Point3d & heliocentric,const UniversalCoord & starPosition)331 UniversalCoord astro::universalPosition(const Point3d& heliocentric,
332 const UniversalCoord& starPosition)
333 {
334 return starPosition +
335 Vec3d(kilometersToMicroLightYears(heliocentric.x),
336 kilometersToMicroLightYears(heliocentric.y),
337 kilometersToMicroLightYears(heliocentric.z));
338 }
339
340
341 // Convert equatorial coordinates to Cartesian celestial (or ecliptical)
342 // coordinates.
equatorialToCelestialCart(float ra,float dec,float distance)343 Point3f astro::equatorialToCelestialCart(float ra, float dec, float distance)
344 {
345 double theta = ra / 24.0 * PI * 2 + PI;
346 double phi = (dec / 90.0 - 1.0) * PI / 2;
347 double x = cos(theta) * sin(phi) * distance;
348 double y = cos(phi) * distance;
349 double z = -sin(theta) * sin(phi) * distance;
350
351 return (Point3f((float) x, (float) y, (float) z) * equatorialToCelestial);
352 }
353
354
355 // Convert equatorial coordinates to Cartesian celestial (or ecliptical)
356 // coordinates.
equatorialToCelestialCart(double ra,double dec,double distance)357 Point3d astro::equatorialToCelestialCart(double ra, double dec, double distance)
358 {
359 double theta = ra / 24.0 * PI * 2 + PI;
360 double phi = (dec / 90.0 - 1.0) * PI / 2;
361 double x = cos(theta) * sin(phi) * distance;
362 double y = cos(phi) * distance;
363 double z = -sin(theta) * sin(phi) * distance;
364
365 return (Point3d(x, y, z) * equatorialToCelestiald);
366 }
367
368
anomaly(double meanAnomaly,double eccentricity,double & trueAnomaly,double & eccentricAnomaly)369 void astro::anomaly(double meanAnomaly, double eccentricity,
370 double& trueAnomaly, double& eccentricAnomaly)
371 {
372 double e, delta, err;
373 double tol = 0.00000001745;
374 int iterations = 20; // limit while() to maximum of 20 iterations.
375
376 e = meanAnomaly - 2*PI * (int) (meanAnomaly / (2*PI));
377 err = 1;
378 while(abs(err) > tol && iterations > 0)
379 {
380 err = e - eccentricity*sin(e) - meanAnomaly;
381 delta = err / (1 - eccentricity * cos(e));
382 e -= delta;
383 iterations--;
384 }
385
386 trueAnomaly = 2*atan(sqrt((1+eccentricity)/(1-eccentricity))*tan(e/2));
387 eccentricAnomaly = e;
388 }
389
390
391 /*! Return the angle between the mean ecliptic plane and mean equator at
392 * the specified Julian date.
393 */
394 // TODO: replace this with a better precession model
meanEclipticObliquity(double jd)395 double astro::meanEclipticObliquity(double jd)
396 {
397 double t, de;
398
399 jd -= 2451545.0;
400 t = jd / 36525;
401 de = (46.815 * t + 0.0006 * t * t - 0.00181 * t * t * t) / 3600;
402
403 return J2000Obliquity - de;
404 }
405
406
407 /*! Return a quaternion giving the transformation from the J2000 ecliptic
408 * coordinate system to the J2000 Earth equatorial coordinate system.
409 */
eclipticToEquatorial()410 Quatd astro::eclipticToEquatorial()
411 {
412 return ECLIPTIC_TO_EQUATORIAL_ROTATION;
413 }
414
415
416 /*! Rotate a vector in the J2000 ecliptic coordinate system to
417 * the J2000 Earth equatorial coordinate system.
418 */
eclipticToEquatorial(const Vec3d & v)419 Vec3d astro::eclipticToEquatorial(const Vec3d& v)
420 {
421 return v * ECLIPTIC_TO_EQUATORIAL_MATRIX;
422 }
423
424
425 /*! Return a quaternion giving the transformation from the J2000 Earth
426 * equatorial coordinate system to the galactic coordinate system.
427 */
equatorialToGalactic()428 Quatd astro::equatorialToGalactic()
429 {
430 return EQUATORIAL_TO_GALACTIC_ROTATION;
431 }
432
433
434 /*! Rotate a vector int the J2000 Earth equatorial coordinate system to
435 * the galactic coordinate system.
436 */
equatorialToGalactic(const Vec3d & v)437 Vec3d astro::equatorialToGalactic(const Vec3d& v)
438 {
439 return v * EQUATORIAL_TO_GALACTIC_MATRIX;
440 }
441
442
443
Date()444 astro::Date::Date()
445 {
446 year = 0;
447 month = 0;
448 day = 0;
449 hour = 0;
450 minute = 0;
451 seconds = 0.0;
452 wday = 0;
453 utc_offset = 0;
454 tzname = "UTC";
455 }
456
Date(int Y,int M,int D)457 astro::Date::Date(int Y, int M, int D)
458 {
459 year = Y;
460 month = M;
461 day = D;
462 hour = 0;
463 minute = 0;
464 seconds = 0.0;
465 wday = 0;
466 utc_offset = 0;
467 tzname = "UTC";
468 }
469
Date(double jd)470 astro::Date::Date(double jd)
471 {
472 int64 a = (int64) floor(jd + 0.5);
473 wday = (a + 1) % 7;
474 double c;
475 if (a < 2299161)
476 {
477 c = (double) (a + 1524);
478 }
479 else
480 {
481 double b = (double) ((int64) floor((a - 1867216.25) / 36524.25));
482 c = a + b - (int64) floor(b / 4) + 1525;
483 }
484
485 int64 d = (int64) floor((c - 122.1) / 365.25);
486 int64 e = (int64) floor(365.25 * d);
487 int64 f = (int64) floor((c - e) / 30.6001);
488
489 double dday = c - e - (int64) floor(30.6001 * f) + ((jd + 0.5) - a);
490
491 // This following used to be 14.0, but gcc was computing it incorrectly, so
492 // it was changed to 14
493 month = (int) (f - 1 - 12 * (int64) (f / 14));
494 year = (int) (d - 4715 - (int64) ((7.0 + month) / 10.0));
495 day = (int) dday;
496
497 double dhour = (dday - day) * 24;
498 hour = (int) dhour;
499
500 double dminute = (dhour - hour) * 60;
501 minute = (int) dminute;
502
503 seconds = (dminute - minute) * 60;
504 utc_offset = 0;
505 tzname = "UTC";
506 }
507
toCStr(Format format) const508 const char* astro::Date::toCStr(Format format) const
509 {
510 static char date[255];
511
512 // MinGW's libraries don't have the tm_gmtoff and tm_zone fields for
513 // struct tm.
514 #if defined(__GNUC__) && !defined(_WIN32)
515 struct tm cal_time;
516 memset(&cal_time, 0, sizeof(cal_time));
517 cal_time.tm_year = year-1900;
518 cal_time.tm_mon = month-1;
519 cal_time.tm_mday = day;
520 cal_time.tm_hour = hour;
521 cal_time.tm_min = minute;
522 cal_time.tm_sec = (int)seconds;
523 cal_time.tm_wday = wday;
524 cal_time.tm_gmtoff = utc_offset;
525 #if defined(TARGET_OS_MAC) || defined(__FreeBSD__)
526 // tm_zone is a non-const string field on the Mac and FreeBSD (why?)
527 cal_time.tm_zone = const_cast<char*>(tzname.c_str());
528 #else
529 cal_time.tm_zone = tzname.c_str();
530 #endif
531
532 const char* strftime_format;
533 switch(format)
534 {
535 case Locale:
536 strftime_format = "%c";
537 break;
538 case TZName:
539 strftime_format = "%Y %b %d %H:%M:%S %Z";
540 break;
541 default:
542 strftime_format = "%Y %b %d %H:%M:%S %z";
543 break;
544 }
545
546 strftime(date, sizeof(date), strftime_format, &cal_time);
547 #else
548 switch(format)
549 {
550 case Locale:
551 case TZName:
552 snprintf(date, sizeof(date), "%04d %s %02d %02d:%02d:%02d %s",
553 year, _(MonthAbbrList[month-1]), day,
554 hour, minute, (int)seconds, tzname.c_str());
555 break;
556 case UTCOffset:
557 {
558 int sign = utc_offset < 0 ? -1:1;
559 int h_offset = sign * utc_offset / 3600;
560 int m_offset = (sign * utc_offset - h_offset * 3600) / 60;
561 snprintf(date, sizeof(date), "%04d %s %02d %02d:%02d:%02d %c%02d%02d",
562 year, _(MonthAbbrList[month-1]), day,
563 hour, minute, (int)seconds, (sign==1?'+':'-'), h_offset, m_offset);
564 }
565 break;
566 }
567 #endif
568
569 return date;
570
571 }
572
573 // Convert a calendar date to a Julian date
operator double() const574 astro::Date::operator double() const
575 {
576 int y = year, m = month;
577 if (month <= 2)
578 {
579 y = year - 1;
580 m = month + 12;
581 }
582
583 // Correct for the lost days in Oct 1582 when the Gregorian calendar
584 // replaced the Julian calendar.
585 int B = -2;
586 if (year > 1582 || (year == 1582 && (month > 10 || (month == 10 && day >= 15))))
587 {
588 B = y / 400 - y / 100;
589 }
590
591 return (floor(365.25 * y) +
592 floor(30.6001 * (m + 1)) + B + 1720996.5 +
593 day + hour / 24.0 + minute / 1440.0 + seconds / 86400.0);
594 }
595
596
597 // TODO: need option to parse UTC times (with leap seconds)
parseDate(const string & s,astro::Date & date)598 bool astro::parseDate(const string& s, astro::Date& date)
599 {
600 int year = 0;
601 unsigned int month = 1;
602 unsigned int day = 1;
603 unsigned int hour = 0;
604 unsigned int minute = 0;
605 double second = 0.0;
606
607 if (sscanf(s.c_str(), " %d %u %u %u:%u:%lf ",
608 &year, &month, &day, &hour, &minute, &second) == 6 ||
609 sscanf(s.c_str(), " %d %u %u %u:%u ",
610 &year, &month, &day, &hour, &minute) == 5 ||
611 sscanf(s.c_str(), " %d %u %u ", &year, &month, &day) == 3)
612 {
613 if (month < 1 || month > 12)
614 return false;
615 if (hour > 23 || minute > 59 || second >= 60.0 || second < 0.0)
616 return false;
617
618 // Days / month calculation . . .
619 int maxDay = 31 - ((0xa50 >> month) & 0x1);
620 if (month == 2)
621 {
622 // Check for a leap year
623 if (year % 4 == 0 && (year % 100 != 0 || year % 400 == 0))
624 maxDay = 29;
625 else
626 maxDay = 28;
627 }
628 if (day > (unsigned int) maxDay || day < 1)
629 return false;
630
631 date.year = year;
632 date.month = month;
633 date.day = day;
634 date.hour = hour;
635 date.minute = minute;
636 date.seconds = second;
637
638 return true;
639 }
640
641 return false;
642 }
643
644
645 astro::Date
systemDate()646 astro::Date::systemDate()
647 {
648 time_t t = time(NULL);
649 struct tm *gmt = gmtime(&t);
650 astro::Date d;
651 d.year = gmt->tm_year + 1900;
652 d.month = gmt->tm_mon + 1;
653 d.day = gmt->tm_mday;
654 d.hour = gmt->tm_hour;
655 d.minute = gmt->tm_min;
656 d.seconds = (int) gmt->tm_sec;
657
658 return d;
659 }
660
661
operator <<(ostream & s,const astro::Date d)662 ostream& operator<<(ostream& s, const astro::Date d)
663 {
664 s << d.toCStr();
665 return s;
666 }
667
668
669 /********* Time scale conversion functions ***********/
670
671 // Convert from Atomic Time to UTC
672 astro::Date
TAItoUTC(double tai)673 astro::TAItoUTC(double tai)
674 {
675 unsigned int nRecords = sizeof(LeapSeconds) / sizeof(LeapSeconds[0]);
676 double dAT = LeapSeconds[0].seconds;
677 /*double dD = 0.0; Unused*/
678 int extraSecs = 0;
679
680 for (unsigned int i = nRecords - 1; i > 0; i--)
681 {
682 if (tai - secsToDays(LeapSeconds[i].seconds) >= LeapSeconds[i].t)
683 {
684 dAT = LeapSeconds[i].seconds;
685 break;
686 }
687 else if (tai - secsToDays(LeapSeconds[i - 1].seconds) >= LeapSeconds[i].t)
688 {
689 dAT = LeapSeconds[i].seconds;
690 extraSecs = LeapSeconds[i].seconds - LeapSeconds[i - 1].seconds;
691 break;
692 }
693 }
694
695 Date utcDate(tai - secsToDays(dAT));
696 utcDate.seconds += extraSecs;
697
698 return utcDate;
699 }
700
701
702 // Convert from UTC to Atomic Time
703 double
UTCtoTAI(const astro::Date & utc)704 astro::UTCtoTAI(const astro::Date& utc)
705 {
706 unsigned int nRecords = sizeof(LeapSeconds) / sizeof(LeapSeconds[0]);
707 double dAT = LeapSeconds[0].seconds;
708 double utcjd = (double) Date(utc.year, utc.month, utc.day);
709
710 for (unsigned int i = nRecords - 1; i > 0; i--)
711 {
712 if (utcjd >= LeapSeconds[i].t)
713 {
714 dAT = LeapSeconds[i].seconds;
715 break;
716 }
717 }
718
719 double tai = utcjd + secsToDays(utc.hour * 3600.0 + utc.minute * 60.0 + utc.seconds + dAT);
720
721 return tai;
722 }
723
724
725 // Convert from Terrestrial Time to Atomic Time
726 double
TTtoTAI(double tt)727 astro::TTtoTAI(double tt)
728 {
729 return tt - secsToDays(dTA);
730 }
731
732
733 // Convert from Atomic Time to Terrestrial TIme
734 double
TAItoTT(double tai)735 astro::TAItoTT(double tai)
736 {
737 return tai + secsToDays(dTA);
738 }
739
740
741 // Correction for converting from Terrestrial Time to Barycentric Dynamical
742 // Time. Constants and algorithm from "Time Routines in CSPICE",
743 // http://sohowww.nascom.nasa.gov/solarsoft/stereo/gen/exe/icy/doc/time.req
744 static const double K = 1.657e-3;
745 static const double EB = 1.671e-2;
746 static const double M0 = 6.239996;
747 static const double M1 = 1.99096871e-7;
748
749 // Input is a TDB Julian Date; result is in seconds
TDBcorrection(double tdb)750 double TDBcorrection(double tdb)
751 {
752 // t is seconds from J2000.0
753 double t = astro::daysToSecs(tdb - astro::J2000);
754
755 // Approximate calculation of Earth's mean anomaly
756 double M = M0 + M1 * t;
757
758 // Compute the eccentric anomaly
759 double E = M + EB * sin(M);
760
761 return K * sin(E);
762 }
763
764
765 // Convert from Terrestrial Time to Barycentric Dynamical Time
766 double
TTtoTDB(double tt)767 astro::TTtoTDB(double tt)
768 {
769 return tt + secsToDays(TDBcorrection(tt));
770 }
771
772
773 // Convert from Barycentric Dynamical Time to Terrestrial Time
774 double
TDBtoTT(double tdb)775 astro::TDBtoTT(double tdb)
776 {
777 return tdb - secsToDays(TDBcorrection(tdb));
778 }
779
780
781 // Convert from Coordinated Universal time to Barycentric Dynamical Time
782 astro::Date
TDBtoUTC(double tdb)783 astro::TDBtoUTC(double tdb)
784 {
785 return TAItoUTC(TTtoTAI(TDBtoTT(tdb)));
786 }
787
788 // Convert from Barycentric Dynamical Time to local calendar if possible
789 // otherwise convert to UTC
790 astro::Date
TDBtoLocal(double tdb)791 astro::TDBtoLocal(double tdb)
792 {
793 double tai = astro::TTtoTAI(astro::TDBtoTT(tdb));
794 double jdutc = astro::TAItoJDUTC(tai);
795 if (jdutc < 2465442 &&
796 jdutc > 2415733)
797 {
798 time_t time = (int) astro::julianDateToSeconds(jdutc - 2440587.5);
799 struct tm *localt = localtime(&time);
800 if (localt != NULL)
801 {
802 astro::Date d;
803 d.year = localt->tm_year + 1900;
804 d.month = localt->tm_mon + 1;
805 d.day = localt->tm_mday;
806 d.hour = localt->tm_hour;
807 d.minute = localt->tm_min;
808 d.seconds = (int) localt->tm_sec;
809 d.wday = localt->tm_wday;
810 #if defined(__GNUC__) && !defined(_WIN32)
811 d.utc_offset = localt->tm_gmtoff;
812 d.tzname = localt->tm_zone;
813 #else
814 {
815 astro::Date utcDate = astro::TAItoUTC(tai);
816 int daydiff = d.day - utcDate.day;
817 int hourdiff;
818 if (daydiff == 0)
819 hourdiff = 0;
820 else if (daydiff > 1 || daydiff == -1)
821 hourdiff = -24;
822 else
823 hourdiff = 24;
824 d.utc_offset = (hourdiff + d.hour - utcDate.hour) * 3600
825 + (d.minute - utcDate.minute) * 60;
826 }
827 d.tzname = localt->tm_isdst ? _("DST"): _("STD");
828 #endif
829 return d;
830 }
831 }
832 return astro::TDBtoUTC(tdb);
833 }
834
835 // Convert from Barycentric Dynamical Time to UTC
836 double
UTCtoTDB(const astro::Date & utc)837 astro::UTCtoTDB(const astro::Date& utc)
838 {
839 return TTtoTDB(TAItoTT(UTCtoTAI(utc)));
840 }
841
842
843 // Convert from TAI to Julian Date UTC. The Julian Date UTC functions should
844 // generally be avoided because there's no provision for dealing with leap
845 // seconds.
846 double
JDUTCtoTAI(double utc)847 astro::JDUTCtoTAI(double utc)
848 {
849 unsigned int nRecords = sizeof(LeapSeconds) / sizeof(LeapSeconds[0]);
850 double dAT = LeapSeconds[0].seconds;
851
852 for (unsigned int i = nRecords - 1; i > 0; i--)
853 {
854 if (utc > LeapSeconds[i].t)
855 {
856 dAT = LeapSeconds[i].seconds;
857 break;
858 }
859 }
860
861 return utc + secsToDays(dAT);
862 }
863
864
865 // Convert from Julian Date UTC to TAI
866 double
TAItoJDUTC(double tai)867 astro::TAItoJDUTC(double tai)
868 {
869 unsigned int nRecords = sizeof(LeapSeconds) / sizeof(LeapSeconds[0]);
870 double dAT = LeapSeconds[0].seconds;
871
872 for (unsigned int i = nRecords - 1; i > 0; i--)
873 {
874 if (tai - secsToDays(LeapSeconds[i - 1].seconds) > LeapSeconds[i].t)
875 {
876 dAT = LeapSeconds[i].seconds;
877 break;
878 }
879 }
880
881 return tai - secsToDays(dAT);
882 }
883
884
885