InsTime.hpp

Go to the documentation of this file.

// This file is part of INSTINCT, the INS Toolkit for Integrated
// Navigation Concepts and Training by the Institute of Navigation of
// the University of Stuttgart, Germany.
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at https://mozilla.org/MPL/2.0/.

 
#pragma once
 
#include <cstdint>
#include <string>
#include <array>
#include <limits>
#include <iostream>
#include <chrono>
 
#include <fmt/ostream.h>
#include "gcem.hpp"
 
#include <nlohmann/json.hpp>
using json = nlohmann::json; 
 
#include "TimeSystem.hpp"
#include "Navigation/Math/Math.hpp"
 
namespace NAV
{
namespace InsTimeUtil
{
constexpr int32_t END_OF_THE_CENTURY_MJD = 400000;   
constexpr int32_t WEEKS_PER_GPS_CYCLE = 1024;        
constexpr int32_t DIFF_TO_6_1_1980_MJD = 44244;      
constexpr int32_t DIFF_TO_1_1_1970_MJD = 40587;      
constexpr int32_t DIFF_BDT_WEEK_TO_GPST_WEEK = 1356; 
 
constexpr int32_t DIFF_MJD_TO_JD_DAYS = 2400000;  
constexpr long double DIFF_MJD_TO_JD_FRAC = 0.5L; 
 
constexpr int32_t MONTHS_PER_YEAR = 12;                                                    
constexpr int32_t DAYS_PER_YEAR = 365;                                                     
constexpr int32_t DAYS_PER_WEEK = 7;                                                       
constexpr int32_t HOURS_PER_DAY = 24;                                                      
constexpr int32_t HOURS_PER_WEEK = HOURS_PER_DAY * DAYS_PER_WEEK;                          
constexpr int32_t MINUTES_PER_HOUR = 60;                                                   
constexpr int32_t MINUTES_PER_DAY = MINUTES_PER_HOUR * HOURS_PER_DAY;                      
constexpr int32_t MINUTES_PER_WEEK = MINUTES_PER_DAY * DAYS_PER_WEEK;                      
constexpr int32_t SECONDS_PER_MINUTE = 60;                                                 
constexpr int32_t SECONDS_PER_HOUR = SECONDS_PER_MINUTE * MINUTES_PER_HOUR;                
constexpr int32_t SECONDS_PER_DAY = SECONDS_PER_MINUTE * MINUTES_PER_HOUR * HOURS_PER_DAY; 
constexpr int32_t SECONDS_PER_WEEK = SECONDS_PER_DAY * DAYS_PER_WEEK;                      

constexpr long double EPSILON = std::max(1e-17L, 10 * std::numeric_limits<long double>::epsilon());

constexpr std::array<int32_t, 20> GPS_LEAP_SEC_MJD = {
    0,     // + 0 at 1 Jan 1980 and before
    44786, // + 1 at 1 Jul 1981 = diff UTC-TAI: 20
    45151, // + 2 at 1 Jul 1982 = diff UTC-TAI: 21
    45516, // + 3 at 1 Jul 1983 = diff UTC-TAI: 22
    46247, // + 4 at 1 Jul 1985 = diff UTC-TAI: 23
    47161, // + 5 at 1 Jan 1988 = diff UTC-TAI: 24
    47892, // + 6 at 1 Jan 1990 = diff UTC-TAI: 25
    48257, // + 7 at 1 Jan 1991 = diff UTC-TAI: 26
    48804, // + 8 at 1 Jul 1992 = diff UTC-TAI: 27
    49169, // + 9 at 1 Jul 1993 = diff UTC-TAI: 28
    49534, // +10 at 1 Jul 1994 = diff UTC-TAI: 29
    50083, // +11 at 1 Jan 1996 = diff UTC-TAI: 30
    50630, // +12 at 1 Jul 1997 = diff UTC-TAI: 31
    51179, // +13 at 1 Jan 1999 = diff UTC-TAI: 32
    53736, // +14 at 1 Jan 2006 = diff UTC-TAI: 33
    54832, // +15 at 1 Jan 2009 = diff UTC-TAI: 34
    56109, // +16 at 1 Jul 2012 = diff UTC-TAI: 35
    57204, // +17 at 1 Jul 2015 = diff UTC-TAI: 36
    57754, // +18 at 1 Jan 2017 = diff UTC-TAI: 37
    99999, // future
};

constexpr bool isLeapYear(int32_t year)
{
    return (((year % 4 == 0) && (year % 100 != 0)) || (year % 400 == 0));
}

constexpr int32_t daysInMonth(int32_t month, int32_t year)
{
    --month; // Make month 0 based
    if (month >= InsTimeUtil::MONTHS_PER_YEAR)
    {
        year += month / InsTimeUtil::MONTHS_PER_YEAR;
        month %= InsTimeUtil::MONTHS_PER_YEAR;
    }
    while (month < 0)
    {
        month += MONTHS_PER_YEAR;
        year--;
    }
    ++month; // Make month 1 based
 
    switch (month)
    {
    case 1:
        return 31;
    case 2:
        if (isLeapYear(year))
        {
            return 29;
        }
        return 28;
    case 3:
        return 31;
    case 4:
        return 30;
    case 5:
        return 31;
    case 6:
        return 30;
    case 7:
        return 31;
    case 8:
        return 31;
    case 9:
        return 30;
    case 10:
        return 31;
    case 11:
        return 30;
    case 12:
        return 31;
    default:
        return 366;
    }
}

constexpr std::string formatDuration(double seconds, size_t digits = 9)
{
    std::string str;
    if (seconds >= SECONDS_PER_DAY)
    {
        auto days = static_cast<int32_t>(seconds) / SECONDS_PER_DAY;
        str += fmt::format("{}{:d}d", str.empty() ? "" : " ", days);
        seconds -= days * SECONDS_PER_DAY;
    }
    if (seconds >= SECONDS_PER_HOUR)
    {
        auto hours = static_cast<int32_t>(seconds) / SECONDS_PER_HOUR;
        str += fmt::format("{}{:d}h", str.empty() ? "" : " ", hours);
        seconds -= hours * SECONDS_PER_HOUR;
    }
    if (seconds >= SECONDS_PER_MINUTE)
    {
        auto minutes = static_cast<int32_t>(seconds) / SECONDS_PER_MINUTE;
        str += fmt::format("{}{:d}m", str.empty() ? "" : " ", minutes);
        seconds -= minutes * SECONDS_PER_MINUTE;
    }
    if (seconds > 1e-9)
    {
        size_t d = 0;
        for (; d < digits + 1; d++)
        {
            if (2.0 * std::abs(std::round(seconds * std::pow(10.0, d)) - seconds * std::pow(10.0, d))
                < std::pow(10.0, -(static_cast<double>(digits) - static_cast<double>(d))))
            {
                break;
            }
        }
        str += fmt::format("{}{:.{}f}s", str.empty() ? "" : " ", seconds, d);
    }
    return str;
}
 
} // namespace InsTimeUtil

struct InsTime_MJD
{
    int32_t mjd_day = 0;         
    long double mjd_frac = 0.0L; 

