diff --git a/cxx/mcc_spdlog.h b/cxx/mcc_spdlog.h index 34da56e..19f7325 100644 --- a/cxx/mcc_spdlog.h +++ b/cxx/mcc_spdlog.h @@ -141,7 +141,14 @@ protected: std::shared_ptr _loggerSPtr; - // helper method + // helper methods + + auto getThreadId() const + { + std::ostringstream st; + st << std::this_thread::get_id(); + return st.str(); + } // 'after_idx' is 0-based index! void addMarkToPatternIdx(const traits::mcc_input_char_range auto& mark, size_t after_idx = 1) diff --git a/cxx/mount.h b/cxx/mount.h index 62dfa4b..5c6bcce 100644 --- a/cxx/mount.h +++ b/cxx/mount.h @@ -151,78 +151,33 @@ public: addMarkToPatternIdx(logger_mark); + logDebug("Create MccMount class instance: thread = {}", getThreadId()); - logDebug("Create MccMount class instance: thread = {}", std::this_thread::get_id()); - - - // init time scales related databases to default (pre-defined) state - - logInfo("initializing leap seconds database to default state ..."); - strst.str(defaults::MCC_DEFAULT_LEAP_SECONDS_FILE); - _leapSecondsDB = astro::mcc_parse_leapsecs(strst); - logInfo("leap seconds default database expired date: {}", _leapSecondsDB.expireDate); - - logInfo("initializing Earth orientation (pole coordinates, UT1-UTC) database to default state ..."); - strst.clear(); - strst.str(defaults::MCC_DEFAULT_IERS_BULLETIN_A_FILE); - _earthOrientDB = astro::mcc_parse_bulletinA(strst); - logInfo("Earth orientation default database (Bulletin A) date: {}", _earthOrientDB.bulletinDate); - - // load time scales relates databases from files - std::ifstream fst; logInfo("Load leap seconds and Earth orientation databases ..."); - auto time_db_loader = [&fst, this](const std::string& filename, std::string_view type, auto& db) { + auto time_db_loader = [this](const std::string& filename, std::string_view type, auto& db) { if (filename.empty()) { logWarn("An empty {} filename! Skip and keep default values!", type); return; } - fst.open(filename); - if (!fst.is_open()) { - logError("CANNOT open {} file '{}'!", type, filename); - logWarn("Keep {} database in default state!", type); - return; - } - - if constexpr (std::same_as>) { - db = astro::mcc_parse_leapsecs(fst); - } else if constexpr (std::same_as>) { - db = astro::mcc_parse_bulletinA(fst); - } else { - static_assert(false, "INVALID DATABASE TYPE!!!"); - } - - if (db.state != astro::IERS_DB_STATE_OK) { - logError("CANNOT parse {} file '{}'!", type, filename); + bool ok = db.load(filename); + if (!ok) { + logError("CANNOT parse {} file '{}' or it is not accessable!", type, filename); logWarn("Keep {} database in default state!", type); } else { logInfo("{} database was successfully loaded from '{}' file", type, filename); } - - fst.close(); }; - astro::leapsecond_db_t ldb; - astro::earth_orient_db_t edb; - - time_db_loader(_mountCurrentConfig.leap_seconds_filename, "leap seconds", ldb); - if (ldb.state == astro::IERS_DB_STATE_OK) { - _leapSecondsDB = std::move(ldb); - logInfo("leap seconds default database expired date: {}", _leapSecondsDB.expireDate); - } - - time_db_loader(_mountCurrentConfig.earth_orient_filename, "Earth orientation", edb); - if (edb.state == astro::IERS_DB_STATE_OK) { - _earthOrientDB = std::move(edb); - logInfo("Earth orientation default database (Bulletin A) date: {}", _earthOrientDB.bulletinDate); - } + time_db_loader(_mountCurrentConfig.leap_seconds_filename, "leap seconds", _leapSecondsDB); + time_db_loader(_mountCurrentConfig.earth_orient_filename, "Earth orientation", _earthOrientDB); } virtual ~MccMount() { - logDebug("Delete MccMount class instance: thread = {}", std::this_thread::get_id()); + logDebug("Delete MccMount class instance: thread = {}", getThreadId()); } @@ -257,8 +212,8 @@ protected: std::function _exitCurrentState; // time scales related databases - astro::leapsecond_db_t _leapSecondsDB; - astro::earth_orient_db_t _earthOrientDB; + astrom::MccLeapSeconds _leapSecondsDB; + astrom::MccIersBulletinA _earthOrientDB; std::atomic _currentMountOrient; diff --git a/cxx/mount_astrom.h b/cxx/mount_astrom.h index f71949e..35e6195 100644 --- a/cxx/mount_astrom.h +++ b/cxx/mount_astrom.h @@ -33,9 +33,10 @@ concept mcc_real_or_char_range = } // namespace mcc::traits -namespace mcc::astro +namespace mcc::astrom { + // modified Julian date (based on ERFA eraCal2jd) template static int mcc_julday(const std::chrono::system_clock::time_point& start_time, @@ -64,11 +65,15 @@ static int mcc_julday(const std::chrono::system_clock::time_point& start_time, return -2; } + int64_t im = (unsigned)ymd.month(); + int64_t id = (unsigned)ymd.day(); + int64_t iy = (int)ymd.year(); + // my = (im - 14) / 12; // iypmy = (long) (iy + my); - int64_t my = -(14 - (unsigned)ymd.month()) / 12; - int64_t iypmy = (int)ymd.year() + my; + int64_t my = (im - 14LL) / 12LL; + int64_t iypmy = iy + my; // (1461L * (iypmy + 4800L)) / 4L // + (367L * (long) (im - 2 - 12 * my)) / 12L @@ -76,17 +81,14 @@ static int mcc_julday(const std::chrono::system_clock::time_point& start_time, // + (long) id - 2432076L // integer part of result MJD - int64_t mjd_int = 1461LL * (iypmy + 480LL) / 4LL + - (367LL * ((int64_t)(unsigned)ymd.month() - 2LL - 12LL * my)) / 12LL - - (3LL * (iypmy + 4900LL) / 100LL) / 4LL + (int64_t)(unsigned)ymd.day() - 2432076LL; + int64_t mjd_int = (1461LL * (iypmy + 4800LL)) / 4LL + (367LL * (im - 2LL - 12LL * my)) / 12LL - + (3LL * ((iypmy + 4900LL) / 100LL)) / 4LL + id - 2432076LL; - constexpr double nanosec = 1.0 / 24.0 / 3600.0 / 1.0E-9; // 1 nanosecond in days - - double mjd_float = static_cast(mjd_int) + - std::chrono::duration_cast(start_time - dd).count() * nanosec; + using fd_t = std::chrono::duration>; // fraction of day + double mjd_float = static_cast(mjd_int) + std::chrono::duration_cast(start_time - dd).count(); if constexpr (std::ranges::range) { - double d_step = std::chrono::duration_cast(step).count() * nanosec; + double d_step = std::chrono::duration_cast(step).count(); size_t i = 0; @@ -211,200 +213,20 @@ double mcc_time_to_alt_limit(const AT& alt_limit, } -/* IERS bulletins + +/* Classes to represent IERS bulletins * * BULLETIN A: https://datacenter.iers.org/data/latestVersion/bulletinA.txt * leapseconds: https://hpiers.obspm.fr/iers/bul/bulc/Leap_Second.dat * */ -enum iers_db_state_t { IERS_DB_STATE_UNINITIALIZED, IERS_DB_STATE_OK, IERS_DB_STATE_PARSE_ERROR }; -struct leapsecond_db_elem_t { - double mjd; - unsigned day, month; - int year; - // std::chrono::year_month_day ymd; - double tai_utc; // TAI-UTC in seconds -}; - -struct leapsecond_db_t { - iers_db_state_t state = IERS_DB_STATE_UNINITIALIZED; - std::chrono::system_clock::time_point expireDate{}; - std::vector db{}; -}; - - - -struct earth_orient_db_elem_t { - int year; - unsigned month, day; - double mjd; - double x, y; // Polar coordinates in arcsecs - double dut1; // UT1-UTC in seconds -}; - - -struct earth_orient_db_t { - iers_db_state_t state = IERS_DB_STATE_UNINITIALIZED; - std::chrono::system_clock::time_point bulletinDate{}; - double tt_tai = 0.0; // TT-TAI - std::vector db{}; -}; - - - -static earth_orient_db_t mcc_parse_bulletinA(std::derived_from> auto& stream, - char comment_sym = '*') -{ - const std::regex bull_date_rx{ - "^ *[0-9]{1,2} +(January|February|March|April|May|June|July|August|September|October|November|December) " - "+[0-9]{4,} +Vol\\. +[XMLCDVI]+ +No\\. +[0-9]+ *$"}; - - const std::regex bull_tt_tai_rx{"^ *TT += +TAI +\\+ +[0-9]+\\.[0-9]+ +seconds *$"}; - - const std::regex bull_tab_title_rx{"^ *MJD +x\\(arcsec\\) +y\\(arcsec\\) +UT1-UTC\\(sec\\) *$"}; - - // 2025 3 7 60741 0.0663 0.3341 0.04348 - const std::regex bull_tab_vals_rx{ - "^ *[0-9]{4,} +[0-9]{1,2} +[0-9]{1,2} +[0-9]{5,} +[0-9]+\\.[0-9]+ +[0-9]+\\.[0-9]+ +[0-9]+\\.[0-9]+ *$"}; - - earth_orient_db_t db; - enum { TAB_STATE_SEEK, TAB_STATE_START }; - int tab_state = TAB_STATE_SEEK; - - int year; - unsigned month, day; - double mjd, x, y, dut1; - std::istringstream is; - - for (std::string line; std::getline(stream, line);) { - if (line.empty()) { - continue; - } - - auto sv = utils::trimSpaces(line, utils::TrimType::TRIM_LEFT); - - if (sv.size()) { - if (sv[0] == comment_sym) { // comment string - continue; - } - - if (tab_state == TAB_STATE_START) { - if (std::regex_match(sv.begin(), sv.end(), bull_tab_vals_rx)) { - // is.str({sv.begin(), sv.end()}); - is.str(line); - is >> year >> month >> day >> mjd >> x >> y >> dut1; - db.db.emplace_back(year, month, day, mjd, x, y, dut1); - is.clear(); - } else { // end of the table - just stop parsing - break; - } - - continue; - } - - if (std::regex_match(sv.begin(), sv.end(), bull_date_rx)) { - is.str({sv.begin(), sv.end()}); - is >> std::chrono::parse("%d %B %Y", db.bulletinDate); - continue; - } - - if (std::regex_match(sv.begin(), sv.end(), bull_tt_tai_rx)) { - is.str({sv.begin(), sv.end()}); - std::string dummy; - is >> dummy >> dummy >> dummy >> dummy >> db.tt_tai; - continue; - } - - if (std::regex_match(sv.begin(), sv.end(), bull_tab_title_rx)) { - tab_state = TAB_STATE_START; - continue; - } - - } else { // empty string (only spaces) - continue; - } - } - - if (db.db.empty()) { - db.state = IERS_DB_STATE_PARSE_ERROR; - } else { - db.state = IERS_DB_STATE_OK; - } - - return db; -} - - -static leapsecond_db_t mcc_parse_leapsecs(std::derived_from> auto& stream, - char comment_sym = '#') -{ - // # File expires on 28 December 2025 - const std::regex expr_date_rx{ - "^ *# *File +expires +on +[0-8]{1,2} " - "+(January|February|March|April|May|June|July|August|September|October|November|December) +[0-9]{4} *$"}; - - const std::regex data_rx{"^ *[0-9]{5,}(\\.?[0-9]+) +[0-9]{1,2} +[0-9]{1,2} +[0-9]{4} +[0-9]{1,} *$"}; - - std::istringstream is; - double mjd; - unsigned day, month; - int year; - double tai_utc; - - - leapsecond_db_t db; - - for (std::string line; std::getline(stream, line);) { - auto sv = utils::trimSpaces(line, utils::TrimType::TRIM_LEFT); - - if (sv.size()) { - if (sv[0] == comment_sym) { // comment string - if (std::regex_match(line, expr_date_rx)) { - auto pos = line.find("on"); - sv = utils::trimSpaces(std::string_view{line.begin() + pos + 2, line.end()}, - utils::TrimType::TRIM_LEFT); - is.str({sv.begin(), sv.end()}); - is >> std::chrono::parse("%d %B %Y", db.expireDate); - is.clear(); - } - continue; - } - } else { - continue; - } - - if (std::regex_match(line, data_rx)) { - is.str(line); - is >> mjd >> day >> month >> year >> tai_utc; - db.db.emplace_back(mjd, day, month, year, tai_utc); - is.clear(); - continue; - } - } - - if (db.db.empty()) { - db.state = IERS_DB_STATE_PARSE_ERROR; - } else { - db.state = IERS_DB_STATE_OK; - } - - return db; -} - - - -} // namespace mcc::astro - - -namespace mcc::astrom -{ - - -class MccLeapSeconds +class MccLeapSeconds final { public: + typedef std::chrono::system_clock::time_point time_point_t; + MccLeapSeconds() { // create default values @@ -416,7 +238,7 @@ public: ~MccLeapSeconds() = default; - std::chrono::system_clock::time_point expireDate() const + time_point_t expireDate() const { return _expireDate; } @@ -494,7 +316,7 @@ public: } - std::optional operator[](const std::chrono::system_clock::time_point& tp) const + std::optional operator[](const time_point_t& tp) const { if (tp > _expireDate) { // ???????!!!!!!!!!!! return std::nullopt; @@ -516,7 +338,7 @@ public: { double e_mjd; - astro::mcc_julday(_expireDate, e_mjd); + astrom::mcc_julday(_expireDate, e_mjd); if (mjd > e_mjd) { // ???????!!!!!!!!!!! return std::nullopt; // return _db.back().tai_utc; @@ -531,6 +353,16 @@ public: return std::nullopt; } + + void dump(std::derived_from> auto& stream) const + { + stream << std::format("Leap seconds database expire date: {}", _expireDate) << '\n'; + for (auto const& el : _db) { + stream << std::format("{} {} {}", el.mjd, el.ymd, el.tai_utc) << '\n'; + } + } + + private: inline static const std::regex expr_date_rx{ "^ *# *File +expires +on +[0-8]{1,2} " @@ -538,7 +370,7 @@ private: inline static const std::regex data_rx{"^ *[0-9]{5,}(\\.?[0-9]+) +[0-9]{1,2} +[0-9]{1,2} +[0-9]{4} +[0-9]{1,} *$"}; - std::chrono::system_clock::time_point _expireDate{}; + time_point_t _expireDate{}; struct leapsecond_db_elem_t { double mjd; @@ -551,37 +383,74 @@ private: -class MccIersBulletinA +class MccIersBulletinA final { public: + typedef std::chrono::system_clock::time_point time_point_t; + struct pole_pos_t { double x, y; }; + struct date_range_t { + std::chrono::year_month_day begin; + std::chrono::year_month_day end; + }; + + struct date_range_mjd_t { + double begin; + double end; + }; + MccIersBulletinA() { + // create pre-defined (default-state) database std::istringstream ist(defaults::MCC_DEFAULT_IERS_BULLETIN_A_FILE); load(ist); } + ~MccIersBulletinA() = default; + std::chrono::system_clock::time_point bulletinDate() const { return _date; } + + date_range_t dateRange() const + { + return {_db.front().ymd, _db.back().ymd}; + } + + date_range_mjd_t dateRangeMJD() const + { + return {_db.front().mjd, _db.back().mjd}; + } + + double TT_TAI() const { return _tt_tai; } // DUT1 = UT1 - UTC - std::optional DUT1(const std::chrono::system_clock::time_point& tp) const + std::optional DUT1(const time_point_t& tp) const { - std::chrono::year_month_day ymd{std::chrono::floor(tp)}; + // use of the closest date + std::chrono::year_month_day ymd{std::chrono::round(tp)}; - for (auto const& el : _db | std::views::reverse) { - if (ymd >= el.ymd) { + if (ymd < _db.front().ymd) { + return std::nullopt; + } + + if (ymd > _db.back().ymd) { + return std::nullopt; + } + + auto el_prev = _db.front(); + for (auto const& el : _db) { + if (ymd <= el.ymd) { return el.dut1; } } @@ -589,10 +458,20 @@ public: return std::nullopt; } - std::optional DUT1(const double& mjd) const + std::optional DUT1(double mjd) const { - for (auto const& el : _db | std::views::reverse) { - if (mjd >= el.mjd) { + mjd = std::round(mjd); // round to closest integer MJD + + if (mjd < _db.front().mjd) { + return std::nullopt; + } + + if (mjd > _db.back().mjd) { + return std::nullopt; + } + + for (auto const& el : _db) { + if (mjd <= el.mjd) { return el.dut1; } } @@ -600,12 +479,20 @@ public: return std::nullopt; } - std::optional polePos(const std::chrono::system_clock::time_point& tp) const + std::optional polarCoords(const time_point_t& tp) const { - std::chrono::year_month_day ymd{std::chrono::floor(tp)}; + std::chrono::year_month_day ymd{std::chrono::round(tp)}; - for (auto const& el : _db | std::views::reverse) { - if (ymd >= el.ymd) { + if (ymd < _db.front().ymd) { + return std::nullopt; + } + + if (ymd > _db.back().ymd) { + return std::nullopt; + } + + for (auto const& el : _db) { + if (ymd <= el.ymd) { return pole_pos_t{el.x, el.y}; } } @@ -613,10 +500,20 @@ public: return std::nullopt; } - std::optional polePos(const double& mjd) const + std::optional polarCoords(double mjd) const { - for (auto const& el : _db | std::views::reverse) { - if (mjd >= el.mjd) { + mjd = std::round(mjd); // round to closest integer MJD + + if (mjd < _db.front().mjd) { + return std::nullopt; + } + + if (mjd > _db.back().mjd) { + return std::nullopt; + } + + for (auto const& el : _db) { + if (mjd <= el.mjd) { return pole_pos_t{el.x, el.y}; } } @@ -714,6 +611,16 @@ public: return ok; } + + void dump(std::derived_from> auto& stream) const + { + stream << std::format("Bulletin A issue date: {}", _date) << '\n'; + stream << std::format("TT-TAI: {}", _tt_tai) << '\n'; + for (auto const& el : _db) { + stream << std::format("{} {} {:6.4f} {:6.4f} {:7.5f}", el.mjd, el.ymd, el.x, el.y, el.dut1) << '\n'; + } + } + private: inline static const std::regex bull_date_rx{ "^ *[0-9]{1,2} +(January|February|March|April|May|June|July|August|September|October|November|December) " @@ -727,7 +634,7 @@ private: "^ *[0-9]{4,} +[0-9]{1,2} +[0-9]{1,2} +[0-9]{5,} +[0-9]+\\.[0-9]+ +[0-9]+\\.[0-9]+ +[0-9]+\\.[0-9]+ *$"}; - std::chrono::system_clock::time_point _date; + time_point_t _date; double _tt_tai; struct earth_orient_db_elem_t { diff --git a/cxx/tests/astrom_test.cpp b/cxx/tests/astrom_test.cpp index 4309860..4ebae73 100644 --- a/cxx/tests/astrom_test.cpp +++ b/cxx/tests/astrom_test.cpp @@ -4,51 +4,6 @@ #include "../mount_astrom.h" -static std::string leap_secs_file = R"--( -# Value of TAI-UTC in second valid beetween the initial value until -# the epoch given on the next line. The last line reads that NO -# leap second was introduced since the corresponding date -# Updated through IERS Bulletin 69 issued in January 2025 -# -# -# File expires on 28 December 2025 -# -# -# MJD Date TAI-UTC (s) -# day month year -# --- -------------- ------ -# - 41317.0 1 1 1972 10 - 41499.0 1 7 1972 11 - 41683.0 1 1 1973 12 - 42048.0 1 1 1974 13 - 42413.0 1 1 1975 14 - 42778.0 1 1 1976 15 - 43144.0 1 1 1977 16 - 43509.0 1 1 1978 17 - 43874.0 1 1 1979 18 - 44239.0 1 1 1980 19 - 44786.0 1 7 1981 20 - 45151.0 1 7 1982 21 - 45516.0 1 7 1983 22 - 46247.0 1 7 1985 23 - 47161.0 1 1 1988 24 - 47892.0 1 1 1990 25 - 48257.0 1 1 1991 26 - 48804.0 1 7 1992 27 - 49169.0 1 7 1993 28 - 49534.0 1 7 1994 29 - 50083.0 1 1 1996 30 - 50630.0 1 7 1997 31 - 51179.0 1 1 1999 32 - 53736.0 1 1 2006 33 - 54832.0 1 1 2009 34 - 56109.0 1 7 2012 35 - 57204.0 1 7 2015 36 - 57754.0 1 1 2017 37 -)--"; - - int main(int argc, char* argv[]) { int ecode = 0; @@ -64,37 +19,39 @@ int main(int argc, char* argv[]) exit(1); } - auto db_a = mcc::astro::mcc_parse_bulletinA(ist); - if (db_a.state != mcc::astro::IERS_DB_STATE_OK) { - std::cout << "Cannot parse input IERS Bulletin A file!\n"; - ecode = 1; - } else { - std::cout << "IERS Bulletin A data:\n"; - std::cout << "Date: " << db_a.bulletinDate << "\n"; - std::cout << "TT-TAI: " << db_a.tt_tai << "\n"; - for (auto& el : db_a.db) { - std::cout << "MJD: " << el.mjd << "\tDUT1 = " << el.dut1 << "\n"; - } - } + mcc::astrom::MccIersBulletinA bullA; + bullA.dump(std::cout); + + std::cout << "\n\n"; + + mcc::astrom::MccLeapSeconds lps; + lps.dump(std::cout); + + std::cout << "\n\n"; + + bool ok = bullA.load(ist); + bullA.dump(std::cerr); ist.close(); - std::cout << "\n\n\n"; + std::cout << "\n\n\n------------------\n"; - std::istringstream isst(leap_secs_file); - auto db_ls = mcc::astro::mcc_parse_leapsecs(isst); + auto now = std::chrono::system_clock::now(); + std::cout << "LEAP SECS for now: " << lps[now].value_or(std::numeric_limits::quiet_NaN()) << "\n"; + + std::cout << "DUT1 for now: " << bullA.DUT1(now).value_or(std::numeric_limits::quiet_NaN()) << " (" << now + << ")\n"; + + double mjd = 61077.4; + + std::cout << "DUT1 for MJD = " << mjd << ": " << bullA.DUT1(mjd).value_or(std::numeric_limits::quiet_NaN()) + << "\n"; + + + auto st = mcc::astrom::mcc_julday(now, mjd); + std::cout << "MJD for now: " << std::setprecision(19) << mjd << " (" << now << ")\n"; - if (db_ls.state != mcc::astro::IERS_DB_STATE_OK) { - std::cout << "Cannot parse input IERS leap seconds file!\n"; - ecode = 1; - } else { - std::cout << "IERS leap seconds data:\n"; - std::cout << "Expire date: " << db_ls.expireDate << "\n"; - for (auto& el : db_ls.db) { - std::cout << "MJD: " << el.mjd << "\tTAI-UTC = " << el.tai_utc << "\n"; - } - } return ecode; }