From 67340d1926af0f085eb51466aca55412babeea14 Mon Sep 17 00:00:00 2001 From: "Timur A. Fatkhullin" Date: Mon, 21 Apr 2025 18:36:47 +0300 Subject: [PATCH] ... --- cxx/mcc_coord.h | 84 +++++++++++++++++++++++++++++++++++++- cxx/mount_astrom.h | 86 ++++++++++++++++++++++++++++++++++++++- cxx/mount_pz.h | 37 +++++++++++++++++ cxx/tests/astrom_test.cpp | 24 +++++++++-- 4 files changed, 225 insertions(+), 6 deletions(-) diff --git a/cxx/mcc_coord.h b/cxx/mcc_coord.h index b54b4b4..ddc14fd 100644 --- a/cxx/mcc_coord.h +++ b/cxx/mcc_coord.h @@ -57,10 +57,17 @@ static constexpr MccHMSTag mcc_hms{}; class MccCoordinate { +protected: double _angleInRads{0.0}; int _precision{2}; public: + enum norm_kind_t { + NORM_KIND_0_360, // [0,360] + NORM_KIND_0_180, // [0,180] + NORM_KIND_90_90, // [-90,90] + }; + MccCoordinate() = default; // by default angle is in radians @@ -102,6 +109,45 @@ public: } } + virtual ~MccCoordinate() = default; + + // normalize coordinate + template + MccCoordinate& normalize() + { + _angleInRads = std::fmod(_angleInRads, std::numbers::pi * 2.0); + + if constexpr (KIND == NORM_KIND_0_360) { + if (_angleInRads < 0.0) { + _angleInRads += 2.0 * std::numbers::pi; + } + } else if constexpr (KIND == NORM_KIND_0_180) { + if (_angleInRads < -std::numbers::pi) { + _angleInRads = 2.0 * std::numbers::pi + _angleInRads; + } else if (_angleInRads < 0.0) { + _angleInRads = -_angleInRads; + } + } else if constexpr (KIND == NORM_KIND_90_90) { + if (_angleInRads >= 1.5 * std::numbers::pi) { + _angleInRads = _angleInRads - 2.0 * std::numbers::pi; + } else if (_angleInRads >= std::numbers::pi / 2.0) { + _angleInRads = std::numbers::pi - _angleInRads; + } else if (_angleInRads <= -1.5 * std::numbers::pi) { + _angleInRads += 2.0 * std::numbers::pi; + } else if (_angleInRads <= -std::numbers::pi / 2.0) { + _angleInRads = -(std::numbers::pi + _angleInRads); + } + } + + return *this; + } + + MccCoordinate& normalize() + { + return normalize(); + } + + operator double() const { return _angleInRads; @@ -132,7 +178,10 @@ public: friend MccCoordinate operator+(std::convertible_to auto&& left, std::convertible_to auto&& right) { - return std::forward(left)._angleInRads + std::forward(right)._angleInRads; + // return std::forward(left)._angleInRads + + // std::forward(right)._angleInRads; + return MccCoordinate{std::forward(left)}._angleInRads + + MccCoordinate{std::forward(right)}._angleInRads; } // friend MccCoordinate operator+(const MccCoordinate& left, const MccCoordinate& right) // { @@ -143,7 +192,10 @@ public: friend MccCoordinate operator-(std::convertible_to auto&& left, std::convertible_to auto&& right) { - return std::forward(left)._angleInRads - std::forward(right)._angleInRads; + // return std::forward(left)._angleInRads - + // std::forward(right)._angleInRads; + return MccCoordinate{std::forward(left)}._angleInRads - + MccCoordinate{std::forward(right)}._angleInRads; } // friend MccCoordinate operator-(const MccCoordinate& left, const MccCoordinate& right) // { @@ -220,5 +272,33 @@ public: } }; +enum class MccNormCoordKind { + NORM_KIND_0_360, // [0,360] + NORM_KIND_0_180, // [0,180] + NORM_KIND_90_90, // [-90,90] +}; + +template +class MccNormCoordinate : public MccCoordinate +{ +public: + static constexpr MccNormCoordKind normKind = KIND; + +protected: + static constexpr double pi2 = 2.0 * std::numbers::pi; + static constexpr double normScale = KIND == MccNormCoordKind::NORM_KIND_0_360 ? std::numbers::pi * 2.0 + : KIND == MccNormCoordKind::NORM_KIND_0_180 ? std::numbers::pi + : std::numbers::pi / 2.0; + + void norm() + { + _angleInRads = std::fmod(_angleInRads, normScale); + if constexpr (KIND == MccNormCoordKind::NORM_KIND_0_360) { + _angleInRads += normScale; + } else if constexpr (KIND == MccNormCoordKind::NORM_KIND_0_180) { + _angleInRads = -_angleInRads; + } + } +}; } // namespace mcc diff --git a/cxx/mount_astrom.h b/cxx/mount_astrom.h index 3b9772a..d8cc8f9 100644 --- a/cxx/mount_astrom.h +++ b/cxx/mount_astrom.h @@ -8,6 +8,7 @@ #include #include +#include "mcc_coord.h" #ifdef VEC_XSIMD #include @@ -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> 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(std::numeric_limits::quiet_NaN()); + auto inf_dur = std::chrono::duration(std::numeric_limits::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(d1 * rad2secs), std::chrono::duration(d2 * rad2secs)}; + + // return std::chrono::duration(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 * diff --git a/cxx/mount_pz.h b/cxx/mount_pz.h index 3dc0ad2..58cbe68 100644 --- a/cxx/mount_pz.h +++ b/cxx/mount_pz.h @@ -4,6 +4,43 @@ namespace mcc { + +namespace traits +{ + +template +concept mcc_prohibited_zone_c = requires(T t, const T const_t) { + typename T::pzcoords_kind_t; + + { const_t.name() } -> std::same_as; + { const_t.desc() } -> std::same_as; + + // check if given coordinates are within the zone + { + t.inZone(std::declval(), std::declval(), + std::declval()) + } -> std::convertible_to; + + // a time duration to reach the zone (0 - if already in the zone, chrono::duration<>::max() if never + // reach the zone) + { + t.timeTo(std::declval(), std::declval(), + std::declval(), + std::declval()) + } -> mcc_time_duration_c; + + // a time duration to exit the zone (0 - if already out of the zone, chrono::duration<>::max() if + // never exit the zone) + { + t.timeFrom(std::declval(), std::declval(), + std::declval(), + std::declval()) + } -> mcc_time_duration_c; +}; + +} // namespace traits + + class MccProhibitedZone { public: diff --git a/cxx/tests/astrom_test.cpp b/cxx/tests/astrom_test.cpp index 295f6f2..814fa7e 100644 --- a/cxx/tests/astrom_test.cpp +++ b/cxx/tests/astrom_test.cpp @@ -74,11 +74,29 @@ int main(int argc, char* argv[]) std::cout << "\n\nTimes to object (RA = " << ra_str << ", DEC = " << dec_str << ") sets to given altitude (" << alt_lim << " degrees):\n"; using namespace std::chrono_literals; - double stm = mcc::astrom::mcc_time_to_alt_limit(alt_lim, ra_str, dec_str, 43.646711, 41.440732, - std::chrono::system_clock::now(), 0.041s, 32.184s, 37.0s); + // auto stm = mcc::astrom::mcc_time_to_alt_limit(alt_lim, ra_str, dec_str, 43.646711, 41.440732, + // std::chrono::system_clock::now(), 0.041s, 32.184s, 37.0s); // auto stm_d = mcc::astrom::mcc_chrono_radians{stm}; - std::cout << "STM: " << stm * 12.0 / std::numbers::pi * 60.0 << " minutes\n"; + // std::cout << "STM: " << stm * 12.0 / std::numbers::pi * 60.0 << " minutes\n"; + // std::cout << "STM: " << std::chrono::duration_cast(stm) << " minutes\n"; + + alt_lim = 85.0; + ra_str = "02:30:00.0", dec_str = "45:00:00.0"; + std::cout << "\n\nTimes to object (RA = " << ra_str << ", DEC = " << dec_str << ") sets to given altitude (" + << alt_lim << " degrees):\n"; + + auto stm = + mcc::astrom::mcc_time_to_alt({alt_lim, mcc::mcc_degrees}, {ra_str, mcc::mcc_hms}, dec_str, 43.646711_degs, + 41.440732_degs, std::chrono::system_clock::now(), 0.041s, 32.184s, 37.0s); + // stm = mcc::astrom::mcc_time_to_alt_limit(alt_lim, ra_str, dec_str, 43.646711, 41.440732, + // std::chrono::system_clock::now(), 0.041s, 32.184s, 37.0s); + + // auto stm_d = mcc::astrom::mcc_chrono_radians{stm}; + // std::cout << "STM: " << stm * 12.0 / std::numbers::pi * 60.0 << " minutes\n"; + std::cout << "STM: " << stm.first << ", " << stm.second << " seconds\n"; + std::cout << "STM: " << std::chrono::duration_cast(stm.first) << ", " + << std::chrono::duration_cast(stm.second) << " minutes\n"; std::cout << "\n\n\n";