    constexpr InsTime_MJD() = default;

    constexpr InsTime_MJD(int32_t mjd_day, long double mjd_frac)
        : mjd_day(mjd_day), mjd_frac(mjd_frac)
    {
        if (this->mjd_frac >= 1.0L)
        {
            this->mjd_day += static_cast<int32_t>(this->mjd_frac);
            this->mjd_frac -= static_cast<int32_t>(this->mjd_frac);
        }
        while (this->mjd_frac < 0.0L)
        {
            this->mjd_frac += 1.0L;
            this->mjd_day--;
        }
    }

    constexpr bool operator==(const InsTime_MJD& rhs) const
    {
        if (mjd_day == rhs.mjd_day)
        {
            auto difference = gcem::abs(mjd_frac - rhs.mjd_frac);
            return difference <= InsTimeUtil::EPSILON;
        }
        if (auto diffDays = mjd_day - rhs.mjd_day;
            gcem::abs(diffDays) == 1)
        {
            auto difference = 1 + diffDays * (mjd_frac - rhs.mjd_frac);
            return difference <= InsTimeUtil::EPSILON;
        }
        return false;
    }

    constexpr bool operator!=(const InsTime_MJD& rhs) const
    {
        return !(*this == rhs);
    }

    constexpr bool operator<=(const InsTime_MJD& rhs) const
    {
        return *this < rhs || *this == rhs;
    }

    constexpr bool operator>=(const InsTime_MJD& rhs) const
    {
        return *this > rhs || *this == rhs;
    }

    constexpr bool operator<(const InsTime_MJD& rhs) const
    {
        return (mjd_day < rhs.mjd_day || (mjd_day == rhs.mjd_day && mjd_frac < rhs.mjd_frac))
               && *this != rhs;
    }

    constexpr bool operator>(const InsTime_MJD& rhs) const
    {
        return !(*this <= rhs);
    }

    explicit operator std::string() const;
};

struct InsTime_JD
{
    int32_t jd_day{};      
    long double jd_frac{}; 

    constexpr InsTime_JD(int32_t jd_day, long double jd_frac)
        : jd_day(jd_day), jd_frac(jd_frac)
    {
        if (this->jd_frac >= 1.0L)
        {
            this->jd_day += static_cast<int32_t>(this->jd_frac);
            this->jd_frac -= static_cast<int32_t>(this->jd_frac);
        }
        while (this->jd_frac < 0.0L)
        {
            this->jd_frac += 1.0L;
            this->jd_day--;
        }
    }

    constexpr bool operator==(const InsTime_JD& rhs) const
    {
        if (jd_day == rhs.jd_day)
        {
            auto difference = gcem::abs(jd_frac - rhs.jd_frac);
            return difference <= InsTimeUtil::EPSILON;
        }
        if (auto diffDays = jd_day - rhs.jd_day;
            gcem::abs(diffDays) == 1)
        {
            auto difference = 1 + diffDays * (jd_frac - rhs.jd_frac);
            return difference <= InsTimeUtil::EPSILON;
        }
        return false;
    }

    constexpr bool operator!=(const InsTime_JD& rhs) const
    {
        return !(*this == rhs);
    }

