...

cxx/CMakeLists.txt

@@ -94,9 +94,10 @@ ExternalProject_Add(erfa_lib
     BUILD_COMMAND ninja -C <BINARY_DIR>
     INSTALL_COMMAND meson install -C <BINARY_DIR>
 )
-add_library(erfa STATIC IMPORTED)
+add_library(ERFA_LIB STATIC IMPORTED)
-set_target_properties(erfa PROPERTIES IMPORTED_LOCATION ${CMAKE_BINARY_DIR}/erfa_lib/liberfa.a)
+set_target_properties(ERFA_LIB PROPERTIES IMPORTED_LOCATION ${CMAKE_BINARY_DIR}/erfa_lib/liberfa.a)
-
+set(ERFA_INCLUDE_DIR ${CMAKE_BINARY_DIR}/erfa_lib)
+include_directories(${ERFA_INCLUDE_DIR})
 
 option(WITH_TESTS "Build tests" ON)
 
@@ -118,7 +119,8 @@ set(MOUNT_SERVER_APP_SRC mount.h mount_state.h mount_server.cpp comm_server.h co
     mount_astrom_default.h)
 set(MOUNT_SERVER_APP mount_server)
 add_executable(${MOUNT_SERVER_APP} ${MOUNT_SERVER_APP_SRC})
-target_link_libraries(${MOUNT_SERVER_APP} ${CNTR_PROTO_LIB} spdlog::spdlog_header_only)
+# target_include_directories(${MOUNT_SERVER_APP} PUBLIC ${ERFA_INCLUDE_DIR})
+target_link_libraries(${MOUNT_SERVER_APP} ${CNTR_PROTO_LIB} spdlog::spdlog_header_only ${ERFA_LIB})
 
 if (WITH_TESTS)
     set(CNTR_PROTO_TEST_APP cntr_proto_test)

cxx/comm_server.h

@@ -33,6 +33,7 @@
 
 #include "comm_server_endpoint.h"
 #include "control_proto.h"
+#include "mcc_traits.h"
 #include "mount.h"
 
 namespace mcc
@@ -75,14 +76,6 @@ static constexpr bool is_local_seqpack_proto = std::derived_from<T, asio::local:
                                                std::derived_from<T, asio::local::seq_packet_protocol::socket>;
 
 
-template <typename T>
-concept mcc_time_duration_c = requires {
-    []<class RT, class PT>(std::type_identity<std::chrono::duration<RT, PT>>) {
-
-    }(std::type_identity<std::remove_cvref_t<T>>());
-};
-
-
 }  // namespace traits
 
 

cxx/mcc_traits.h

@@ -1,5 +1,6 @@
 #pragma once
 
+#include <chrono>
 #include <format>
 #include <ranges>
 
@@ -38,6 +39,15 @@ concept mcc_formattable =
     requires(T v, std::format_context ctx) { std::formatter<std::remove_cvref_t<T>>().format(v, ctx); };
 
 
+template <typename T>
+// from https://stackoverflow.com/questions/74383254/concept-that-models-only-the-stdchrono-duration-types
+concept mcc_time_duration_c = requires {
+    []<class Rep, class Period>(std::type_identity<std::chrono::duration<Rep, Period>>) {
+
+    }(std::type_identity<std::remove_cvref_t<T>>());
+};
+
+
 /*  a callable concept and its signature traits   */
 
 template <typename T>

cxx/mount_astrom.h

@@ -14,6 +14,12 @@
 #include "mount_astrom_default.h"
 #include "utils.h"
 
