From ca74d1dd73dc0e106bed41449297669b111778ad Mon Sep 17 00:00:00 2001 From: "Timur A. Fatkhullin" Date: Mon, 18 Aug 2025 12:22:45 +0300 Subject: [PATCH] ... --- mcc/CMakeLists.txt | 2 +- mcc/mcc_angle.h | 532 ++++++++++++++++++++++++++++++++++++++++++++ mcc/mcc_ccte_erfa.h | 468 ++++++++++++++++++++++++++++++-------- mcc/mcc_generics.h | 43 ++-- mcc/mcc_pzone.h | 78 +++++++ 5 files changed, 1007 insertions(+), 116 deletions(-) create mode 100644 mcc/mcc_angle.h create mode 100644 mcc/mcc_pzone.h diff --git a/mcc/CMakeLists.txt b/mcc/CMakeLists.txt index 6279612..9ef4230 100644 --- a/mcc/CMakeLists.txt +++ b/mcc/CMakeLists.txt @@ -65,7 +65,7 @@ include_directories(${ERFA_INCLUDE_DIR}) -set(MCC_LIBRARY_SRC1 mcc_generics.h mcc_defaults.h mcc_traits.h mcc_utils.h mcc_ccte_iers.h mcc_ccte_iers_default.h) +set(MCC_LIBRARY_SRC1 mcc_generics.h mcc_defaults.h mcc_traits.h mcc_utils.h mcc_ccte_iers.h mcc_ccte_iers_default.h mcc_ccte_erfa.h) set(MCC_LIBRARY1 mcc1) add_library(${MCC_LIBRARY1} INTERFACE ${MCC_LIBRARY_SRC1}) target_compile_features(${MCC_LIBRARY1} INTERFACE cxx_std_23) diff --git a/mcc/mcc_angle.h b/mcc/mcc_angle.h new file mode 100644 index 0000000..ec253d0 --- /dev/null +++ b/mcc/mcc_angle.h @@ -0,0 +1,532 @@ +#pragma once + +#include "mcc_traits.h" +#include "mcc_utils.h" + +constexpr double operator""_rads(long double val) // angle in radians (no conversion) +{ + return val; +} + +constexpr double operator""_degs(long double val) // angle in degrees +{ + return val * std::numbers::pi / 180.0; +} + +constexpr double operator""_arcmins(long double val) // angle in arc minutes +{ + return val * std::numbers::pi / 180.0 / 60.0; +} + +constexpr double operator""_arcsecs(long double val) // angle in arc seconds +{ + return val * std::numbers::pi / 180.0 / 3600.0; +} + +constexpr double operator""_dms(const char* s, size_t size) // as a string "DEGREES:MINUTES:SECONDS" +{ + auto res = mcc::utils::parsAngleString(std::span{s, size}); + if (res.has_value()) { + return res.value() * std::numbers::pi / 180.0; + } else { + throw std::invalid_argument("invalid sexagesimal representation"); + } +} + +constexpr double operator""_hms(const char* s, size_t len) // as a string "HOURS:MINUTES:SECONDS" +{ + auto res = mcc::utils::parsAngleString(std::span{s, len}, true); + if (res.has_value()) { + return res.value() * std::numbers::pi / 180.0; + } else { + throw std::invalid_argument("invalid sexagesimal representation"); + } +} + + +namespace mcc +{ + +/* MCC-LIBRARY COORDINATES REPRESENTATION DEFINITIONS AND CLASS */ + + + +// tags for MccAngle class construction + +struct MccRadianTag { +}; +struct MccDegreeTag { +}; +static constexpr MccDegreeTag mcc_degrees{}; +struct MccHMSTag { +}; +static constexpr MccHMSTag mcc_hms{}; + + +class MccAngle +{ +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_180_180, // [-180,180] + NORM_KIND_90_90, // [-90,90] + }; + + MccAngle() = default; + + // by default angle is in radians + // MccAngle(std::convertible_to auto const& val, const MccRadianTag = MccRadianTag{}) : _angleInRads(val) {} + + // construct angle in degrees, e.g.: + // auto ang = MccAngle{180.0, mcc_degrees}; + // MccAngle(std::convertible_to auto const& val, const MccDegreeTag) + // { + // _angleInRads = val * utils::deg2radCoeff; + // } + + // by default angle is in radians + // template + // MccAngle(const T& val, const MccRadianTag = MccRadianTag{}) + // requires std::is_arithmetic_v + // : _angleInRads(val) + // { + // } + + // // construct angle in degrees, e.g.: + // // auto ang = MccAngle{180.0, mcc_degrees}; + // template + // MccAngle(const T& val, const MccDegreeTag) + // requires std::is_arithmetic_v + // : _angleInRads(val * utils::deg2radCoeff) + // { + // } + + + // by default angle is in radians + constexpr MccAngle(const double& val, const MccRadianTag = MccRadianTag{}) : _angleInRads(val) {} + + // construct angle in degrees, e.g.: + // auto ang = MccAngle{180.0, mcc_degrees}; + constexpr MccAngle(const double& val, const MccDegreeTag) : _angleInRads(val * utils::deg2radCoeff) {} + + + + // constuct angle from sexagesimal representation or floating-point number of degrees, e.g.: + // auto ang = MccAngle{"-12:34:56.789"}; // from degrees:minutes:seconds + // auto ang = MccAngle{"123.574698"}; // from degrees + constexpr MccAngle(traits::mcc_input_char_range auto const& val) + { + auto res = utils::parsAngleString(val); + if (res.has_value()) { + _angleInRads = res.value() * utils::deg2radCoeff; + } else { + throw std::invalid_argument("invalid sexagesimal representation"); + } + } + + // construct angle from sexagesimal representation or floating-point number of degrees, e.g.: + // auto ang = MccAngle{"01:23:45.6789", mcc_hms}; // from hours:minutes:seconds + // auto ang = MccAngle{"123.574698"}; // from degrees + constexpr MccAngle(traits::mcc_input_char_range auto const& val, const MccHMSTag) + { + auto res = utils::parsAngleString(val, true); + if (res.has_value()) { + _angleInRads = res.value() * utils::deg2radCoeff; + } else { + throw std::invalid_argument("invalid sexagesimal representation"); + } + } + + virtual ~MccAngle() = default; + + template T> + constexpr auto& operator=(this T&& self, const T& other) + { + std::forward(self)._angleInRads = other._angleInRads; + + return self; + } + + template T> + constexpr auto& operator=(this T&& self, T&& other) + { + std::forward(self)._angleInRads = std::move(other._angleInRads); + + return self; + } + + template + constexpr auto& operator=(this auto&& self, const T& val) + requires std::is_arithmetic_v + { + std::forward(self)._angleInRads = val; + + return self; + } + + // normalize coordinate + template + MccAngle& 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_180_180) { + if (_angleInRads > std::numbers::pi) { + _angleInRads = 2.0 * std::numbers::pi - _angleInRads; + } else if (_angleInRads < -std::numbers::pi) { + _angleInRads += 2.0 * std::numbers::pi; + } + } 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; + } + + MccAngle& normalize() + { + return normalize(); + } + + + // template + // operator T() const + // requires std::is_arithmetic_v + // { + // return _angleInRads; + // } + + constexpr operator double() const + { + return _angleInRads; + } + + + template + constexpr T degrees() const + { + return _angleInRads * 180.0 / std::numbers::pi; + } + + constexpr double degrees() const + { + return degrees(); + } + + + template + T arcmins() const + { + return _angleInRads * 10800.0 / std::numbers::pi; + } + + double arcmins() const + { + return arcmins(); + } + + + template + T arcsecs() const + { + return _angleInRads * 648000.0 / std::numbers::pi; + } + + double arcsecs() const + { + return arcsecs(); + } + + template + T hours() const + { + return _angleInRads * 12.0 / std::numbers::pi; + } + + double hours() const + { + return hours(); + } + + template + T minutes() const + { + return _angleInRads * 720.0 / std::numbers::pi; + } + + double minutes() const + { + return minutes(); + } + + template + T seconds() const + { + return _angleInRads * 43200.0 / std::numbers::pi; + } + + double seconds() const + { + return seconds(); + } + + + template + T sexagesimal(bool hms = false, int prec = 2) const + { + return utils::rad2sxg(_angleInRads, hms, prec >= 0 ? prec : _precision); + } + + std::string sexagesimal(bool hms = false, int prec = 2) const + { + return sexagesimal(hms, prec); + } + + + // arithmetics + + template T> + SelfT& operator+=(this SelfT& self, const T& v) + { + if constexpr (std::derived_from) { + static_assert(std::derived_from, "INCOMPATIBLE TYPES!"); + + self._angleInRads += v._angleInRads; + } else if constexpr (std::is_arithmetic_v) { + self._angleInRads += v; + } else { + self._angleInRads += MccAngle(v)._angleInRads; + } + + return self; + } + + template T> + SelfT& operator-=(this SelfT& self, const T& v) + { + if constexpr (std::derived_from) { + static_assert(std::derived_from, "INCOMPATIBLE TYPES!"); + + self._angleInRads -= v._angleInRads; + } else if constexpr (std::is_arithmetic_v) { + self._angleInRads -= v; + } else { + self._angleInRads -= MccAngle(v)._angleInRads; + } + + return self; + } + + template + SelfT& operator*=(this SelfT& self, const T& v) + requires std::is_arithmetic_v + { + self._angleInRads *= v; + + return self; + } + + template + SelfT& operator/=(this SelfT& self, const T& v) + requires std::is_arithmetic_v + { + self._angleInRads /= v; + + return self; + } + + + // unary '-' and '+' + template + SelfT operator-(this SelfT& self) + { + SelfT res = -self._angleInRads; + + return res; + } + + template + SelfT operator+(this SelfT& self) + { + return self; + } +}; + + +// binary arithmetic operations + +template T1, std::convertible_to T2> +static auto operator+(const T1& v1, const T2& v2) +{ + static_assert(std::convertible_to || std::convertible_to, "INCOMPATIBLE TYPES!"); + + using res_t = std::conditional_t && std::derived_from, T1, T2>; + + return res_t{(double)v1 + (double)v2}; +} + +template T1, std::convertible_to T2> +static auto operator-(const T1& v1, const T2& v2) +{ + static_assert(std::convertible_to || std::convertible_to, "INCOMPATIBLE TYPES!"); + + using res_t = std::conditional_t && std::derived_from, T1, T2>; + + return res_t{(double)v1 - (double)v2}; +} + + +template T1, std::convertible_to T2> +static auto operator*(const T1& v1, const T2& v2) +{ + if constexpr (std::is_arithmetic_v) { + return v2 *= v1; + } else if constexpr (std::is_arithmetic_v) { + return v1 *= v2; + } else { + using res_t = std::conditional_t && std::derived_from, T1, T2>; + return res_t{(double)v1 * (double)v2}; + // static_assert(false, "INCOMPATIBLE TYPES!"); + } +} + +template T1, typename T2> +static auto operator/(const T1& v1, const T2& v2) + requires std::is_arithmetic_v +{ + return v1 /= v2; +} + + +/* */ + +class MccAngleRA_ICRS : public MccAngle +{ + using MccAngle::MccAngle; +}; + +class MccAngleDEC_ICRS : public MccAngle +{ + using MccAngle::MccAngle; +}; + +class MccAngleRA_APP : public MccAngle +{ + using MccAngle::MccAngle; +}; + +class MccAngleDEC_APP : public MccAngle +{ + using MccAngle::MccAngle; +}; + +class MccAngleHA : public MccAngle +{ + using MccAngle::MccAngle; +}; + +class MccAngleAZ : public MccAngle +{ + using MccAngle::MccAngle; +}; + +class MccAngleZD : public MccAngle +{ +public: + using MccAngle::MccAngle; +}; + +class MccAngleALT : public MccAngle +{ +public: + using MccAngle::MccAngle; +}; + + +class MccAngleX : public MccAngle // some co-longitude coordinate +{ + using MccAngle::MccAngle; +}; + +class MccAngleY : public MccAngle // some co-latitude coordinate +{ + using MccAngle::MccAngle; +}; + + + +class MccAngleLAT : public MccAngle +{ + using MccAngle::MccAngle; +}; + +class MccAngleLON : public MccAngle +{ + using MccAngle::MccAngle; +}; + + +enum class MccCoordKind : size_t { + COORDS_KIND_GENERIC = traits::mcc_type_hash, + COORDS_KIND_RA_ICRS = traits::mcc_type_hash, + COORDS_KIND_DEC_ICRS = traits::mcc_type_hash, + COORDS_KIND_RA_APP = traits::mcc_type_hash, + COORDS_KIND_DEC_APP = traits::mcc_type_hash, + COORDS_KIND_HA = traits::mcc_type_hash, + COORDS_KIND_AZ = traits::mcc_type_hash, + COORDS_KIND_ZD = traits::mcc_type_hash, + COORDS_KIND_ALT = traits::mcc_type_hash, + COORDS_KIND_X = traits::mcc_type_hash, + COORDS_KIND_Y = traits::mcc_type_hash, + COORDS_KIND_LAT = traits::mcc_type_hash, + COORDS_KIND_LON = traits::mcc_type_hash +}; + +enum class MccCoordPairKind : size_t { + COORDS_KIND_GENERIC = traits::mcc_type_pair_hash(), + COORDS_KIND_RADEC_ICRS = traits::mcc_type_pair_hash(), + COORDS_KIND_RADEC_APP = traits::mcc_type_pair_hash(), + COORDS_KIND_HADEC_APP = traits::mcc_type_pair_hash(), + COORDS_KIND_AZZD = traits::mcc_type_pair_hash(), + COORDS_KIND_AZALT = traits::mcc_type_pair_hash(), + COORDS_KIND_XY = traits::mcc_type_pair_hash(), + COORDS_KIND_LATLON = traits::mcc_type_pair_hash() +}; + + +template +static constexpr std::string_view MccCoordPairKindStr = + KIND == MccCoordPairKind::COORDS_KIND_RADEC_ICRS ? "RADEC-IRCS" + : KIND == MccCoordPairKind::COORDS_KIND_RADEC_APP ? "RADEC-APP" + : KIND == MccCoordPairKind::COORDS_KIND_HADEC_APP ? "HADEC-APP" + : KIND == MccCoordPairKind::COORDS_KIND_AZALT ? "Azimuth-Altitude" + : KIND == MccCoordPairKind::COORDS_KIND_AZZD ? "Azimuth-Zendist" + : KIND == MccCoordPairKind::COORDS_KIND_XY ? "X-Y" + : KIND == MccCoordPairKind::COORDS_KIND_LATLON ? "Latitude-Longitude" + : "UNKNOWN"; + + +} // namespace mcc diff --git a/mcc/mcc_ccte_erfa.h b/mcc/mcc_ccte_erfa.h index 619b0f3..abbcdd8 100644 --- a/mcc/mcc_ccte_erfa.h +++ b/mcc/mcc_ccte_erfa.h @@ -10,6 +10,7 @@ #include +#include "mcc_angle.h" #include "mcc_ccte_iers.h" #include "mcc_defaults.h" @@ -116,6 +117,8 @@ inline std::error_code make_error_code(MccCCTE_ERFAErrorCode ec) class MccCCTE_ERFA { + static constexpr double PI_2 = std::numbers::pi / 2.0; + public: static constexpr double DEFAULT_WAVELENGTH = 0.55; // default observed wavelength in mkm @@ -160,6 +163,14 @@ public: MccCCTE_ERFA() : _stateMutex(new std::mutex) {} + MccCCTE_ERFA(engine_state_t state) : _currentState(std::move(state)), _stateMutex(new std::mutex) {} + + MccCCTE_ERFA(const MccCCTE_ERFA&) = delete; + MccCCTE_ERFA& operator=(const MccCCTE_ERFA&) = delete; + + MccCCTE_ERFA(MccCCTE_ERFA&&) = default; + MccCCTE_ERFA& operator=(MccCCTE_ERFA&&) = default; + virtual ~MccCCTE_ERFA() = default; // engine state related methods @@ -344,60 +355,12 @@ public: from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP || from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { // - ret = timepointToJulday(from_pt.time_point, &jd); - if (ret) { - return ret; - } - - std::lock_guard lock{*_stateMutex}; - - auto dut1 = _currentState._bulletinA.DUT1(jd.mjd); - - if (!dut1.has_value()) { - return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; - } - - auto pol_pos = _currentState._bulletinA.polarCoords(jd.mjd); - if (!pol_pos.has_value()) { - return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; - } - const auto arcsec2rad = std::numbers::pi / 180 / 3600; - pol_pos->x *= arcsec2rad; - pol_pos->y *= arcsec2rad; - - std::string type; - switch (from_pt.pair_kind) { - case mcc::MccCoordPairKind::COORDS_KIND_AZZD: - type = "A"; - break; - case mcc::MccCoordPairKind::COORDS_KIND_AZALT: - from_pt.Y = std::numbers::pi / 2.0 - from_pt.Y; // altitude to zenithal distance - type = "A"; - break; - case mcc::MccCoordPairKind::COORDS_KIND_HADEC_APP: - type = "H"; - break; - case mcc::MccCoordPairKind::COORDS_KIND_RADEC_APP: - type = "R"; - break; - default: - return MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; - } - - int err = eraAtoc13(type.c_str(), from_pt.X, from_pt.Y, jd.MJD0, jd.mjd, dut1->count(), - _currentState.lon, _currentState.lat, _currentState.elev, pol_pos->x, pol_pos->y, - _currentState.meteo.pressure, _currentState.meteo.temperature, - _currentState.meteo.humidity, _currentState.wavelength, &to_pt->X, &to_pt->Y); - - if (err == 1) { - return MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR; - } else if (err == -1) { - return MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE; - } } else { ret = MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; } + ret = obs2icrs(from_pt, to_pt); + return ret; } @@ -409,53 +372,24 @@ public: to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP || to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { // - ret = timepointToJulday(to_pt->time_point, &jd); - if (ret) { - return ret; - } + MccEqtHrzCoords eqhrz; - std::lock_guard lock{*_stateMutex}; + ret = icrs2obs(from_pt, &eqhrz); - auto dut1 = _currentState._bulletinA.DUT1(jd.mjd); - - if (!dut1.has_value()) { - return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; - } - - auto pol_pos = _currentState._bulletinA.polarCoords(jd.mjd); - if (!pol_pos.has_value()) { - return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; - } - const auto arcsec2rad = std::numbers::pi / 180 / 3600; - pol_pos->x *= arcsec2rad; - pol_pos->y *= arcsec2rad; - - double oaz, ozd, oha, odec, ora, eo_; - - int err = eraAtco13(from_pt.X, from_pt.Y, 0.0, 0.0, 0.0, 0.0, jd.MJD0, jd.mjd, dut1->count(), - _currentState.lon, _currentState.lat, _currentState.elev, pol_pos->x, pol_pos->y, - _currentState.meteo.pressure, _currentState.meteo.temperature, - _currentState.meteo.humidity, _currentState.wavelength, &oaz, &ozd, &oha, &odec, - &ora, &eo_); - - if (err == 1) { - return MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR; - } else if (err == -1) { - return MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE; - } - - if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZALT) { - to_pt->X = oaz; - to_pt->Y = std::numbers::pi / 2.0 - ozd; - } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZZD) { - to_pt->X = oaz; - to_pt->Y = ozd; - } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { - to_pt->X = ora; - to_pt->Y = odec; - } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { - to_pt->X = oha; - to_pt->Y = odec; + if (!ret) { + if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZALT) { + to_pt->X = eqhrz.AZ; + to_pt->Y = eqhrz.ALT; + } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZZD) { + to_pt->X = eqhrz.AZ; + to_pt->Y = eqhrz.ZD; + } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { + to_pt->X = eqhrz.RA_APP; + to_pt->Y = eqhrz.DEC_APP; + } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { + to_pt->X = eqhrz.HA; + to_pt->Y = eqhrz.DEC_APP; + } } } else { ret = MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; @@ -463,6 +397,202 @@ public: return ret; } + + + // from apparent to apparent + if (from_pt.time_point != to_pt->time_point) { // first, convert 'from' coordinates to ICRS + MccCelestialPoint pt{.time_point = to_pt->time_point, + .pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_ICRS}; + ret = obs2icrs(from_pt, &pt); + + if (!ret) { // from ICRS to required + return transformCoordinates(pt, to_pt); + } + } else { // the same time points + double eo, lst, ha, dec, az, alt; + + auto lst_eo = [&]() { + ret = eqOrigins(from_pt.time_point, &eo); + if (!ret) { + ret = timepointToAppSideral(from_pt.time_point, &lst, true); + } + }; + + if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { + // first, compute HA from CIO-based RA!!! + lst_eo(); + if (!ret) { + ha = lst - from_pt.X + eo; + } else { + return ret; + } + + if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZZD) { + eraHd2ae(ha, from_pt.Y, _currentState.lat, &az, &alt); + to_pt->X = az; + to_pt->Y = PI_2 - alt; + } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZALT) { + eraHd2ae(ha, from_pt.Y, _currentState.lat, &az, &alt); + to_pt->X = az; + to_pt->Y = alt; + } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { + to_pt->X = ha; + to_pt->Y = from_pt.Y; + } else { + ret = MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; + } + } else if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { + if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZZD) { + eraHd2ae(ha, from_pt.Y, _currentState.lat, &az, &alt); + to_pt->X = az; + to_pt->Y = PI_2 - alt; + } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZALT) { + eraHd2ae(ha, from_pt.Y, _currentState.lat, &az, &alt); + to_pt->X = az; + to_pt->Y = alt; + } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { + lst_eo(); + if (!ret) { + to_pt->X = lst - from_pt.X + eo; + to_pt->Y = from_pt.Y; + } else { + return ret; + } + } else { + ret = MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; + } + + } else if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_AZZD) { + from_pt.Y = PI_2 - from_pt.Y; + from_pt.pair_kind = MccCoordPairKind::COORDS_KIND_AZALT; + + ret = transformCoordinates(std::move(from_pt), to_pt); + } else if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_AZALT) { + if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { + eraAe2hd(from_pt.X, from_pt.Y, _currentState.lat, &ha, &dec); + to_pt->X = ha; + to_pt->Y = dec; + } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZZD) { + to_pt->X = from_pt.X; + to_pt->Y = PI_2 - from_pt.Y; + } else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { + eraAe2hd(from_pt.X, from_pt.Y, _currentState.lat, &ha, &dec); + lst_eo(); + if (!ret) { + to_pt->X = lst - ha + eo; + to_pt->Y = dec; + } + } else { + ret = MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; + } + + } else { + ret = MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; + } + } + + return ret; + } + + + error_t transformCoordinates(mcc_celestial_point_c auto from_pt, mcc_eqt_hrz_coord_c auto* to_pt) + { + error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK; + + MccJulianDay jd; + + if (to_pt == nullptr) { + return ret; + } + + using tp_pt_t = std::remove_cvref_t; + MccGenericCelestialPoint cpt; + cpt.time_point = from_pt.time_point; + + // the main scenario: from ICRS to apparent + if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { + return icrs2obs(from_pt, to_pt); + } + + // from apparent: copy corresponded coordinates and compute other ones + if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { + to_pt->RA_APP = from_pt.X; + to_pt->DEC_APP = from_pt.Y; + + cpt.pair_kind = MccCoordPairKind::COORDS_KIND_HADEC_APP; + + ret = transformCoordinates(from_pt, &cpt); + if (!ret) { + to_pt->HA = cpt.X; + + auto cpt1 = cpt; + cpt1.pair_kind = MccCoordPairKind::COORDS_KIND_AZALT; + + ret = transformCoordinates(cpt, &cpt1); + if (!ret) { + to_pt->AZ = cpt1.X; + to_pt->ALT = cpt1.Y; + to_pt->ZD = PI_2 - cpt1.Y; + } + } + } else if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { + to_pt->HA = from_pt.X; + to_pt->DEC_APP = from_pt.Y; + + cpt.pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_APP; + ret = transformCoordinates(from_pt, &cpt); + if (!ret) { + to_pt->RA_APP = cpt.X; + + cpt.pair_kind = MccCoordPairKind::COORDS_KIND_AZALT; + ret = transformCoordinates(from_pt, &cpt); + if (!ret) { + to_pt->AZ = cpt.X; + to_pt->ALT = cpt.Y; + to_pt->ZD = PI_2 - cpt.Y; + } + } + } else if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_AZZD) { + to_pt->AZ = from_pt.X; + to_pt->ZD = from_pt.Y; + to_pt->ALT = PI_2 - from_pt.Y; + + cpt.pair_kind = MccCoordPairKind::COORDS_KIND_HADEC_APP; + ret = transformCoordinates(from_pt, &cpt); + if (!ret) { + to_pt->HA = cpt.X; + to_pt->DEC_APP = cpt.Y; + + auto cpt1 = cpt; + cpt1.pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_APP; + ret = transformCoordinates(cpt, &cpt1); + if (!ret) { + to_pt->RA_APP = cpt1.X; + } + } + } else if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_AZALT) { + to_pt->AZ = from_pt.X; + to_pt->ALT = from_pt.Y; + to_pt->ZD = PI_2 - from_pt.Y; + + cpt.pair_kind = MccCoordPairKind::COORDS_KIND_HADEC_APP; + ret = transformCoordinates(from_pt, &cpt); + if (!ret) { + to_pt->HA = cpt.X; + to_pt->DEC_APP = cpt.Y; + + auto cpt1 = cpt; + cpt1.pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_APP; + ret = transformCoordinates(cpt, &cpt1); + if (!ret) { + to_pt->RA_APP = cpt1.X; + } + } + } else { + ret = MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; + } + + return ret; } @@ -536,6 +666,154 @@ protected: engine_state_t _currentState{}; std::unique_ptr _stateMutex; + + error_t icrs2obs(mcc_celestial_point_c auto const& pt, mcc_eqt_hrz_coord_c auto* result) + { + if (result == nullptr) { + return MccCCTE_ERFAErrorCode::ERROR_OK; + } + + MccJulianDay jd; + + auto ret = timepointToJulday(result->time_point, &jd); + if (ret) { + return ret; + } + + std::lock_guard lock{*_stateMutex}; + + auto dut1 = _currentState._bulletinA.DUT1(jd.mjd); + + if (!dut1.has_value()) { + return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; + } + + auto pol_pos = _currentState._bulletinA.polarCoords(jd.mjd); + if (!pol_pos.has_value()) { + return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; + } + + // const auto arcsec2rad = std::numbers::pi / 180 / 3600; + const auto arcsec2rad = 1.0_arcsecs; + pol_pos->x *= arcsec2rad; + pol_pos->y *= arcsec2rad; + + double az, zd, ha, ra, dec, eo; + + int err = eraAtco13(pt.X, pt.Y, 0.0, 0.0, 0.0, 0.0, jd.MJD0, jd.mjd, dut1->count(), _currentState.lon, + _currentState.lat, _currentState.elev, pol_pos->x, pol_pos->y, _currentState.meteo.pressure, + _currentState.meteo.temperature, _currentState.meteo.humidity, _currentState.wavelength, + &az, &zd, &ha, &dec, &ra, &eo); + + if (err == 1) { + ret = MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR; + } else if (err == -1) { + ret = MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE; + } + + result->AZ = az; + result->ZD = zd; + result->HA = ha; + result->RA_APP = ra; + result->DEC_APP = dec; + result->ALT = std::numbers::pi / 2.0 - zd; + + return ret; + } + + + error_t obs2icrs(mcc_celestial_point_c auto const& from_pt, mcc_celestial_point_c auto* to_pt) + { + if (to_pt == nullptr) { + return MccCCTE_ERFAErrorCode::ERROR_OK; + } + + MccJulianDay jd; + + auto ret = timepointToJulday(from_pt.time_point, &jd); + if (ret) { + return ret; + } + + std::lock_guard lock{*_stateMutex}; + + auto dut1 = _currentState._bulletinA.DUT1(jd.mjd); + + if (!dut1.has_value()) { + return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; + } + + auto pol_pos = _currentState._bulletinA.polarCoords(jd.mjd); + if (!pol_pos.has_value()) { + return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; + } + + // const auto arcsec2rad = std::numbers::pi / 180 / 3600; + const auto arcsec2rad = 1.0_arcsecs; + pol_pos->x *= arcsec2rad; + pol_pos->y *= arcsec2rad; + + std::string type; + switch (from_pt.pair_kind) { + case mcc::MccCoordPairKind::COORDS_KIND_AZZD: + type = "A"; + break; + case mcc::MccCoordPairKind::COORDS_KIND_AZALT: + from_pt.Y = std::numbers::pi / 2.0 - from_pt.Y; // altitude to zenithal distance + type = "A"; + break; + case mcc::MccCoordPairKind::COORDS_KIND_HADEC_APP: + type = "H"; + break; + case mcc::MccCoordPairKind::COORDS_KIND_RADEC_APP: + type = "R"; + break; + default: + return MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; + } + + int err = eraAtoc13(type.c_str(), from_pt.X, from_pt.Y, jd.MJD0, jd.mjd, dut1->count(), _currentState.lon, + _currentState.lat, _currentState.elev, pol_pos->x, pol_pos->y, _currentState.meteo.pressure, + _currentState.meteo.temperature, _currentState.meteo.humidity, _currentState.wavelength, + &to_pt->X, &to_pt->Y); + + if (err == 1) { + return MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR; + } else if (err == -1) { + return MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE; + } + } + + error_t eqOrigins(mcc_time_point_c auto const& tp, double* eo) + { + if (eo == nullptr) { + return MccCCTE_ERFAErrorCode::ERROR_OK; + } + + MccJulianDay jd, tt; + auto ret = timepointToJulday(tp, &jd); + if (ret) { + return ret; + } + + std::lock_guard lock{*_stateMutex}; + + using real_days_t = std::chrono::duration>; + + auto tai_utc = _currentState._leapSeconds[jd.mjd]; + if (tai_utc.has_value()) { + tt.mjd = jd.mjd + std::chrono::duration_cast(tai_utc.value()).count(); + + auto tt_tai = _currentState._bulletinA.TT_TAI(); + tt.mjd = jd.mjd + std::chrono::duration_cast(tt_tai).count(); + + *eo = eraEo06a(tt.MJD0, tt.mjd); + } else { + ret = MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE; + } + + return ret; + } }; } // namespace mcc::ccte::erfa diff --git a/mcc/mcc_generics.h b/mcc/mcc_generics.h index a473c5b..3c8fc39 100644 --- a/mcc/mcc_generics.h +++ b/mcc/mcc_generics.h @@ -11,7 +11,8 @@ #include -#include "mcc_traits.h" +// #include "mcc_traits.h" +#include "mcc_angle.h" namespace mcc { @@ -20,16 +21,16 @@ namespace mcc // mount construction type (only the most common ones) enum class MccMountType : uint8_t { GERMAN_TYPE, FORK_TYPE, CROSSAXIS_TYPE, ALTAZ_TYPE }; -enum MccCoordPairKind : size_t { - COORDS_KIND_GENERIC, - COORDS_KIND_RADEC_ICRS, - COORDS_KIND_RADEC_APP, - COORDS_KIND_HADEC_APP, - COORDS_KIND_AZZD, - COORDS_KIND_AZALT, - COORDS_KIND_XY, - COORDS_KIND_LATLON -}; +// enum MccCoordPairKind : size_t { +// COORDS_KIND_GENERIC, +// COORDS_KIND_RADEC_ICRS, +// COORDS_KIND_RADEC_APP, +// COORDS_KIND_HADEC_APP, +// COORDS_KIND_AZZD, +// COORDS_KIND_AZALT, +// COORDS_KIND_XY, +// COORDS_KIND_LATLON +// }; /* FLOATING-POINT LIKE CLASS CONCEPT */ @@ -165,14 +166,13 @@ protected: template -concept mcc_ast_engine_c = - std::derived_from> && requires(const T t_const, T t) { - { t_const.name() } -> std::formattable; +concept mcc_ccte_c = std::derived_from> && requires(const T t_const, T t) { + { t_const.name() } -> std::formattable; - requires mcc_refract_model_c; + requires mcc_refract_model_c; - { t.refractionModel(std::declval()) } -> std::same_as; - }; + { t.refractionModel(std::declval()) } -> std::same_as; +}; @@ -303,7 +303,8 @@ struct mcc_pzone_interface_t { template SelfT, typename InputT> RetT inPZone(this SelfT&& self, InputT coords, bool* result) - requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) + requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) && + requires { self.inPZone(coords, result); } { return std::forward(self).InPZone(std::move(coords), result); } @@ -311,14 +312,16 @@ struct mcc_pzone_interface_t { template SelfT, typename InputT> RetT timeToPZone(this SelfT&& self, InputT coords, traits::mcc_time_duration_c auto* res_time) - requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) + requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) && + requires { self.timeToPZone(coords, res_time); } { return std::forward(self).timeToPZone(std::move(coords), res_time); } template SelfT, typename InputT> RetT timeFromPZone(this SelfT&& self, InputT coords, traits::mcc_time_duration_c auto* res_time) - requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) + requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) && + requires { self.timeFromPZone(coords, res_time); } { return std::forward(self).timeFromPZone(std::move(coords), res_time); } diff --git a/mcc/mcc_pzone.h b/mcc/mcc_pzone.h new file mode 100644 index 0000000..dd9af90 --- /dev/null +++ b/mcc/mcc_pzone.h @@ -0,0 +1,78 @@ +#pragma once + +/* MOUNT CONTROL COMPONENTS LIBRARY */ + + +/* IMPLEMENTATION OF SOME SIMPLE PROHIBITED ZONES */ + +#include "mcc_generics.h" + +namespace mcc +{ + +static constexpr double mcc_sideral_to_UT1_ratio = 1.002737909350795; // sideral/UT1 + + +/* MINIMAL OR MAXIMAL ALTITUDE PROHIBITED ZONES */ + +enum class MccAltLimitKind { MIN_ALT_LIMIT, MAX_ALT_LIMIT }; + +template +class MccAltLimitPZ : public mcc_pzone_interface_t +{ +protected: + double _altLimit, _cosALim, _sinAlim; + double _cosLat, _sinLat, _absLat; + + CCTE_T* _ccteEngine; + +public: + typedef std::error_code error_t; + + MccAltLimitPZ(mcc_angle_c auto const& alt_limit, mcc_angle_c auto const& latitude, CCTE_T* ccte_engine) + : _altLimit(alt_limit), + _cosALim(cos(_altLimit)), + _sinAlim(sin(_altLimit)), + _cosLat(cos(latitude)), + _sinLat(sin(latitude)), + _absLat(abs(latitude)), + _ccteEngine(ccte_engine) + { + } + + consteval std::string_view name() const + { + return KIND == MccAltLimitKind::MIN_ALT_LIMIT ? "MINALT-ZONE" + : KIND == MccAltLimitKind::MAX_ALT_LIMIT ? "MAXALT-ZONE" + : "ALTLIMIT-UNKNOWN"; + } + + + template + error_t inPZone(InputT coords, bool* result) + requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) + { + // auto ret = + // return ret; + } + + template + error_t timeToPZone(InputT coords, traits::mcc_time_duration_c auto* res_time) + requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) + { + } + + template + error_t timeFromPZone(InputT coords, traits::mcc_time_duration_c auto* res_time) + requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) + { + } + + template + error_t intersectPZone(InputT coords, mcc_celestial_point_c auto* point) + requires(mcc_eqt_hrz_coord_c || mcc_celestial_point_c) + { + } +}; + +} // namespace mcc