    constexpr bool operator<=(const InsTime_JD& rhs) const
    {
        return *this < rhs || *this == rhs;
    }

    constexpr bool operator>=(const InsTime_JD& rhs) const
    {
        return *this > rhs || *this == rhs;
    }

    constexpr bool operator<(const InsTime_JD& rhs) const
    {
        return (jd_day < rhs.jd_day || (jd_day == rhs.jd_day && jd_frac < rhs.jd_frac))
               && *this != rhs;
    }

    constexpr bool operator>(const InsTime_JD& rhs) const
    {
        return !(*this <= rhs);
    }

    explicit operator std::string() const;
};

struct InsTime_GPSweekTow
{
    int32_t gpsCycle{}; 
    int32_t gpsWeek{};  
    long double tow{};  

    constexpr InsTime_GPSweekTow(int32_t gpsCycle, int32_t gpsWeek, long double tow)
        : gpsCycle(gpsCycle), gpsWeek(gpsWeek), tow(tow)
    {
        if (this->tow >= InsTimeUtil::SECONDS_PER_WEEK)
        {
            this->gpsWeek += static_cast<int32_t>(this->tow / InsTimeUtil::SECONDS_PER_WEEK);
            this->tow = gcem::fmod(this->tow, InsTimeUtil::SECONDS_PER_WEEK);
        }
        while (this->tow < 0.0L)
        {
            this->tow += InsTimeUtil::SECONDS_PER_WEEK;
            this->gpsWeek--;
        }
 
        if (this->gpsWeek >= InsTimeUtil::WEEKS_PER_GPS_CYCLE)
        {
            this->gpsCycle += this->gpsWeek / InsTimeUtil::WEEKS_PER_GPS_CYCLE;
            this->gpsWeek %= InsTimeUtil::WEEKS_PER_GPS_CYCLE;
        }
        while (this->gpsWeek < 0)
        {
            this->gpsWeek += InsTimeUtil::WEEKS_PER_GPS_CYCLE;
            this->gpsCycle--;
        }
    };

    constexpr bool operator==(const InsTime_GPSweekTow& rhs) const
    {
        if (gpsCycle == rhs.gpsCycle && gpsWeek == rhs.gpsWeek)
        {
            return gcem::abs(tow - rhs.tow) <= InsTimeUtil::EPSILON;
        }
        if (auto diffWeeks = gpsCycle * InsTimeUtil::WEEKS_PER_GPS_CYCLE + gpsWeek - (rhs.gpsCycle * InsTimeUtil::WEEKS_PER_GPS_CYCLE + rhs.gpsWeek);
            gcem::abs(diffWeeks) == 1)
        {
            return gcem::abs(tow + diffWeeks * InsTimeUtil::SECONDS_PER_WEEK - rhs.tow) <= InsTimeUtil::EPSILON;
        }
        return false;
    }

    constexpr bool operator!=(const InsTime_GPSweekTow& rhs) const
    {
        return !(*this == rhs);
    }

    constexpr bool operator<=(const InsTime_GPSweekTow& rhs) const
    {
        return *this < rhs || *this == rhs;
    }

    constexpr bool operator>=(const InsTime_GPSweekTow& rhs) const
    {
        return *this > rhs || *this == rhs;
    }

    constexpr bool operator<(const InsTime_GPSweekTow& rhs) const
    {
        return (gpsCycle < rhs.gpsCycle
                || (gpsCycle == rhs.gpsCycle && gpsWeek < rhs.gpsWeek)
                || (gpsCycle == rhs.gpsCycle && gpsWeek == rhs.gpsWeek && tow < rhs.tow))
               && *this != rhs;
    }

    constexpr bool operator>(const InsTime_GPSweekTow& rhs) const
    {
        return !(*this <= rhs);
    }

    explicit operator std::string() const;
};

struct InsTime_YMDHMS
{
    int32_t year{};    
    int32_t month{};   
    int32_t day{};     
    int32_t hour{};    
    int32_t min{};     
    long double sec{}; 

