From 0aa0113be35152d6afa6d7fa4662f9f1e69c7e70 Mon Sep 17 00:00:00 2001 From: "Timur A. Fatkhullin" Date: Thu, 15 Jan 2026 23:21:28 +0300 Subject: [PATCH] ... --- mcc/mcc_ccte_erfa_new.h | 53 ++++++++++++++++- mcc/mcc_coord.h | 110 +++++++++++++++++++++++------------ mcc/tests/mcc_coord_test.cpp | 18 ++++-- 3 files changed, 140 insertions(+), 41 deletions(-) diff --git a/mcc/mcc_ccte_erfa_new.h b/mcc/mcc_ccte_erfa_new.h index e2fe512..fd0afad 100644 --- a/mcc/mcc_ccte_erfa_new.h +++ b/mcc/mcc_ccte_erfa_new.h @@ -250,6 +250,51 @@ public: } + // apparent sideral time (Greenwitch or local) + error_t apparentSideralTime(mcc_coord_epoch_c auto const& epoch, mcc_angle_c auto* st, bool islocal = false) + { + error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK; + + if (st == nullptr) { + return MccCCTE_ERFAErrorCode::ERROR_NULLPTR; + } + + using real_days_t = std::chrono::duration>; + + double ut1 = epoch.MJD(); + double tt = epoch.MJD(); + + std::lock_guard lock{*_stateMutex}; + + auto dut1 = _currentState._bulletinA.DUT1(epoch.MJD()); + + if (dut1.has_value()) { + ut1 += std::chrono::duration_cast(dut1.value()).count(); + } else { // out of range + return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; + } + + auto tai_utc = _currentState._leapSeconds[epoch.MJD()]; + if (tai_utc.has_value()) { + tt += std::chrono::duration_cast(tai_utc.value()).count(); + } else { + return MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE; + } + + + auto tt_tai = _currentState._bulletinA.TT_TAI(); + tt += std::chrono::duration_cast(tt_tai).count(); + + *st = eraGst06a(ERFA_DJM0, ut1, ERFA_DJM0, tt); + + if (islocal) { + *st = eraAnp(*st + _currentState.lon); + } + + return ret; + } + + // ICRS to observed // returned azimuth is counted from the South through the West error_t icrsToObs(mcc_angle_c auto const& ra_icrs, @@ -401,7 +446,7 @@ public: } - error_t equationOrigins(const mcc_julday_c auto& mjd, mcc_angle_c auto* eo) + error_t equationOrigins(mcc_coord_epoch_c auto const& epoch, mcc_angle_c auto* eo) { if (eo == nullptr) { return MccCCTE_ERFAErrorCode::ERROR_NULLPTR; @@ -413,6 +458,8 @@ public: using real_days_t = std::chrono::duration>; + double mjd = epoch.MJD(); + auto tai_utc = _currentState._leapSeconds[mjd]; if (tai_utc.has_value()) { double tt = mjd; @@ -658,6 +705,10 @@ protected: *dec = d; } + if (ha) { + *ha = h; + } + if (az) { // NOTE: according to definition of astronomical azimuth it is counted from the South through the West, but // in the ERFA the azimuth is counted from the North through the East!!! diff --git a/mcc/mcc_coord.h b/mcc/mcc_coord.h index 99afc60..9a2e1d6 100644 --- a/mcc/mcc_coord.h +++ b/mcc/mcc_coord.h @@ -31,13 +31,13 @@ public: ? MccCoordPairKind::COORDS_KIND_RADEC_APP // observed RA and DEC : (std::same_as && std::same_as) - ? MccCoordPairKind::COORDS_KIND_RADEC_APP + ? MccCoordPairKind::COORDS_KIND_RADEC_OBS // apparent HA and DEC : (std::same_as && std::same_as) ? MccCoordPairKind::COORDS_KIND_HADEC_APP // observed HA and DEC : (std::same_as && std::same_as) - ? MccCoordPairKind::COORDS_KIND_HADEC_APP + ? MccCoordPairKind::COORDS_KIND_HADEC_OBS // apparent AZ and ZD : (std::same_as && std::same_as) ? MccCoordPairKind::COORDS_KIND_AZZD @@ -132,6 +132,11 @@ template class MccNamedCoordPair : public MccCoordPair { public: + MccNamedCoordPair() : MccCoordPair(CO_LON_T{0.0}, CO_LAT_T{0.0}, MccCelestialCoordEpoch::now()) + { + } + + template requires(std::is_arithmetic_v && std::is_arithmetic_v) MccNamedCoordPair(CxT const& x, CyT const& y, EpT const& epoch = EpT::now()) @@ -166,7 +171,7 @@ struct MccSkyRADEC_ICRS : MccNamedCoordPair { MccSkyRADEC_ICRS(MccAngle const& x, MccAngle const& y) : MccSkyRADEC_ICRS((double)x, (double)y) {} // ignore epoch setting (it is always J2000.0) - void setEpoch(mcc_coord_epoch_c auto const& ep) + void setEpoch(mcc_coord_epoch_c auto const&) { static_assert(false, "CANNOT SET EPOCH FOR ICRS-KIND COORDINATE PAIR!!!"); } @@ -403,33 +408,33 @@ protected: return; } else { // from ICRS to apparent or observed - if constexpr (mccIsAppCoordPairKind) { - ccte_err = cctEngine.icrsToApp(_x, _y, cpair.epoch(), &ra, &dec, &ha, &az, &zd); - } else if constexpr (mccIsObsCoordPairKind) { - ccte_err = cctEngine.icrsToObs(_x, _y, cpair.epoch(), &ra, &dec, &ha, &az, &zd); - } else { - static_assert(true, "UNSUPPORTED SKY POINT TRANSFORMATION!"); - } + // if constexpr (mccIsAppCoordPairKind) { + // ccte_err = cctEngine.icrsToApp(_x, _y, cpair.epoch(), &ra, &dec, &ha, &az, &zd); + // } else if constexpr (mccIsObsCoordPairKind) { + // ccte_err = cctEngine.icrsToObs(_x, _y, cpair.epoch(), &ra, &dec, &ha, &az, &zd); + // } else { + // static_assert(true, "UNSUPPORTED SKY POINT TRANSFORMATION!"); + // } - if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP || - PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { - cpair.setX(ra); - cpair.setY(dec); - } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP || - PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { - cpair.setX(ha); - cpair.setY(dec); - } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) { - cpair.setX(az); - cpair.setY(zd); - } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) { - cpair.setX(az); - cpair.setY(half_pi - zd); - } else { - static_assert(true, "UNSUPPORTED SKY POINT TRANSFORMATION!"); - } + // if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP || + // PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { + // cpair.setX(ra); + // cpair.setY(dec); + // } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP || + // PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { + // cpair.setX(ha); + // cpair.setY(dec); + // } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) { + // cpair.setX(az); + // cpair.setY(zd); + // } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) { + // cpair.setX(az); + // cpair.setY(half_pi - zd); + // } else { + // static_assert(true, "UNSUPPORTED SKY POINT TRANSFORMATION!"); + // } - return; + // return; } } @@ -476,11 +481,17 @@ protected: // here, the input coordinates and stored one are at the same epoch - ccte_err = cctEngine.equationOrigins(cpair.MJD(), &eo); + ccte_err = cctEngine.equationOrigins(cpair.epoch(), &eo); if (ccte_err) { return; } + ccte_err = cctEngine.apparentSideralTime(cpair.epoch(), &lst, true); + if (ccte_err) { + return; + } + + if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_APP || pkind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { ra = _x; dec = _y; @@ -494,16 +505,43 @@ protected: } else if (pkind == MccCoordPairKind::COORDS_KIND_AZALT) { az = _x; alt = _y; - } else if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { - ra_icrs = _x; - dec_icrs = _y; - } else { // unsupported transformation!!! - return; } + // else if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { + // ra_icrs = _x; + // dec_icrs = _y; + // } else { // unsupported transformation!!! + // return; + // } // coordinate transformation lambda (possibly recursive!!!) - auto comp_func = [&, this](this auto&& self, MccCoordPairKind cp_kind) -> void { - if (cp_kind == MccCoordPairKind::COORDS_KIND_AZALT) { + auto comp_func = [&](this auto&& self, MccCoordPairKind cp_kind) -> void { + if (cp_kind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { + if constexpr (mccIsAppCoordPairKind) { + ccte_err = cctEngine.icrsToApp(ra_icrs, dec_icrs, cpair.epoch(), &ra, &dec, &ha, &az, &zd); + } else if constexpr (mccIsObsCoordPairKind) { + ccte_err = cctEngine.icrsToObs(ra_icrs, dec_icrs, cpair.epoch(), &ra, &dec, &ha, &az, &zd); + } else { + static_assert(true, "UNSUPPORTED SKY POINT TRANSFORMATION!"); + } + + if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP || + PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { + cpair.setX(ra); + cpair.setY(dec); + } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP || + PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { + cpair.setX(ha); + cpair.setY(dec); + } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) { + cpair.setX(az); + cpair.setY(zd); + } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) { + cpair.setX(az); + cpair.setY(half_pi - zd); + } else { + static_assert(true, "UNSUPPORTED SKY POINT TRANSFORMATION!"); + } + } else if (cp_kind == MccCoordPairKind::COORDS_KIND_AZALT) { if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) { zd = half_pi - alt; cpair.setX(az); diff --git a/mcc/tests/mcc_coord_test.cpp b/mcc/tests/mcc_coord_test.cpp index 3dbf83c..0c79d4e 100644 --- a/mcc/tests/mcc_coord_test.cpp +++ b/mcc/tests/mcc_coord_test.cpp @@ -23,15 +23,23 @@ int main() skypt_t::cctEngine.setStateERFA(saoras); skypt_t pt; - MccSkyRADEC_ICRS icrs(0.0, 0.0); + MccSkyRADEC_ICRS icrs(0.0, 30.0_degs); pt = icrs; MccSkyRADEC_OBS radec_obs{0.0, 0.0}; + MccSkyRADEC_APP radec_app; MccSkyAZALT azalt{0, 0}; MccSkyAZZD azzd{0, 0}; + MccSkyHADEC_OBS hadec_obs; + MccAngle eo, lst; - pt.to(radec_obs, azalt, azzd); + skypt_t::cctEngine.equationOrigins(radec_obs.epoch(), &eo); + skypt_t::cctEngine.apparentSideralTime(radec_obs.epoch(), &lst, true); + std::cout << "EO = " << eo.sexagesimal(true) << "\n"; + std::cout << "LST = " << lst.sexagesimal(true) << "\n\n"; + + pt.to(radec_obs, azalt, azzd, radec_app); std::cout << "FROM ICRS TO OBSERVED:\n"; std::cout << "RA_ICRS = " << icrs.x().sexagesimal(true) << "\n"; @@ -42,11 +50,13 @@ int main() std::cout << "AZ = " << azalt.x().sexagesimal() << "\n"; std::cout << "ALT = " << azalt.y().sexagesimal() << "\n"; std::cout << "ZD = " << azzd.y().sexagesimal() << "\n"; + std::cout << "RA_APP = " << radec_app.x().sexagesimal(true) << "\n"; + std::cout << "DEC_APP = " << radec_app.y().sexagesimal() << "\n"; // radec_obs = {10.2387983_degs, "43:21:34.5465"_dms}; pt = radec_obs; - pt.to(icrs); + pt.to(icrs, azzd, hadec_obs); std::cout << "\n\nFROM OBSERVED TO ICRS:\n"; std::cout << "OBS COORD EPOCH: " << radec_obs.epoch().UTC() << "\n"; @@ -55,11 +65,11 @@ int main() std::cout << "RA_ICRS = " << icrs.x().sexagesimal(true) << "\n"; std::cout << "DEC_ICRS = " << icrs.y().sexagesimal() << "\n"; - pt.to(azzd); std::cout << "\n\nFROM OBSERVED TO OBSERVED:\n"; std::cout << "OBS COORD EPOCH: " << radec_obs.epoch().UTC() << "\n"; std::cout << "RA_OBS = " << radec_obs.x().sexagesimal(true) << "\n"; std::cout << "DEC_OBS = " << radec_obs.y().sexagesimal() << "\n"; + std::cout << "HA_OBS = " << hadec_obs.x().sexagesimal(true) << "\n"; std::cout << "AZ = " << azzd.x().sexagesimal() << "\n"; std::cout << "ZD = " << azzd.y().sexagesimal() << "\n";