From 18628aff281e7518404017012ce3c85af0a13210 Mon Sep 17 00:00:00 2001
From: "Timur A. Fatkhullin" <amirych@gmail.com>
Date: Tue, 29 Apr 2025 00:11:51 +0300
Subject: [PATCH] ...

---
 cxx/mcc_coord.h    |  6 ++--
 cxx/mount_astrom.h | 75 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 79 insertions(+), 2 deletions(-)

diff --git a/cxx/mcc_coord.h b/cxx/mcc_coord.h
index 6a496bc..c589f3f 100644
--- a/cxx/mcc_coord.h
+++ b/cxx/mcc_coord.h
@@ -260,7 +260,8 @@ public:
 
     auto operator==(const MccAngle& other) const
     {
-        return _angleInRads == other._angleInRads;
+        return utils::isEqual(_angleInRads, other._angleInRads);
+        // return _angleInRads == other._angleInRads;
     }
 
     auto operator<=>(double other) const
@@ -270,7 +271,8 @@ public:
 
     auto operator==(double other) const
     {
-        return _angleInRads == other;
+        return utils::isEqual(_angleInRads, other);
+        // return _angleInRads == other;
     }
 };
 
diff --git a/cxx/mount_astrom.h b/cxx/mount_astrom.h
index 3119e5b..8a63c77 100644
--- a/cxx/mount_astrom.h
+++ b/cxx/mount_astrom.h
@@ -167,6 +167,9 @@ static int mcc_julday(traits::mcc_systime_c auto const& start_time,
  *    the function returns a std::pair:
  *        .first = time duration to reach given altitude BEFORE object upper culmination
  *        .second = time duration to reach given altitude AFTER object upper culmination
+ *
+ *
+ *    input RA and DEC are apparent!
  */
 std::pair<std::chrono::duration<double>, std::chrono::duration<double>> mcc_time_to_alt(
     const MccAngle& alt,
@@ -370,6 +373,78 @@ double mcc_time_to_alt_limit(traits::mcc_real_or_char_range auto const& alt_limi
 
 */
 
+
+/*
+ *    Computes a time duration to reach given azimuth
+ *
+ *    the function returns a std::pair:
+ *        .first = time duration to reach given altitude BEFORE object upper culmination
+ *        .second = time duration to reach given altitude AFTER object upper culmination
+ *
+ *
+ *    input RA and DEC are apparent!
+ *    input azimuth is assumed to be started from the South through the West
+ *    (the South is az = 0, the West is az = pi/2 ...)
+ */
+
+std::pair<std::chrono::duration<double>, std::chrono::duration<double>> mcc_time_to_az(
+    const MccAngle& az,
+    const MccAngle& ra,
+    const MccAngle& dec,
+    const MccAngle& lat,
+    const MccAngle& lon,
+    traits::mcc_systime_c auto const& now,
+    traits::mcc_time_duration_c auto const& dut1,    // UT1-UTC
+    traits::mcc_time_duration_c auto const& tt_tai,  // TT-TAI
+    // TAI-UTC (leap seconds)
+    traits::mcc_time_duration_c auto const& tai_utc)
+{
+    auto nan_dur = std::chrono::duration<double>(std::numeric_limits<double>::quiet_NaN());
+    auto inf_dur = std::chrono::duration<double>(std::numeric_limits<double>::infinity());
+
+
+    auto cs_d = cos(dec);
+    auto cos_az = cos(az);
+
+    auto C1 = cs_d / tan(az);
+    auto C2 = -sin(dec) * cos(lat);
+    auto C3 = cs_d * sin(lat);
+
+    auto C2C3 = C2 * C3;
+    auto C1_2 = C1 * C1;
+    auto C2_2 = C2 * C2;
+    auto C3_2 = C2 * C3;
+
+    auto C1C3 = C1_2 + C3_2;
+    auto D = C2_2 * C3_2 - C1C3 * (C2_2 - C1_2);
+    if (D < 0.0) {
+        return {inf_dur, inf_dur};
+    }
+
+    auto cos_t1 = -(C2C3 + sqrt(D)) / C1C3;
+    auto cos_t2 = -(C2C3 - sqrt(D)) / C1C3;
+
+    auto sin_z1 = (C2 + C3 * cos_t1) / cos_az;
+    auto sin_z2 = (C2 + C3 * cos_t2) / cos_az;
+
+    auto ut1 = now + dut1;
+    auto tt = now + tai_utc + tt_tai;
+
+    double ut1_mjd, tt_mjd;
+    int ret = mcc_julday(ut1, ut1_mjd);
+    if (ret) {
+        return {nan_dur, nan_dur};
+    }
+    ret = mcc_julday(tt, tt_mjd);
+    if (ret) {
+        return {nan_dur, nan_dur};
+    }
+
+    double lst_now = erfa::eraGst06a(ERFA_DJM0, ut1_mjd, ERFA_DJM0, tt_mjd);
+    lst_now += lon;
+}
+
+
 /* Classes to represent IERS bulletins
  *
  * BULLETIN A: https://datacenter.iers.org/data/latestVersion/bulletinA.txt