    constexpr InsTime_YMDHMS(int32_t year, int32_t month, int32_t day, int32_t hour, int32_t min, long double sec, int digits = -1)
        : year(year), month(month), day(day), hour(hour), min(min), sec(sec)
    {
        if (digits >= 0) { this->sec = math::round(this->sec, static_cast<size_t>(digits)); }
        if (this->sec >= InsTimeUtil::SECONDS_PER_MINUTE)
        {
            this->min += static_cast<int32_t>(this->sec / InsTimeUtil::SECONDS_PER_MINUTE);
            this->sec = gcem::fmod(this->sec, InsTimeUtil::SECONDS_PER_MINUTE);
        }
        while (this->sec < 0.0L)
        {
            this->sec += InsTimeUtil::SECONDS_PER_MINUTE;
            this->min--;
        }
        if (digits >= 0) { this->sec = math::round(this->sec, static_cast<size_t>(digits)); }
 
        if (this->min >= InsTimeUtil::MINUTES_PER_HOUR)
        {
            this->hour += this->min / InsTimeUtil::MINUTES_PER_HOUR;
            this->min %= InsTimeUtil::MINUTES_PER_HOUR;
        }
        while (this->min < 0)
        {
            this->min += InsTimeUtil::MINUTES_PER_HOUR;
            this->hour--;
        }
 
        if (this->hour >= InsTimeUtil::HOURS_PER_DAY)
        {
            this->day += this->hour / InsTimeUtil::HOURS_PER_DAY;
            this->hour %= InsTimeUtil::HOURS_PER_DAY;
        }
        while (this->hour < 0)
        {
            this->hour += InsTimeUtil::HOURS_PER_DAY;
            this->day--;
        }
 
        while (this->day >= InsTimeUtil::DAYS_PER_YEAR + static_cast<int32_t>(InsTimeUtil::isLeapYear(this->year)))
        {
            this->day -= InsTimeUtil::DAYS_PER_YEAR + static_cast<int32_t>(InsTimeUtil::isLeapYear(this->year));
            this->year++;
        }
 
        while (this->day > InsTimeUtil::daysInMonth(this->month, this->year))
        {
            this->day -= InsTimeUtil::daysInMonth(this->month, this->year);
            this->month++;
        }
        while (this->day < 1)
        {
            this->month--;
            this->day += InsTimeUtil::daysInMonth(this->month, this->year);
        }
 
        while (this->month > InsTimeUtil::MONTHS_PER_YEAR)
        {
            this->month -= InsTimeUtil::MONTHS_PER_YEAR;
            this->year++;
        }
        while (this->month < 1)
        {
            this->month += InsTimeUtil::MONTHS_PER_YEAR;
            this->year--;
        }
    }

    constexpr bool operator==(const InsTime_YMDHMS& rhs) const
    {
        if (year == rhs.year && month == rhs.month && day == rhs.day && hour == rhs.hour && min == rhs.min)
        {
            return gcem::abs(sec - rhs.sec) <= InsTimeUtil::EPSILON;
        }
        InsTime_YMDHMS this_plus = InsTime_YMDHMS(this->year, this->month, this->day, this->hour, this->min, this->sec + 10);
        InsTime_YMDHMS rhs_plus = InsTime_YMDHMS(rhs.year, rhs.month, rhs.day, rhs.hour, rhs.min, rhs.sec + 10);
        if (this_plus.year == rhs_plus.year && this_plus.month == rhs_plus.month && this_plus.day == rhs_plus.day && this_plus.hour == rhs_plus.hour && this_plus.min == rhs_plus.min)
        {
            return gcem::abs(this_plus.sec - rhs_plus.sec) <= InsTimeUtil::EPSILON;
        }
        return false;
    }

    constexpr bool operator!=(const InsTime_YMDHMS& rhs) const
    {
        return !(*this == rhs);
    }

    constexpr bool operator<=(const InsTime_YMDHMS& rhs) const
    {
        return *this < rhs || *this == rhs;
    }

    constexpr bool operator>=(const InsTime_YMDHMS& rhs) const
    {
        return *this > rhs || *this == rhs;
    }

    constexpr bool operator<(const InsTime_YMDHMS& rhs) const
    {
        return (year < rhs.year
                || (year == rhs.year && month < rhs.month)
                || (year == rhs.year && month == rhs.month && day < rhs.day)
                || (year == rhs.year && month == rhs.month && day == rhs.day && hour < rhs.hour)
                || (year == rhs.year && month == rhs.month && day == rhs.day && hour == rhs.hour && min < rhs.min)
                || (year == rhs.year && month == rhs.month && day == rhs.day && hour == rhs.hour && min == rhs.min && sec < rhs.sec))
               && *this != rhs;
    }

    constexpr bool operator>(const InsTime_YMDHMS& rhs) const
    {
        return !(*this <= rhs);
    }

    explicit operator std::string() const;
};

struct InsTime_YDoySod
{
    int32_t year{};    
    int32_t doy{};     
    long double sod{}; 

    constexpr InsTime_YDoySod(int32_t year, int32_t doy, long double sod)
        : year(year), doy(doy), sod(sod)
    {
        if (this->sod >= InsTimeUtil::SECONDS_PER_DAY)
        {
            this->doy += static_cast<int32_t>(this->sod / InsTimeUtil::SECONDS_PER_DAY);
            this->sod = gcem::fmod(this->sod, InsTimeUtil::SECONDS_PER_DAY);
        }
        while (this->sod < 0)
        {
            this->sod += InsTimeUtil::SECONDS_PER_DAY;
            this->doy--;
        }
 
        while (this->doy > InsTimeUtil::DAYS_PER_YEAR + static_cast<int32_t>(InsTimeUtil::isLeapYear(this->year)))
        {
            this->doy -= InsTimeUtil::DAYS_PER_YEAR + static_cast<int32_t>(InsTimeUtil::isLeapYear(this->year));
            this->year++;
        }
        while (this->doy < 1)
        {
            this->doy += InsTimeUtil::DAYS_PER_YEAR + static_cast<int32_t>(InsTimeUtil::isLeapYear(this->year - 1));
            this->year--;
        }
    }

