From 95fcca2deb659adfc661f5d41766ff6e28a349f8 Mon Sep 17 00:00:00 2001 From: "Timur A. Fatkhullin" Date: Fri, 7 Mar 2025 00:28:17 +0300 Subject: [PATCH] ... --- cxx/mount_astrom.h | 213 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 210 insertions(+), 3 deletions(-) diff --git a/cxx/mount_astrom.h b/cxx/mount_astrom.h index 7a93ac3..497e130 100644 --- a/cxx/mount_astrom.h +++ b/cxx/mount_astrom.h @@ -7,7 +7,9 @@ *********************************/ #include +#include #include +#include "utils.h" namespace mcc::traits @@ -16,6 +18,12 @@ namespace mcc::traits template concept mcc_scalar_or_simd_c = xsimd::is_batch::value || std::is_arithmetic_v; + +template +concept mcc_real_or_char_range = + std::floating_point || + (std::ranges::contiguous_range && std::same_as>>); + } // namespace mcc::traits @@ -93,19 +101,218 @@ static int mcc_julday(const std::chrono::system_clock::time_point& start_time, return 0; } -template -T mcc_time_to_alt_limit(const T& alt_limit, const T& RA, const T& DEC, const std::chrono::time_point& now) +/* + * angles are in degrees or sexagimal string form + */ +template +double mcc_time_to_alt_limit(const AT& alt_limit, + const RT& RA, + const DT& DEC, + const LT& LAT, + const std::chrono::system_clock::time_point& now) { // sin(alt) = sin(DEC)*sin(phi) + cos(DEC)*cos(phi)*cos(HA) // HA = LST - RA + + double ra, dec, lat, alt; + + if constexpr (std::floating_point) { + alt = alt_limit * utils::deg2radCoeff; + } else { + alt = utils::parsAngleString(alt_limit); + alt *= utils::deg2radCoeff; + } + + if constexpr (std::floating_point) { + ra = RA * utils::deg2radCoeff; + } else { + ra = utils::parsAngleString(RA, true); + ra *= utils::deg2radCoeff; + } + + if constexpr (std::floating_point
) { + dec = DEC * utils::deg2radCoeff; + } else { + dec = utils::parsAngleString(DEC); + dec *= utils::deg2radCoeff; + } + + if constexpr (std::floating_point) { + lat = LAT * utils::deg2radCoeff; + } else { + lat = utils::parsAngleString(LAT); + lat *= utils::deg2radCoeff; + } + + double cos_ha = (std::sin(alt) - std::sin(dec) * std::sin(lat)) / std::cos(dec) / std::cos(lat); + if (std::abs(cos_ha) > 1.0) { // circumpolar (it never sets below horizon) + return std::numeric_limits::infinity(); + } + + double lst = std::acos(cos_ha) + ra; + + + return 0.0; } -/* IERS bulletines +/* IERS bulletins * * BULLETIN A: https://datacenter.iers.org/data/latestVersion/bulletinA.txt * leapseconds: https://hpiers.obspm.fr/iers/bul/bulc/Leap_Second.dat * */ +struct leapsecond_db_elem_t { + double mjd; + unsigned day, month; + int year; + double tai_utc; // TAI-UTC in seconds +}; + +typedef std::vector leapsecond_db_t; + + +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 +}; + + +typedef std::vector earth_orient_db_t; + +// init to some known state +static leapsecond_db_t CURRENT_LEAPSECONDS_DB = { + {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} + +}; + +bool mcc_parse_bulletinA(std::derived_from> auto& stream, char comment_sym = '*') +{ + for (std::string line; std::getline(stream, line);) { + if (line.empty()) { + continue; + } + + auto sv = utils::trimSpaces(line, utils::TrimType::TRIM_BOTH); + + if (sv.size()) { + if (sv[0] == comment_sym) { // comment string + continue; + } + } else { + continue; + } + } + + return true; +} + + +bool mcc_parse_leapsecs(std::derived_from> auto& stream, char comment_sym = '#') +{ + size_t n; + std::optional vd; + std::optional vui; + std::optional vi; + + leapsecond_db_t db; + leapsecond_db_elem_t db_elem; + + for (std::string line; std::getline(stream, line);) { + if (line.empty()) { + continue; + } + + auto sv = utils::trimSpaces(line, utils::TrimType::TRIM_BOTH); + + if (sv.size()) { + if (sv[0] == comment_sym) { // comment string + continue; + } + } else { + continue; + } + + n = 0; + for (auto const& el : std::views::split(sv, std::string_view(" "))) { + if (std::ranges::size(el) == 0) { + continue; + } + + switch (n) { + case 0: + vd = utils::numFromStr(el); + if (!vd) { + return false; + } + db_elem.mjd = vd.value(); + break; + case 1: + vui = utils::numFromStr(el); + if (!vui) { + return false; + } + db_elem.day = vui.value(); + break; + case 2: + vui = utils::numFromStr(el); + if (!vui) { + return false; + } + db_elem.month = vui.value(); + break; + case 3: + vi = utils::numFromStr(el); + if (!vi) { + return false; + } + db_elem.year = vi.value(); + break; + case 4: + vd = utils::numFromStr(el); + if (!vd) { + return false; + } + db_elem.tai_utc = vd.value(); + break; + default: + // just ignore??!!! + break; + } + + ++n; + + if (n == 5) { + break; + } + } + + if (n < 5) { // not enough elements + return false; + } + + db.emplace_back(std::move(db_elem)); + } + + CURRENT_LEAPSECONDS_DB = std::move(db); + + return true; +} + } // namespace mcc::astro