+namespace erfa
+{
+#include <erfa.h>
+#include <erfam.h>
+}  // namespace erfa
+
 namespace mcc::traits
 {
 
@@ -36,6 +42,14 @@ concept mcc_real_or_char_range =
 namespace mcc::astrom
 {
 
+// a time duration represented in radians (the precision is near 1 nanosecond)
+// typedef std::chrono::duration<double, std::ratio<134670467, 10000000000000000>> mcc_radsec_duration_t;
+typedef std::chrono::duration<double, std::ratio<134670467, 10000000000000000>> mcc_chrono_radians;
+
+
+// https://gssc.esa.int/navipedia/index.php?title=CEP_to_ITRF
+static constexpr double mcc_UT1_to_sideral_ratio = 1.002737909350795;  // UT1/sideral
+
 
 // modified Julian date (based on ERFA eraCal2jd)
 template <traits::mcc_real_scalar_or_real_range_c ResT>
@@ -132,24 +146,34 @@ static int mcc_julday(const std::chrono::system_clock::time_point& start_time,
 
 /*
  *  angles are in degrees or sexagimal string form
+ *  RA and DEC are apparent!
  *
  *  returns
  *      NaN if object is non-rising or "alt_limit" < 0, Inf is circumpolar
  */
-template <traits::mcc_real_or_char_range AT,
+
-          traits::mcc_real_or_char_range RT,
+double mcc_time_to_alt_limit(traits::mcc_real_or_char_range auto const& alt_limit,
-          traits::mcc_real_or_char_range DT,
+                             traits::mcc_real_or_char_range auto const& RA,
-          traits::mcc_real_or_char_range LT>
+                             traits::mcc_real_or_char_range auto const& DEC,
-double mcc_time_to_alt_limit(const AT& alt_limit,
+                             traits::mcc_real_or_char_range auto const& LAT,
-                             const RT& RA,
+                             traits::mcc_real_or_char_range auto const& LON,
-                             const DT& DEC,
+                             const std::chrono::system_clock::time_point& now,
-                             const LT& LAT,
+                             traits::mcc_time_duration_c auto const& dut1,    // UT1-UTC
-                             const std::chrono::system_clock::time_point& now)
+                             traits::mcc_time_duration_c auto const& tt_tai,  // TT-TAI
+                             // TAI-UTC (leap seconds)
+                             traits::mcc_time_duration_c auto const& tai_utc)
 {
     //     sin(alt) = sin(DEC)*sin(phi) + cos(DEC)*cos(phi)*cos(HA)
     //     HA = LST - RA
+    //     cos(HA) = cos(LST)*cos(RA) + sin(LST)*sin(RA)
 
-    double ra, dec, lat, alt;
+    using AT = std::decay_t<decltype(alt_limit)>;
+    using RT = std::decay_t<decltype(RA)>;
+    using DT = std::decay_t<decltype(DEC)>;
+    using LT = std::decay_t<decltype(LAT)>;
+    using LGT = std::decay_t<decltype(LON)>;
+
+    double ra, dec, lat, lon, alt;
 
     if constexpr (std::floating_point<AT>) {
         alt = alt_limit * utils::deg2radCoeff;
@@ -183,25 +207,53 @@ double mcc_time_to_alt_limit(const AT& alt_limit,
         lat *= utils::deg2radCoeff;
     }
 
+    if constexpr (std::floating_point<LGT>) {
+        lon = LON * utils::deg2radCoeff;
+    } else {
+        lon = utils::parsAngleString(LON);
+        lon *= utils::deg2radCoeff;
+    }
+
     if (lat >= 0.0) {                                // north hemisphere
-        if (dec < (lat - std::numbers::pi / 2.0)) {  // never rises above horizont
+        if (dec < (lat - std::numbers::pi / 2.0)) {  // never rises above horizon
             return std::numeric_limits<double>::quiet_NaN();
         }
     } else {                                         // south hemisphere
-        if (dec > (lat + std::numbers::pi / 2.0)) {  // never rises above horizont
+        if (dec > (lat + std::numbers::pi / 2.0)) {  // never rises above horizon
             return std::numeric_limits<double>::quiet_NaN();
         }
     }
 
     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)
+    if (std::abs(cos_ha) > 1.0) {  // it never sets below given altitude
         return std::numeric_limits<double>::infinity();
     }
 
     double lst = std::acos(cos_ha) + ra;
 
+    auto ut1 = now + dut1;
+    auto tt = now + tai_utc + tt_tai;
 
-    return 0.0;
+    double ut1_mjd, tt_mjd, result;
+    int ret = mcc_julday(ut1, ut1_mjd);
+    if (ret) {
+        return std::numeric_limits<double>::quiet_NaN();
+    }
+    ret = mcc_julday(tt, tt_mjd);
+    if (ret) {
+        return std::numeric_limits<double>::quiet_NaN();
+    }
+
+    double lst_now = erfa::eraGst06a(ERFA_DJM0, ut1_mjd, ERFA_DJM0, tt_mjd) + lon;
+
+    result = lst - lst_now;
+    if (result < 0.0) {  // the next sideral day
+        result += 2.0 * std::numbers::pi;
+    }
+
+    result *= mcc_UT1_to_sideral_ratio;  // to UT1 scale
+
+    return result;
 }
 
 
@@ -218,6 +270,7 @@ class MccLeapSeconds final
 {
 public:
     typedef std::chrono::system_clock::time_point time_point_t;
+    typedef std::chrono::duration<double> real_secs_t;  // seconds duration in double
 
     MccLeapSeconds()
     {
@@ -308,7 +361,8 @@ public:
     }
 
 
-    std::optional<double> operator[](const time_point_t& tp) const
+    // std::optional<double> operator[](const time_point_t& tp) const
+    std::optional<real_secs_t> operator[](const time_point_t& tp) const
     {
         if (tp > _expireDate) {  // ???????!!!!!!!!!!!
             return std::nullopt;
@@ -319,14 +373,16 @@ public:
 
         for (auto const& el : _db | std::views::reverse) {
             if (ymd >= el.ymd) {
-                return el.tai_utc;
+                // return el.tai_utc;
+                return real_secs_t{el.tai_utc};
             }
         }
 
         return std::nullopt;
     }
 
-    std::optional<double> operator[](const double& mjd) const
+    // std::optional<double> operator[](const double& mjd) const
+    std::optional<real_secs_t> operator[](const double& mjd) const
     {
         double e_mjd;
 
@@ -338,7 +394,7 @@ public:
 
         for (auto const& el : _db | std::views::reverse) {
             if (mjd >= el.mjd) {
-                return el.tai_utc;
+                return real_secs_t{el.tai_utc};
             }
         }
 
@@ -379,6 +435,7 @@ class MccIersBulletinA final
 {
 public:
     typedef std::chrono::system_clock::time_point time_point_t;
+    typedef std::chrono::duration<double> real_secs_t;  // seconds duration in double
 
     struct pole_pos_t {
         double x, y;
@@ -421,13 +478,15 @@ public:
     }
 
 
-    double TT_TAI() const
+    // double TT_TAI() const
+    real_secs_t TT_TAI() const
     {
-        return _tt_tai;
+        return real_secs_t{_tt_tai};
     }
 
     // DUT1 = UT1 - UTC
-    std::optional<double> DUT1(const time_point_t& tp) const
+    // std::optional<double> DUT1(const time_point_t& tp) const
+    std::optional<real_secs_t> DUT1(const time_point_t& tp) const
     {
         // use of the closest date
         std::chrono::year_month_day ymd{std::chrono::round<std::chrono::days>(tp)};
@@ -443,14 +502,15 @@ public:
         auto el_prev = _db.front();
         for (auto const& el : _db) {
             if (ymd <= el.ymd) {
-                return el.dut1;
+                return real_secs_t{el.dut1};
             }
         }
 
         return std::nullopt;
     }
 
-    std::optional<double> DUT1(double mjd) const
+    // std::optional<double> DUT1(double mjd) const
+    std::optional<real_secs_t> DUT1(double mjd) const
     {
         mjd = std::round(mjd);  // round to closest integer MJD
 
@@ -464,7 +524,7 @@ public:
 
         for (auto const& el : _db) {
             if (mjd <= el.mjd) {
-                return el.dut1;
+                return real_secs_t{el.dut1};
             }
         }
 

cxx/tests/astrom_test.cpp

@@ -4,6 +4,8 @@
 #include "../mount_astrom.h"
 
 
+// BTA coords from maps.yandex.ru: 43.646711, 41.440732
+
 int main(int argc, char* argv[])
 {
     int ecode = 0;
@@ -37,16 +39,21 @@ int main(int argc, char* argv[])
 
     std::cout << "\n\n\n------------------\n";
 
-    auto now = std::chrono::system_clock::now();
+    constexpr auto nnn = mcc::astrom::MccLeapSeconds::real_secs_t{std::numeric_limits<double>::quiet_NaN()};
-    std::cout << "LEAP SECS for now: " << lps[now].value_or(std::numeric_limits<double>::quiet_NaN()) << "\n";
 
-    std::cout << "DUT1 for now: " << bullA.DUT1(now).value_or(std::numeric_limits<double>::quiet_NaN()) << " (" << now
+    auto now = std::chrono::system_clock::now();
-              << ")\n";
+    // std::cout << "LEAP SECS for now: " << lps[now].value_or(std::numeric_limits<double>::quiet_NaN()) << "\n";
+    std::cout << "LEAP SECS for now: " << lps[now].value_or(nnn) << "\n";
+
+    // std::cout << "DUT1 for now: " << bullA.DUT1(now).value_or(std::numeric_limits<double>::quiet_NaN()) << " (" <<
+    // now
+    std::cout << "DUT1 for now: " << bullA.DUT1(now).value_or(nnn) << " (" << now << ")\n";
 
     double mjd = 61077.4;
 
-    std::cout << "DUT1 for MJD = " << mjd << ": " << bullA.DUT1(mjd).value_or(std::numeric_limits<double>::quiet_NaN())
+    // std::cout << "DUT1 for MJD = " << mjd << ": " <<
-              << "\n";
+    // bullA.DUT1(mjd).value_or(std::numeric_limits<double>::quiet_NaN())
+    std::cout << "DUT1 for MJD = " << mjd << ": " << bullA.DUT1(mjd).value_or(nnn) << "\n";
 
 
     auto st = mcc::astrom::mcc_julday(now, mjd);