    constexpr bool operator==(const InsTime_YDoySod& rhs) const
    {
        if (year == rhs.year && doy == rhs.doy)
        {
            return gcem::abs(sod - rhs.sod) <= InsTimeUtil::EPSILON;
        }
        if (auto diffDays = year * InsTimeUtil::DAYS_PER_YEAR + doy - (rhs.year * InsTimeUtil::DAYS_PER_YEAR + rhs.doy);
            gcem::abs(diffDays) == 1)
        {
            auto difference = gcem::abs(sod + diffDays * InsTimeUtil::SECONDS_PER_DAY - rhs.sod);
            return difference <= InsTimeUtil::EPSILON;
        }
        return false;
    }

    constexpr bool operator!=(const InsTime_YDoySod& rhs) const
    {
        return !(*this == rhs);
    }

    constexpr bool operator<=(const InsTime_YDoySod& rhs) const
    {
        return *this < rhs || *this == rhs;
    }

    constexpr bool operator>=(const InsTime_YDoySod& rhs) const
    {
        return *this > rhs || *this == rhs;
    }

    constexpr bool operator<(const InsTime_YDoySod& rhs) const
    {
        return (year < rhs.year
                || (year == rhs.year && doy < rhs.doy)
                || (year == rhs.year && doy == rhs.doy && sod < rhs.sod))
               && *this != rhs;
    }

    constexpr bool operator>(const InsTime_YDoySod& rhs) const
    {
        return !(*this <= rhs);
    }

    explicit operator std::string() const;
};

class InsTime
{
  public:
    /* ------------------------------ Constructors ------------------------------ */

    constexpr InsTime() = default;

    constexpr explicit InsTime(const InsTime_MJD& mjd, TimeSystem timesys = UTC)
        : _mjd(mjd)
    {
        *this -= std::chrono::duration<long double>(differenceToUTC(timesys));
    }

    constexpr explicit InsTime(const InsTime_JD& jd, TimeSystem timesys = UTC)
        : _mjd(jd.jd_day - InsTimeUtil::DIFF_MJD_TO_JD_DAYS, jd.jd_frac - InsTimeUtil::DIFF_MJD_TO_JD_FRAC)
    {
        *this -= std::chrono::duration<long double>(differenceToUTC(timesys));
    }

    constexpr explicit InsTime(const InsTime_GPSweekTow& gpsWeekTow, TimeSystem timesys = GPST)
    {
        auto mjd_day = static_cast<int32_t>((gpsWeekTow.gpsCycle * InsTimeUtil::WEEKS_PER_GPS_CYCLE + gpsWeekTow.gpsWeek) * InsTimeUtil::DAYS_PER_WEEK
                                            + gcem::floor(gpsWeekTow.tow / InsTimeUtil::SECONDS_PER_DAY)
                                            + InsTimeUtil::DIFF_TO_6_1_1980_MJD);
        long double mjd_frac = gcem::fmod(gpsWeekTow.tow, InsTimeUtil::SECONDS_PER_DAY) / InsTimeUtil::SECONDS_PER_DAY;
 
        _mjd = InsTime_MJD(mjd_day, mjd_frac);
 
        *this -= std::chrono::duration<long double>(differenceToUTC(timesys));
    }

    constexpr explicit InsTime(const InsTime_YMDHMS& yearMonthDayHMS, TimeSystem timesys = UTC)
    {
        auto a = static_cast<int32_t>((14 - yearMonthDayHMS.month) / static_cast<double>(InsTimeUtil::MONTHS_PER_YEAR));
        int32_t y = yearMonthDayHMS.year + 4800 - a;
        int32_t m = yearMonthDayHMS.month + InsTimeUtil::MONTHS_PER_YEAR * a - 3;
 
        auto jd_day = static_cast<int32_t>(yearMonthDayHMS.day
                                           + gcem::floor((153.0 * static_cast<double>(m) + 2.0) / 5.0)
                                           + y * InsTimeUtil::DAYS_PER_YEAR
                                           + gcem::floor(static_cast<double>(y) / 4.0)
                                           - gcem::floor(static_cast<double>(y) / 100.0)
                                           + gcem::floor(static_cast<double>(y) / 400.0)
                                           - 32045);
        auto jd_frac = (yearMonthDayHMS.sec
                        + static_cast<long double>(yearMonthDayHMS.min) * InsTimeUtil::SECONDS_PER_MINUTE
                        + static_cast<long double>(yearMonthDayHMS.hour - 12.0) * InsTimeUtil::SECONDS_PER_HOUR)
                       / InsTimeUtil::SECONDS_PER_DAY;
 
        _mjd = InsTime(InsTime_JD(jd_day, jd_frac)).toMJD();
 
        *this -= std::chrono::duration<long double>(differenceToUTC(timesys));
    }

    constexpr explicit InsTime(const InsTime_YDoySod& yearDoySod, TimeSystem timesys = UTC)
    {
        auto year = yearDoySod.year;
        auto doy = yearDoySod.doy;
        auto sod = yearDoySod.sod;
 
        int32_t month = 1;
        while (doy > InsTimeUtil::daysInMonth(month, year))
        {
            doy -= InsTimeUtil::daysInMonth(month, year);
            month++;
        }
 
        _mjd = InsTime(InsTime_YMDHMS(year, month, doy, 0, 0, sod)).toMJD();
 
        *this -= std::chrono::duration<long double>(differenceToUTC(timesys));
    }

