...

cxx/mount_astrom.h

@@ -8,6 +8,7 @@
 
 #include <chrono>
 #include <fstream>
+#include "mcc_coord.h"
 
 #ifdef VEC_XSIMD
 #include <xsimd/xsimd.hpp>
@@ -160,6 +161,88 @@ static int mcc_julday(traits::mcc_systime_c auto const& start_time,
 }
 
 
+/*
+ *    Computes a time duration to reach given altitude
+ *
+ *    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
+ */
+std::pair<std::chrono::duration<double>, std::chrono::duration<double>> mcc_time_to_alt(
+    const MccCoordinate& alt_limit,
+    const MccCoordinate& ra,
+    const MccCoordinate& dec,
+    const MccCoordinate& lat,
+    const MccCoordinate& 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());
+
+    if (alt_limit < 0.0) {
+        return {nan_dur, nan_dur};
+    }
+
+    if (lat >= 0.0) {                                // north hemisphere
+        if (dec < (lat - std::numbers::pi / 2.0)) {  // never rises above horizon
+            return {nan_dur, nan_dur};
+        }
+    } else {                                         // south hemisphere
+        if (dec > (lat + std::numbers::pi / 2.0)) {  // never rises above horizon
+            return {nan_dur, nan_dur};
+        }
+    }
+
+    double cos_ha = (std::sin(alt_limit) - std::sin(dec) * std::sin(lat)) / std::cos(dec) / std::cos(lat);
+    if (std::abs(cos_ha) > 1.0) {  // it never reach given altitude
+        return {inf_dur, inf_dur};
+    }
+
+
+    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;
+
+    auto acos_ha = std::acos(cos_ha);
+
+    // before upper culmination
+    double lst_before = -acos_ha + ra;
+    // after upper culmination
+    double lst_after = acos_ha + ra;
+
+    double d1 = (lst_before - lst_now) * mcc_UT1_to_sideral_ratio,
+           d2 = (lst_after - lst_now) * mcc_UT1_to_sideral_ratio;
+
+    if (d1 < 0.0) {  // the next day
+        d1 += 2.0 * std::numbers::pi;
+    }
+    if (d2 < 0.0) {  // the next day
+        d2 += 2.0 * std::numbers::pi;
+    }
+
+    static constexpr double rad2secs = 12.0 / std::numbers::pi * 3600.0;  // radians to time seconds
+
+    return {std::chrono::duration<double>(d1 * rad2secs), std::chrono::duration<double>(d2 * rad2secs)};
+
+    // return std::chrono::duration<double>(result * 12.0 / std::numbers::pi * 3600.0);  // in seconds
+}
+
 /*
  *  angles are in degrees or sexagimal string form
  *  RA and DEC are apparent!
@@ -168,6 +251,7 @@ static int mcc_julday(traits::mcc_systime_c auto const& start_time,
  *      NaN if object is non-rising or "alt_limit" < 0, Inf is circumpolar
  */
 
+/*
 double mcc_time_to_alt_limit(traits::mcc_real_or_char_range auto const& alt_limit,
                              traits::mcc_real_or_char_range auto const& RA,
                              traits::mcc_real_or_char_range auto const& DEC,
@@ -284,7 +368,7 @@ double mcc_time_to_alt_limit(traits::mcc_real_or_char_range auto const& alt_limi
     return result;
 }
 
-
+*/
 
 /* Classes to represent IERS bulletins
  *