Celestia

1
// date.cpp
2
//
3
// Copyright (C) 2001-2023, the Celestia Development Team
4
// Original version by Chris Laurel <claurel@gmail.com>
5
//
6
// This program is free software; you can redistribute it and/or
7
// modify it under the terms of the GNU General Public License
8
// as published by the Free Software Foundation; either version 2
9
// of the License, or (at your option) any later version.
10

11
#include "date.h"
12

13
#include <array>
14
#include <cmath>
15
#include <cstddef>
16
#include <cstdint>
17
#include <ctime>
18
#include <locale>
19
#include <memory>
20
#include <string_view>
21

22
#include <fmt/format.h>
23

24
#include "astro.h"
25

26
#if defined(__GNUC__) && !defined(_WIN32)
27
#include <fmt/chrono.h>
28
#else
29
#include <iomanip>
30
#include <sstream>
31

32
#include <celutil/gettext.h>
33
#include <celutil/winutil.h>
34
#endif
35

36
using namespace std::string_view_literals;
37

38
namespace celestia::astro
39
{
40

41
namespace
42
{
43

44
// Difference in seconds between Terrestrial Time and International
45
// Atomic Time
46
constexpr double dTA = 32.184;
47

48
// Table of leap second insertions. The leap second always
49
// appears as the last second of the day immediately prior
50
// to the date in the table.
51
constexpr std::array<LeapSecondRecord, 28> LeapSeconds
52
{
53
    LeapSecondRecord{ 10, 2441317.5 }, // 1 Jan 1972
54
    LeapSecondRecord{ 11, 2441499.5 }, // 1 Jul 1972
55
    LeapSecondRecord{ 12, 2441683.5 }, // 1 Jan 1973
56
    LeapSecondRecord{ 13, 2442048.5 }, // 1 Jan 1974
57
    LeapSecondRecord{ 14, 2442413.5 }, // 1 Jan 1975
58
    LeapSecondRecord{ 15, 2442778.5 }, // 1 Jan 1976
59
    LeapSecondRecord{ 16, 2443144.5 }, // 1 Jan 1977
60
    LeapSecondRecord{ 17, 2443509.5 }, // 1 Jan 1978
61
    LeapSecondRecord{ 18, 2443874.5 }, // 1 Jan 1979
62
    LeapSecondRecord{ 19, 2444239.5 }, // 1 Jan 1980
63
    LeapSecondRecord{ 20, 2444786.5 }, // 1 Jul 1981
64
    LeapSecondRecord{ 21, 2445151.5 }, // 1 Jul 1982
65
    LeapSecondRecord{ 22, 2445516.5 }, // 1 Jul 1983
66
    LeapSecondRecord{ 23, 2446247.5 }, // 1 Jul 1985
67
    LeapSecondRecord{ 24, 2447161.5 }, // 1 Jan 1988
68
    LeapSecondRecord{ 25, 2447892.5 }, // 1 Jan 1990
69
    LeapSecondRecord{ 26, 2448257.5 }, // 1 Jan 1991
70
    LeapSecondRecord{ 27, 2448804.5 }, // 1 Jul 1992
71
    LeapSecondRecord{ 28, 2449169.5 }, // 1 Jul 1993
72
    LeapSecondRecord{ 29, 2449534.5 }, // 1 Jul 1994
73
    LeapSecondRecord{ 30, 2450083.5 }, // 1 Jan 1996
74
    LeapSecondRecord{ 31, 2450630.5 }, // 1 Jul 1997
75
    LeapSecondRecord{ 32, 2451179.5 }, // 1 Jan 1999
76
    LeapSecondRecord{ 33, 2453736.5 }, // 1 Jan 2006
77
    LeapSecondRecord{ 34, 2454832.5 }, // 1 Jan 2009
78
    LeapSecondRecord{ 35, 2456109.5 }, // 1 Jul 2012
79
    LeapSecondRecord{ 36, 2457204.5 }, // 1 Jul 2015
80
    LeapSecondRecord{ 37, 2457754.5 }, // 1 Jan 2017
81
};
82

83
celestia::util::array_view<LeapSecondRecord> g_leapSeconds = LeapSeconds; //NOSONAR
84

85
#if !(defined(__GNUC__) && !defined(_WIN32))
86
class MonthAbbreviations
87
{
88
public:
89
    explicit MonthAbbreviations(const std::locale& loc);
90
    std::string_view operator[](int i) const { return abbreviations[static_cast<std::size_t>(i)]; }
91

92
private:
93
    std::array<std::string, 12> abbreviations;
94
};
95

96
MonthAbbreviations::MonthAbbreviations(const std::locale& loc)
97
{
98
    for (int i = 0; i < 12; ++i)
99
    {
100
        std::tm tm;
101
        tm.tm_mon = i;
102
        std::wostringstream stream;
103
        stream.imbue(loc);
104
        stream << std::put_time(&tm, L"%b");
105
        abbreviations[i] = util::WideToUTF8(stream.str());
106
    }
107
}
108
#endif
109

110
inline bool
111
timeToLocal(const std::time_t& time, std::tm& localt)
112
{
113
#ifdef _WIN32
114
    return localtime_s(&localt, &time) == 0;
115
#else
116
    return localtime_r(&time, &localt) != nullptr;
117
#endif
118
}
119

120
inline bool
121
timeToUTC(const std::time_t& time, std::tm& utct)
122
{
123
#ifdef _WIN32
124
    return gmtime_s(&utct, &time) == 0;
125
#else
126
    return gmtime_r(&time, &utct) != nullptr;
127
#endif
128
}
129

130
} // end unnamed namespace
131

132
void
133
setLeapSeconds(celestia::util::array_view<LeapSecondRecord> leapSeconds)
134
{
135
    g_leapSeconds = leapSeconds;
136
}
137

138
Date::Date(int Y, int M, int D) :
139
    year(Y),
140
    month(M),
141
    day(D)
142
{
143
}
144

145
Date::Date(double jd)
146
{
147
    auto a = (std::int64_t) std::floor(jd + 0.5);
148
    wday = (a + 1) % 7;
149
    double c;
150
    if (a < 2299161)
151
    {
152
        c = (double) (a + 1524);
153
    }
154
    else
155
    {
156
        auto b = (double) ((std::int64_t) std::floor((a - 1867216.25) / 36524.25)); //NOSONAR
157
        c = a + b - (std::int64_t) std::floor(b / 4) + 1525; //NOSONAR
158
    }
159

160
    auto d = (std::int64_t) std::floor((c - 122.1) / 365.25);
161
    auto e = (std::int64_t) std::floor(365.25 * d); //NOSONAR
162
    auto f = (std::int64_t) std::floor((c - e) / 30.6001); //NOSONAR
163

164
    double dday = c - e - (std::int64_t) std::floor(30.6001 * f) + ((jd + 0.5) - a); //NOSONAR
165

166
    // This following used to be 14.0, but gcc was computing it incorrectly, so
167
    // it was changed to 14
168
    month = (int) (f - 1 - 12 * (std::int64_t) (f / 14)); //NOSONAR
169
    year = (int) (d - 4715 - (std::int64_t) ((7.0 + month) / 10.0));
170
    day = (int) dday;
171

172
    double dhour = (dday - day) * 24;
173
    hour = (int) dhour;
174

175
    double dminute = (dhour - hour) * 60;
176
    minute = (int) dminute;
177

178
    seconds = (dminute - minute) * 60;
179
    utc_offset = 0;
180
    tzname = "UTC";
181
}
182

183
std::string
184
Date::toString(const std::locale& loc, Format format) const
185
{
186
    if (format == ISO8601)
187
    {
188
        return fmt::format("{:04}-{:02}-{:02}T{:02}:{:02}:{:08.5f}Z"sv,
189
                           year, month, day, hour, minute, seconds);
190
    }
191

192
    // MinGW's libraries don't have the tm_gmtoff and tm_zone fields for
193
    // struct tm.
194
#if defined(__GNUC__) && !defined(_WIN32)
195
    struct tm cal_time {};
196
    cal_time.tm_year = year-1900;
197
    cal_time.tm_mon = month-1;
198
    cal_time.tm_mday = day;
199
    cal_time.tm_hour = hour;
200
    cal_time.tm_min = minute;
201
    cal_time.tm_sec = (int)seconds;
202
    cal_time.tm_wday = wday;
203
    cal_time.tm_gmtoff = utc_offset;
204
#if defined(__APPLE__) || defined(__FreeBSD__)
205
    // tm_zone is a non-const string field on the Mac and FreeBSD (why?)
206
    cal_time.tm_zone = const_cast<char*>(tzname.c_str());
207
#else
208
    cal_time.tm_zone = tzname.c_str();
209
#endif
210

211
    switch(format)
212
    {
213
    case Locale:
214
        return fmt::format(loc, "{:%c}"sv, cal_time);
215
    case TZName:
216
        return fmt::format(loc, "{:%Y %b %d %H:%M:%S %Z}"sv, cal_time);
217
    default:
218
        return fmt::format(loc, "{:%Y %b %d %H:%M:%S %z}"sv, cal_time);
219
    }
220
#else
221
    static const MonthAbbreviations* monthAbbreviations = nullptr;
222
    if (monthAbbreviations == nullptr)
223
        monthAbbreviations = std::make_unique<MonthAbbreviations>(loc).release(); //NOSONAR
224
    switch(format)
225
    {
226
    case Locale:
227
    case TZName:
228
        return fmt::format("{:04} {} {:02} {:02}:{:02}:{:02} {}"sv,
229
                           year, (*monthAbbreviations)[month - 1], day,
230
                           hour, minute, static_cast<int>(seconds), tzname);
231
    default:
232
        {
233
            char sign = utc_offset < 0 ? '-' : '+';
234
            auto offsets = std::div(std::abs(utc_offset), 3600);
235
            return fmt::format("{:04} {} {:02} {:02}:{:02}:{:02} {}{:02}{:02}"sv,
236
                               year, (*monthAbbreviations)[month - 1], day,
237
                               hour, minute, static_cast<int>(seconds),
238
                               sign, offsets.quot, offsets.rem / 60);
239
        }
240
    }
241
#endif
242
}
243

244
// Convert a calendar date to a Julian date
245
Date::operator double() const
246
{
247
    int y = year, m = month;
248
    if (month <= 2)
249
    {
250
        y = year - 1;
251
        m = month + 12;
252
    }
253

254
    // Correct for the lost days in Oct 1582 when the Gregorian calendar
255
    // replaced the Julian calendar.
256
    int B = -2;
257
    if (year > 1582 || (year == 1582 && (month > 10 || (month == 10 && day >= 15))))
258
    {
259
        B = y / 400 - y / 100;
260
    }
261

262
    return (std::floor(365.25 * y) +
263
            std::floor(30.6001 * (m + 1)) + B + 1720996.5 +
264
            day + hour / HOURS_PER_DAY + minute / MINUTES_PER_DAY + seconds / SECONDS_PER_DAY);
265
}
266

267
Date
268
Date::systemDate()
269
{
270
    time_t t = time(nullptr);
271
    tm gmt;
272

273
    Date d;
274
    if (timeToUTC(t, gmt))
275
    {
276
        d.year = gmt.tm_year + 1900;
277
        d.month = gmt.tm_mon + 1;
278
        d.day = gmt.tm_mday;
279
        d.hour = gmt.tm_hour;
280
        d.minute = gmt.tm_min;
281
        d.seconds = gmt.tm_sec;
282
    }
283

284
    return d;
285
}
286

287
// TODO: need option to parse UTC times (with leap seconds)
288
bool
289
parseDate(const std::string& s, Date& date)
290
{
291
    int year = 0;
292
    unsigned int month = 1;
293
    unsigned int day = 1;
294
    unsigned int hour = 0;
295
    unsigned int minute = 0;
296
    double second = 0.0;
297

298
    if (std::sscanf(s.c_str(), "%d-%u-%uT%u:%u:%lf",
299
                    &year, &month, &day, &hour, &minute, &second) != 6 &&
300
        std::sscanf(s.c_str(), " %d %u %u %u:%u:%lf ",
301
                    &year, &month, &day, &hour, &minute, &second) != 6 &&
302
        std::sscanf(s.c_str(), " %d %u %u %u:%u ",
303
                    &year, &month, &day, &hour, &minute) != 5 &&
304
        std::sscanf(s.c_str(), " %d %u %u ", &year, &month, &day) != 3)
305
    {
306
        return false;
307
    }
308

309
    if (month < 1 || month > 12)
310
        return false;
311
    if (hour > 23 || minute > 59 || second >= 60.0 || second < 0.0)
312
        return false;
313

314
    // Days / month calculation . . .
315
    int maxDay = 31 - ((0xa50 >> month) & 0x1);
316
    if (month == 2)
317
    {
318
        // Check for a leap year
319
        if (year % 4 == 0 && (year % 100 != 0 || year % 400 == 0))
320
            maxDay = 29;
321
        else
322
            maxDay = 28;
323
    }
324
    if (day > (unsigned int) maxDay || day < 1)
325
        return false;
326

327
    date.year = year;
328
    date.month = month;
329
    date.day = day;
330
    date.hour = hour;
331
    date.minute = minute;
332
    date.seconds = second;
333

334
    return true;
335
}
336

337
/********* Time scale conversion functions ***********/
338

339
// Convert from Atomic Time to UTC
340
Date
341
TAItoUTC(double tai)
342
{
343
    int dAT = g_leapSeconds[0].seconds;
344
    int extraSecs = 0;
345

346
    for (std::size_t i = g_leapSeconds.size() - 1; i > 0; i--) //NOSONAR
347
    {
348
        if (tai - secsToDays(g_leapSeconds[i].seconds) >= g_leapSeconds[i].t)
349
        {
350
            dAT = g_leapSeconds[i].seconds;
351
            break;
352
        }
353
        if (tai - secsToDays(g_leapSeconds[i - 1].seconds) >= g_leapSeconds[i].t)
354
        {
355
            dAT = g_leapSeconds[i].seconds;
356
            extraSecs = g_leapSeconds[i].seconds - g_leapSeconds[i - 1].seconds;
357
            break;
358
        }
359
    }
360

361
    Date utcDate(tai - secsToDays(dAT));
362
    utcDate.seconds += extraSecs;
363

364
    return utcDate;
365
}
366

367
// Convert from UTC to Atomic Time
368
double
369
UTCtoTAI(const Date& utc)
370
{
371
    double dAT = g_leapSeconds[0].seconds;
372
    auto utcjd = (double) Date(utc.year, utc.month, utc.day);
373

374
    for (std::size_t i = g_leapSeconds.size() - 1; i > 0; i--)
375
    {
376
        if (utcjd >= g_leapSeconds[i].t)
377
        {
378
            dAT = g_leapSeconds[i].seconds;
379
            break;
380
        }
381
    }
382

383
    double tai = utcjd + secsToDays(utc.hour * 3600.0 + utc.minute * 60.0 + utc.seconds + dAT);
384

385
    return tai;
386
}
387

388
// Convert from Terrestrial Time to Atomic Time
389
double
390
TTtoTAI(double tt)
391
{
392
    return tt - secsToDays(dTA);
393
}
394

395
// Convert from Atomic Time to Terrestrial TIme
396
double
397
TAItoTT(double tai)
398
{
399
    return tai + secsToDays(dTA);
400
}
401

402
// Input is a TDB Julian Date; result is in seconds
403
double
404
TDBcorrection(double tdb)
405
{
406
    // Correction for converting from Terrestrial Time to Barycentric Dynamical
407
    // Time. Constants and algorithm from "Time Routines in CSPICE",
408
    // http://sohowww.nascom.nasa.gov/solarsoft/stereo/gen/exe/icy/doc/time.req
409
    constexpr double K  = 1.657e-3;
410
    constexpr double EB = 1.671e-2;
411
    constexpr double M0 = 6.239996;
412
    constexpr double M1 = 1.99096871e-7;
413

414
    // t is seconds from J2000.0
415
    double t = daysToSecs(tdb - J2000);
416

417
    // Approximate calculation of Earth's mean anomaly
418
    double M = M0 + M1 * t;
419

420
    // Compute the eccentric anomaly
421
    double E = M + EB * std::sin(M);
422

423
    return K * std::sin(E);
424
}
425

426
// Convert from Terrestrial Time to Barycentric Dynamical Time
427
double
428
TTtoTDB(double tt)
429
{
430
    return tt + secsToDays(TDBcorrection(tt));
431
}
432

433
// Convert from Barycentric Dynamical Time to Terrestrial Time
434
double
435
TDBtoTT(double tdb)
436
{
437
    return tdb - secsToDays(TDBcorrection(tdb));
438
}
439

440
// Convert from Coordinated Universal time to Barycentric Dynamical Time
441
Date
442
TDBtoUTC(double tdb)
443
{
444
    return TAItoUTC(TTtoTAI(TDBtoTT(tdb)));
445
}
446

447
// Convert from Barycentric Dynamical Time to local calendar if possible
448
// otherwise convert to UTC
449
Date
450
TDBtoLocal(double tdb)
451
{
452
    double tai = TTtoTAI(TDBtoTT(tdb));
453

454
    double jdutc = TAItoJDUTC(tai);
455
    if (jdutc <= 2415733 || jdutc >= 2465442)
456
        return TDBtoUTC(tdb);
457

458
    auto time = static_cast<time_t>(julianDateToSeconds(jdutc - 2440587.5));
459
    struct tm localt;
460

461
    if (!timeToLocal(time, localt)) //NOSONAR
462
        return TDBtoUTC(tdb);
463

464
    Date d;
465
    d.year = localt.tm_year + 1900;
466
    d.month = localt.tm_mon + 1;
467
    d.day = localt.tm_mday;
468
    d.hour = localt.tm_hour;
469
    d.minute = localt.tm_min;
470
    d.seconds = localt.tm_sec;
471
    d.wday = localt.tm_wday;
472
#if defined(__GNUC__) && !defined(_WIN32)
473
    d.utc_offset = static_cast<int>(localt.tm_gmtoff);
474
    d.tzname = localt.tm_zone;
475
#else
476
    Date utcDate = TAItoUTC(tai);
477
    int daydiff = d.day - utcDate.day;
478
    int hourdiff;
479
    if (daydiff == 0)
480
        hourdiff = 0;
481
    else if (daydiff > 1 || daydiff == -1)
482
        hourdiff = -24;
483
    else
484
        hourdiff = 24;
485
    d.utc_offset = (hourdiff + d.hour - utcDate.hour) * 3600
486
                    + (d.minute - utcDate.minute) * 60;
487
    d.tzname = localt.tm_isdst ? _("DST"): _("STD");
488
#endif
489
    return d;
490
}
491

492
// Convert from Barycentric Dynamical Time to UTC
493
double
494
UTCtoTDB(const Date& utc)
495
{
496
    return TTtoTDB(TAItoTT(UTCtoTAI(utc)));
497
}
498

499
// Convert from TAI to Julian Date UTC. The Julian Date UTC functions should
500
// generally be avoided because there's no provision for dealing with leap
501
// seconds.
502
double
503
JDUTCtoTAI(double utc)
504
{
505
    double dAT = g_leapSeconds[0].seconds;
506

507
    for (std::size_t i = g_leapSeconds.size() - 1; i > 0; i--)
508
    {
509
        if (utc > g_leapSeconds[i].t)
510
        {
511
            dAT = g_leapSeconds[i].seconds;
512
            break;
513
        }
514
    }
515

516
    return utc + secsToDays(dAT);
517
}
518

519
// Convert from Julian Date UTC to TAI
520
double
521
TAItoJDUTC(double tai)
522
{
523
    double dAT = g_leapSeconds[0].seconds;
524

525
    for (std::size_t i = g_leapSeconds.size() - 1; i > 0; i--)
526
    {
527
        if (tai - secsToDays(g_leapSeconds[i - 1].seconds) > g_leapSeconds[i].t)
528
        {
529
            dAT = g_leapSeconds[i].seconds;
530
            break;
531
        }
532
    }
533

534
    return tai - secsToDays(dAT);
535
}
536

537
} // end namespace celestia::astro
538