    constexpr InsTime(int32_t gpsCycle, int32_t gpsWeek, long double tow, TimeSystem timesys = GPST)
        : InsTime(InsTime_GPSweekTow(gpsCycle, gpsWeek, tow), timesys) {}

    constexpr InsTime(int32_t year, int32_t month, int32_t day, int32_t hour, int32_t min, long double sec, TimeSystem timesys = UTC)
        : InsTime(InsTime_YMDHMS(year, month, day, hour, min, sec), timesys) {}

    ~InsTime() = default;
    constexpr InsTime(const InsTime&) = default;
    constexpr InsTime(InsTime&&) = default;
    constexpr InsTime& operator=(const InsTime&) = default;
    constexpr InsTime& operator=(InsTime&&) = default;
 
    /* ------------------------ Transformation functions ------------------------ */

    [[nodiscard]] constexpr InsTime_MJD toMJD(TimeSystem timesys = UTC) const
    {
        long double mjdFrac = _mjd.mjd_frac + static_cast<long double>(differenceToUTC(timesys)) / static_cast<long double>(InsTimeUtil::SECONDS_PER_DAY);
        return { _mjd.mjd_day, mjdFrac };
    }

    [[nodiscard]] constexpr InsTime_JD toJD(TimeSystem timesys = UTC) const
    {
        auto jd_day = _mjd.mjd_day + InsTimeUtil::DIFF_MJD_TO_JD_DAYS;
        auto jd_frac = _mjd.mjd_frac + InsTimeUtil::DIFF_MJD_TO_JD_FRAC;
 
        jd_frac += static_cast<long double>(differenceToUTC(timesys)) / static_cast<long double>(InsTimeUtil::SECONDS_PER_DAY);
 
        return { jd_day, jd_frac };
    }

    [[nodiscard]] constexpr InsTime_GPSweekTow toGPSweekTow(TimeSystem timesys = GPST) const
    {
        InsTime_MJD mjd_leap = _mjd;
        // Convert from UTC to intended time system
        mjd_leap.mjd_frac += static_cast<long double>(differenceToUTC(timesys)) / static_cast<long double>(InsTimeUtil::SECONDS_PER_DAY);
 
        // Put everything in the time of week, as it gets splitted in InsTime_GPSweekTow constructor
        auto tow = static_cast<long double>((mjd_leap.mjd_day - InsTimeUtil::DIFF_TO_6_1_1980_MJD)) * InsTimeUtil::SECONDS_PER_DAY
                   + mjd_leap.mjd_frac * InsTimeUtil::SECONDS_PER_DAY;
 
        return { 0, 0, tow };
    }

    [[nodiscard]] constexpr InsTime_YMDHMS toYMDHMS(TimeSystem timesys = UTC, int digits = -1) const
    {
        // transform MJD to JD
        InsTime_JD jd = toJD();
        jd.jd_frac = jd.jd_frac + 0.5L;
        jd.jd_frac += static_cast<long double>(differenceToUTC(timesys)) / static_cast<long double>(InsTimeUtil::SECONDS_PER_DAY);
        if (jd.jd_frac >= 1.0L)
        {
            jd.jd_day += static_cast<int32_t>(jd.jd_frac);
            jd.jd_frac -= gcem::floor(jd.jd_frac);
        }
        // transform JD to YMDHMS
        double a = 32044.0 + jd.jd_day;
        double b = gcem::floor((4.0 * a + 3.0) / 146097.0);
        double c = a - gcem::floor((b * 146097.0) / 4.0);
 
        double d = gcem::floor((4.0 * c + 3.0) / 1461.0);
        double e = c - gcem::floor((1461.0 * d) / 4.0);
        double m = gcem::floor((5.0 * e + 2.0) / 153.0);
 
        auto day = static_cast<uint16_t>(e - gcem::floor((153.0 * m + 2.0) / 5.0) + 1);
        auto month = static_cast<uint16_t>(m + 3 - 12 * gcem::floor(m / 10.0));
        auto year = static_cast<uint16_t>(b * 100 + d - 4800.0 + gcem::floor(m / 10.0));
 
        long double sec = jd.jd_frac * InsTimeUtil::SECONDS_PER_DAY;
 
        return { year, month, day, 0, 0, sec, digits };
    }

    [[nodiscard]] constexpr InsTime_YDoySod toYDoySod(TimeSystem timesys = UTC) const
    {
        InsTime_YMDHMS yearMonthDayHMS = toYMDHMS();
 
        auto year = yearMonthDayHMS.year;
        long double sod = static_cast<long double>(yearMonthDayHMS.hour * InsTimeUtil::SECONDS_PER_HOUR
                                                   + yearMonthDayHMS.min * InsTimeUtil::SECONDS_PER_MINUTE)
                          + yearMonthDayHMS.sec
                          + static_cast<long double>(differenceToUTC(timesys));
 
        int32_t doy = 0;
        for (int32_t i = 1; i < yearMonthDayHMS.month; i++)
        {
            doy += InsTimeUtil::daysInMonth(i, year);
        }
        doy += yearMonthDayHMS.day;
 
        return { year, doy, sod };
    }

    [[nodiscard]] constexpr InsTime toFullDay() const
    {
        return InsTime(InsTime_MJD(_mjd.mjd_day, 0.0L));
    }

    [[nodiscard]] constexpr long double toUnixTime() const
    {
        return static_cast<long double>((_mjd.mjd_day - InsTimeUtil::DIFF_TO_1_1_1970_MJD) * InsTimeUtil::SECONDS_PER_DAY)
               + _mjd.mjd_frac * InsTimeUtil::SECONDS_PER_DAY;
    }
 
    /* ----------------------------- Leap functions ----------------------------- */

    [[nodiscard]] constexpr uint16_t leapGps2UTC() const
    {
        return leapGps2UTC(_mjd);
    }

    static constexpr uint16_t leapGps2UTC(const InsTime& insTime)
    {
        return leapGps2UTC(insTime._mjd);
    }

    static constexpr uint16_t leapGps2UTC(const InsTime_GPSweekTow& gpsWeekTow)
    {
        return leapGps2UTC(InsTime(gpsWeekTow).toMJD());
    }

    static constexpr uint16_t leapGps2UTC(const InsTime_YDoySod& yearDoySod)
    {
        return leapGps2UTC(InsTime(yearDoySod).toMJD());
    }

    static constexpr uint16_t leapGps2UTC(const InsTime_YMDHMS& yearMonthDayHMS)
    {
        return leapGps2UTC(InsTime(yearMonthDayHMS).toMJD());
    }

    static constexpr uint16_t leapGps2UTC(const InsTime_JD& jd)
    {
        return leapGps2UTC(InsTime(jd).toMJD());
    }

    static constexpr uint16_t leapGps2UTC(const InsTime_MJD& mjd_in)
    {
        return static_cast<uint16_t>(std::ranges::upper_bound(InsTimeUtil::GPS_LEAP_SEC_MJD, mjd_in.mjd_day) - InsTimeUtil::GPS_LEAP_SEC_MJD.begin() - 1);
    }

    [[nodiscard]] constexpr bool isLeapYear() const
    {
        return InsTimeUtil::isLeapYear(toYMDHMS().year);
    }
 
    /* --------------------- Comparison operator overloading -------------------- */

    constexpr bool operator==(const InsTime& rhs) const { return _mjd == rhs._mjd; }
    constexpr bool operator!=(const InsTime& rhs) const { return !(*this == rhs); }
    constexpr bool operator<=(const InsTime& rhs) const { return *this < rhs || *this == rhs; }
    constexpr bool operator>=(const InsTime& rhs) const { return *this > rhs || *this == rhs; }
    constexpr bool operator<(const InsTime& rhs) const { return _mjd < rhs._mjd; }
    constexpr bool operator>(const InsTime& rhs) const { return !(*this <= rhs); }
 
    /* --------------------- Arithmetic operator overloading -------------------- */

    constexpr friend std::chrono::duration<long double> operator-(const InsTime& lhs, const InsTime& rhs)
    {
        auto diffDays = lhs._mjd.mjd_day - rhs._mjd.mjd_day;
        auto diffFrac = lhs._mjd.mjd_frac - rhs._mjd.mjd_frac;
        long double diffSec = (diffFrac + static_cast<long double>(diffDays)) * InsTimeUtil::SECONDS_PER_DAY;
        return std::chrono::duration<long double>(diffSec);
    }

    constexpr InsTime& operator+=(const std::chrono::duration<long double>& duration)
    {
        auto duration_mjd_frac = std::chrono::duration<long double, std::ratio<InsTimeUtil::SECONDS_PER_DAY>>(duration).count();
        this->_mjd = InsTime_MJD(this->_mjd.mjd_day,
                                 this->_mjd.mjd_frac + duration_mjd_frac);
        return *this;
    }

    constexpr InsTime& operator-=(const std::chrono::duration<long double>& duration)
    {
        auto duration_mjd_frac = std::chrono::duration<long double, std::ratio<InsTimeUtil::SECONDS_PER_DAY>>(duration).count();
        this->_mjd = InsTime_MJD(this->_mjd.mjd_day,
                                 this->_mjd.mjd_frac - duration_mjd_frac);
        return *this;
    }

    constexpr friend InsTime operator+(const InsTime& time, const std::chrono::duration<long double>& duration)
    {
        return InsTime(time) += duration;
    }

    constexpr friend InsTime operator-(const InsTime& time, const std::chrono::duration<long double>& duration)
    {
        return InsTime(time) -= duration;
    }
 
    /* ---------------------------- Utility Functions --------------------------- */

    explicit operator std::string() const;

    [[nodiscard]] constexpr bool empty() const
    {
        return _mjd.mjd_day == 0 && _mjd.mjd_frac == 0.0L;
    }

    void reset()
    {
        _mjd.mjd_day = 0;
        _mjd.mjd_frac = 0.0L;
    }

    void MakeTimeFromGloOrbit(double UTC_sec)
    {
        auto ymdhms = toYMDHMS();
        // difference between toe (OBRIT-0 last element) and toc (ORBIT-0 first element) in seconds
        long double diff = gcem::fmod(static_cast<long double>(UTC_sec), InsTimeUtil::SECONDS_PER_DAY)
                           - (ymdhms.hour * InsTimeUtil::SECONDS_PER_HOUR
                              + ymdhms.min * InsTimeUtil::SECONDS_PER_MINUTE
                              + ymdhms.sec);
        // std::cout << "orbit diff " << diff << "\n";
        *this += std::chrono::duration<long double>(diff);
    }

    [[nodiscard]] constexpr int differenceToUTC(TimeSystem timesys) const
    {
        switch (TimeSystem_(timesys))
        {
        case GPST: // = GPS Time (UTC + leap_seconds)
            return this->leapGps2UTC();
        case GLNT: // = GLONASS Time (UTC+ 3h)
            return 3 * InsTimeUtil::SECONDS_PER_HOUR;
        case GST: // = GALILEO Time (~ GPS) (UTC = GST - 18) is synchronized with TAI with a nominal offset below 50 ns
            return this->leapGps2UTC();
        case QZSST:
            return this->leapGps2UTC(); // TODO citation for synchronization accuracy
        case IRNSST:
            return this->leapGps2UTC(); // TODO citation for synchronization accuracy
        case BDT:                       // = BeiDou Time (UTC) is synchronized with UTC within 100 ns<
            return this->leapGps2UTC() - 14;
        case UTC:
        case TimeSys_None:
            return 0;
        }
        return 0;
    }
 
  private:
    InsTime_MJD _mjd;
};

std::ostream& operator<<(std::ostream& os, const InsTime_MJD& mjd);
std::ostream& operator<<(std::ostream& os, const InsTime_JD& jd);
std::ostream& operator<<(std::ostream& os, const InsTime_GPSweekTow& gpsWeekTow);
std::ostream& operator<<(std::ostream& os, const InsTime_YMDHMS& ymdhms);
std::ostream& operator<<(std::ostream& os, const InsTime_YDoySod& yDoySod);
std::ostream& operator<<(std::ostream& os, const InsTime& insTime);

void to_json(json& j, const InsTime& insTime);
void from_json(const json& j, InsTime& insTime);
 
} // namespace NAV
 
namespace std
{
template<>
struct hash<NAV::InsTime>
{
    std::size_t operator()(const NAV::InsTime& t) const
    {
        auto hash1 = std::hash<int32_t>{}(t.toMJD().mjd_day);
        auto hash2 = std::hash<long double>{}(t.toMJD().mjd_frac);
 
        return hash1 | (hash2 << 32);
    }
};
} // namespace std
 
#ifndef DOXYGEN_IGNORE

template<>
struct fmt::formatter<NAV::InsTime_MJD> : fmt::formatter<std::string>
{
    template<typename FormatContext>
    auto format(const NAV::InsTime_MJD& mjd, FormatContext& ctx) const
    {
        return fmt::format_to(ctx.out(), "day={}, frac={}", mjd.mjd_day, static_cast<double>(mjd.mjd_frac));
    }
};
template<>
struct fmt::formatter<NAV::InsTime_JD> : fmt::formatter<std::string>
{
    template<typename FormatContext>
    auto format(const NAV::InsTime_JD& jd, FormatContext& ctx) const
    {
        return fmt::format_to(ctx.out(), "day={}, frac={}", jd.jd_day, static_cast<double>(jd.jd_frac));
    }
};
template<>
struct fmt::formatter<NAV::InsTime_GPSweekTow> : fmt::formatter<std::string>
{
    template<typename FormatContext>
    auto format(const NAV::InsTime_GPSweekTow& gpsWeekTow, FormatContext& ctx) const
    {
        return fmt::format_to(ctx.out(), "cycle={}, week={}, tow={}",
                              gpsWeekTow.gpsCycle, gpsWeekTow.gpsWeek, static_cast<double>(gpsWeekTow.tow));
    }
};
template<>
struct fmt::formatter<NAV::InsTime_YMDHMS> : fmt::formatter<std::string>
{
    template<typename FormatContext>
    auto format(const NAV::InsTime_YMDHMS& ymdhms, FormatContext& ctx) const
    {
        return fmt::format_to(ctx.out(), "{}-{}-{} {}:{}:{:.6g}",
                              ymdhms.year, ymdhms.month, ymdhms.day,
                              ymdhms.hour, ymdhms.min, static_cast<double>(ymdhms.sec));
    }
};
template<>
struct fmt::formatter<NAV::InsTime_YDoySod> : fmt::formatter<std::string>
{
    template<typename FormatContext>
    auto format(const NAV::InsTime_YDoySod& yDoySod, FormatContext& ctx) const
    {
        return fmt::format_to(ctx.out(), "year={}, doy={}, sod={:.6g}",
                              yDoySod.year, yDoySod.doy, static_cast<double>(yDoySod.sod));
    }
};
template<>
struct fmt::formatter<NAV::InsTime> : fmt::formatter<std::string>
{
    template<typename FormatContext>
    auto format(const NAV::InsTime& insTime, FormatContext& ctx) const
    {
        return fmt::format_to(ctx.out(), "{} ({})", insTime.toYMDHMS(), insTime.toGPSweekTow());
    }
};
 
#endif