#include "mcc_angle.h" #include "mcc_ccte_iers.h" #include "mcc_defaults.h" #include "mcc_generics.h" #include namespace mcc { template class MccCoordPair { public: typedef CO_LON_T x_t; typedef CO_LAT_T y_t; static constexpr MccCoordPairKind pairKind = !(std::derived_from || std::derived_from) // unknown type (possibly just double or float) ? MccCoordPairKind::COORDS_KIND_GENERIC : (std::same_as || std::same_as) // one of the types is MccAngle ? MccCoordPairKind::COORDS_KIND_GENERIC // ICRS RA and DEC : (std::same_as && std::same_as) ? MccCoordPairKind::COORDS_KIND_RADEC_ICRS // apparent RA and DEC : (std::same_as && std::same_as) ? MccCoordPairKind::COORDS_KIND_RADEC_APP // observed RA and DEC : (std::same_as && std::same_as) ? MccCoordPairKind::COORDS_KIND_RADEC_APP // 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 // apparent AZ and ZD : (std::same_as && std::same_as) ? MccCoordPairKind::COORDS_KIND_AZZD // apparent AZ and ALT : (std::same_as && std::same_as) ? MccCoordPairKind::COORDS_KIND_AZZD // general purpose X and Y : (std::same_as && std::same_as) ? MccCoordPairKind::COORDS_KIND_XY // geographical longitude and latitude : (std::same_as && std::same_as) ? MccCoordPairKind::COORDS_KIND_LONLAT : MccCoordPairKind::COORDS_KIND_UNKNOWN; template MccCoordPair(CO_LON_T const& x, CO_LAT_T const& y, EpT const& epoch = EpT::now()) : _x(x), _y(y), _mjd(epoch.MJD) { } MccCoordPair(const MccCoordPair&) = default; MccCoordPair(MccCoordPair&&) = default; MccCoordPair& operator=(const MccCoordPair&) = default; MccCoordPair& operator=(MccCoordPair&&) = default; virtual ~MccCoordPair() = default; CO_LON_T x() const { return _x; } CO_LAT_T y() const { return _y; } double MJD() const { return _mjd; } // for something like: // auto [ra, dec, mjd] = coord_pair; operator std::tuple() const { return {_x, _y, _mjd}; } void setX(const CO_LON_T& x) { _x = x; } void setY(const CO_LAT_T& y) { _y = y; } void setMJD(double mjd) { _mjd = mjd; } protected: CO_LON_T _x; CO_LAT_T _y; double _mjd; }; template concept mcc_coord_pair_c = requires { typename T::x_t; typename T::y_t; requires std::derived_from>; }; /* predefined coordinate pairs */ template requires(std::derived_from && std::derived_from) class MccNamedCoordPair : public MccCoordPair { public: template requires(std::is_arithmetic_v && std::is_arithmetic_v) MccNamedCoordPair(CxT const& x, CyT const& y, EpT const& epoch = EpT::now()) : MccCoordPair(CO_LON_T{x}, CO_LAT_T{y}, epoch) { } template MccNamedCoordPair(MccAngle const& x, MccAngle const& y, EpT const& epoch = EpT::now()) : MccCoordPair(CO_LON_T{(double)x}, CO_LAT_T{(double)y}, epoch) { } MccNamedCoordPair(const MccNamedCoordPair&) = default; MccNamedCoordPair(MccNamedCoordPair&&) = default; MccNamedCoordPair& operator=(const MccNamedCoordPair&) = default; MccNamedCoordPair& operator=(MccNamedCoordPair&&) = default; virtual ~MccNamedCoordPair() = default; }; using MccSkyRADEC_ICRS = MccNamedCoordPair; using MccSkyRADEC_APP = MccNamedCoordPair; using MccSkyRADEC_OBS = MccNamedCoordPair; using MccSkyHADEC_APP = MccNamedCoordPair; using MccSkyHADEC_OBS = MccNamedCoordPair; using MccSkyAZZD = MccNamedCoordPair; using MccSkyAZALT = MccNamedCoordPair; using MccGenXY = MccNamedCoordPair; using MccGeoLONLAT = MccNamedCoordPair; struct mcc_skypoint_interface_t { virtual ~mcc_skypoint_interface_t() = default; // template SelfT, mcc_angle_c XT, mcc_angle_c YT> // SelfT& from(this SelfT&& self, XT&& x, YT&& y) // { // return std::forward(self).from(std::forward(x), std::forward(y)); // } template SelfT, mcc_coord_pair_c PT> auto from(this SelfT&& self, PT&& cpair) { return std::forward(self).from(std::forward(cpair)); } template SelfT, mcc_coord_pair_c PT> auto operator=(this SelfT&& self, PT&& cpair) { return std::forward(self).operator=(std::forward(cpair)); } template SelfT, mcc_coord_pair_c PT, mcc_coord_pair_c... PTs> auto to(this SelfT&& self, PT& cpair, PTs&... cpairs) { return std::forward(self).to(cpair, cpairs...); } template SelfT, mcc_coord_pair_c PT> operator PT(this SelfT&& self) { return std::forward(self).operator PT(); } template SelfT, mcc_coord_pair_c PT, mcc_coord_pair_c... PTs> operator std::tuple(this SelfT&& self) { return std::forward(self).operator std::tuple(); } }; template concept mcc_skypoint_c = std::derived_from && requires(T t) { typename T::meteo_t; { T::meteo(std::declval()) }; }; /* MCC-LIBRARY DEFAULT SKY POINT CLASS IMPLEMENTATION BASED ON ERFA-LIBRARY */ class MccSkyPoint : public mcc_skypoint_interface_t { public: static constexpr double MJD0 = 2400000.5; struct meteo_t { double temperature; // Temperature in C double humidity; // humidity in % ([0.0, 1.0]) double pressure; // atmospheric presure in hPa=mB }; static ccte::iers::MccLeapSeconds iersLeapSeconds() { return _leapSeconds; } static bool updateLeapSeconds(traits::mcc_input_char_range auto const& filename) { std::lock_guard lock{_leapSecondsMutex}; return _leapSeconds.load(filename); }; static bool updateIersBulletinA(traits::mcc_input_char_range auto const& filename) { std::lock_guard lock{_bulletinAMutex}; return _bulletinA.load(filename); }; static ccte::iers::MccIersBulletinA iersBulletinA() { return _bulletinA; } static void setMeteo(meteo_t meteo) { std::lock_guard lock{_meteoMutex}; _currentMeteo = std::move(meteo); } static meteo_t getMeteo() { return _currentMeteo; } MccSkyPoint() {} template MccSkyPoint(const PT& coord_pair) : MccSkyPoint() { auto self = from(coord_pair); } MccSkyPoint(const MccSkyPoint&) = default; MccSkyPoint(MccSkyPoint&&) = default; MccSkyPoint& operator=(const MccSkyPoint&) = default; MccSkyPoint& operator=(MccSkyPoint&&) = default; virtual ~MccSkyPoint() = default; MccCelestialCoordEpoch epoch() const { return _epoch; } template MccSkyPoint& from(const PT& coord_pair) { _x = coord_pair.x(); _y = coord_pair.y(); _pairKind = coord_pair.pairKind; if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { _epoch = MccCelestialCoordEpoch(); // J2000.0 } else { _epoch.fromMJD(coord_pair.MJD()); } } MccSkyPoint& operator=(mcc_coord_pair_c auto const& coord_pair) { return from(coord_pair); } template auto to(PT& cpair, PTs&... cpairs) { toHelper(cpair); if constexpr (sizeof...(PTs)) { to(cpairs...); } } protected: // IERS related static members static inline ccte::iers::MccLeapSeconds _leapSeconds{}; static inline ccte::iers::MccIersBulletinA _bulletinA{}; static inline std::mutex _leapSecondsMutex{}, _bulletinAMutex{}; // meteo related static members static inline meteo_t _currentMeteo{.temperature = 10.0, .humidity = 0.5, .pressure = 1010.0}; static inline std::mutex _meteoMutex{}; double _x{0.0}, _y{0.0}; MccCoordPairKind _pairKind{MccCoordPairKind::COORDS_KIND_RADEC_ICRS}; MccCelestialCoordEpoch _epoch{}; // J2000.0 template auto toHelper(PT& cpair) { if (this == &cpair) { return; } static constexpr double half_pi = std::numbers::pi / 2.0; // HA, DEC to AZ, ALT (AZ from the South through the West) auto hadec2azalt = [](double ha, double dec, double phi, double& az, double& alt) { const auto cos_phi = std::cos(phi), sin_phi = std::sin(phi); const auto cos_dec = std::cos(dec), sin_dec = std::sin(dec); const auto cos_ha = std::cos(ha), sin_ha = std::sin(ha); auto x = sin_phi * cos_dec * cos_ha - cos_phi * sin_dec; auto y = -cos_dec * sin_ha; auto z = cos_phi * cos_dec * cos_ha + sin_phi * sin_dec; auto xx = x * x, yy = y * y; decltype(x) r; if (xx < yy) { r = yy * sqrt(1.0 + xx / yy); } else { r = xx * sqrt(1.0 + yy / xx); } az = utils::isEqual(r, 0.0) ? 0.0 : std::atan2(y, x); if (az < 0.0) { az += std::numbers::pi * 2.0; // to range of [0, 2*PI] } alt = std::atan2(z, r); }; // AZ, ALT to HA, DEC (AZ from the South through the West) auto azalt2hadec = [](double az, double alt, double phi, double& ha, double& dec) { const auto cos_phi = std::cos(phi), sin_phi = std::sin(phi); const auto cos_az = std::cos(az), sin_az = std::sin(az); const auto cos_alt = std::cos(alt), sin_alt = std::sin(alt); auto x = sin_phi * cos_alt * cos_az + cos_phi * sin_alt; auto y = cos_alt * sin_az; auto z = -cos_phi * cos_alt * cos_az + sin_phi * sin_alt; auto xx = x * x, yy = y * y; decltype(x) r; if (xx < yy) { r = yy * sqrt(1.0 + xx / yy); } else { r = xx * sqrt(1.0 + yy / xx); } ha = utils::isEqual(r, 0.0) ? 0.0 : std::atan2(y, x); dec = std::atan2(z, r); }; double phi; double ra_icrs, dec_icrs, ra, dec, ha, az, zd, alt, lst, eo; static_assert(PT::pairKind == MccCoordPairKind::COORDS_KIND_GENERIC, "UNSUPPORTED SKY POINT TRANSFORMATION!"); static_assert(PT::pairKind == MccCoordPairKind::COORDS_KIND_UNKNOWN, "UNSUPPORTED SKY POINT TRANSFORMATION!"); // from ICRS to ICRS - just copy and exit if (_pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS && PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { cpair = PT(PT::x_t(_x), PT::y_t(_y), _epoch); return; } // just copy coordinates and exit if (_pairKind == PT::pairKind && utils::isEqual(_epoch.MJD(), cpair.MJD())) { cpair = PT(PT::x_t(_x), PT::y_t(_y), _epoch); return; } // if epochs is not the same then // 1) convert stored coordinates to ICRS ones // 2) convert from the computed ICRS coordinates to required ones MccCoordPairKind pkind = _pairKind; if (!utils::isEqual(_epoch.MJD(), cpair.MJD())) { // convert stored pair to ICRS one (ra_icrs, dec_icrs) if (_pairKind != MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { pkind = MccCoordPairKind::COORDS_KIND_RADEC_ICRS; // cct_engine.appToICRS(app_type, app_x, app_y, ra_icrs, dec_icrs) // cct_engine.obsToICRS(obs_type, obs_x, obs_y, ra_icrs, dec_icrs) if (mcc_is_obs_coordpair(_pairKind)) { // cct_engine.obsToICRS(...) } else if (mcc_is_app_coordpair(_pairKind)) { // cct_engine.appToICRS(...) } else { // unsupported transformation!!! return; } } else { ra_icrs = _x; dec_icrs = _y; } } // here, from APP or OBS to ICRS and exit if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS && PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { cpair = PT(PT::x_t(ra_icrs), PT::y_t(dec_icrs), MccCelestialCoordEpoch{}); return; } // here, the input coordinates and stored one are at the same epoch if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { // cct_engine.icrsToObs(ra_icrs, dec_icrs, jd, ra_obs, dec_obs, ha, az, zd) // cct_engine.icrsToApp(ra_icrs, dec_icrs, jd, ra_app, dec_app, ha) if constexpr (mccIsObsCoordPairKind) { // cct_engine.icrsToObs(...) } else if constexpr (mccIsAppCoordPairKind) { // cct_engine.icrsToApp(...) } else { static_assert(false, "UNSUPPORTED SKY POINT TRANSFORMATION!"); } if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { // cct_engine.icrsToApp(...) cpair.setX(ra); cpair.setY(dec); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { // cct_engine.icrsToObs(...) } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { // cct_engine.icrsToApp(...) cpair.setX(ha); cpair.setY(dec); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { // cct_engine.icrsToObs(...) } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) { // cct_engine.icrsToObs(...) cpair.setX(az); cpair.setY(zd); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) { // cct_engine.icrsToObs(...) cpair.setX(az); cpair.setY(half_pi - zd); } else { static_assert(false, "UNSUPPORTED SKY POINT TRANSFORMATION!"); } } else if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { // compute EO, LST, refract model ha = lst + eo - cpair.x(); // from CIO based RA if constexpr (PT::pairKind != MccCoordPairKind::COORDS_KIND_HADEC_APP) { hadec2azalt(ha, _y, phi, az, alt); // to app az, alt // to observed zenithal distance: alt += z_corr } if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { azalt2hadec(az, alt, phi, ha, dec); // to obs ha, dec ra = lst + eo - ha; // CIO based RA cpair.setX(ra); cpair.setY(dec); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { cpair.setX(ha); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { azalt2hadec(az, alt, phi, ha, dec); // to obs ha, dec cpair.setX(ha); cpair.setY(dec); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) { cpair.setX(az); // ????????????!!!!!!!!! cpair.setY(half_pi - alt); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) { cpair.setX(az); // ????????????!!!!!!!!! cpair.setY(alt); } else { static_assert(false, "UNSUPPORTED SKY POINT TRANSFORMATION!"); } } else if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { ha = lst + eo - cpair.x(); // from CIO based RA if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { hadec2azalt(ha, _y, phi, az, alt); // to apparent zenithal distance: alt -= z_corr; azalt2hadec(az, alt, phi, ha, dec); // to app ha,dec ra = lst + eo - ha; // CIO based RA cpair.setX(ra); cpair.setY(dec); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { hadec2azalt(ha, _y, phi, az, alt); // to apparent zenithal distance: alt -= z_corr; azalt2hadec(az, alt, phi, ha, dec); // to app ha,dec cpair.setX(ha); cpair.setY(dec); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { cpair.setX(ha); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) { hadec2azalt(ha, _y, phi, az, alt); cpair.setX(az); cpair.setY(half_pi - alt); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) { hadec2azalt(ha, _y, phi, az, alt); cpair.setX(az); cpair.setY(alt); } else { static_assert(false, "UNSUPPORTED SKY POINT TRANSFORMATION!"); } } else if (pkind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { } else if (pkind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { } else if (pkind == MccCoordPairKind::COORDS_KIND_AZZD) { } else if (pkind == MccCoordPairKind::COORDS_KIND_AZALT) { } else { // unsupported transformation!!! return; } if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_APP || pkind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { ra = _x; dec = _y; } else if (pkind == MccCoordPairKind::COORDS_KIND_HADEC_APP || pkind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { ha = _x; dec = _y; } else if (pkind == MccCoordPairKind::COORDS_KIND_AZZD) { az = _x; zd = _y; } else if (pkind == MccCoordPairKind::COORDS_KIND_AZALT) { az = _x; alt = _y; } else { // unsupported transformation!!! return; } auto comp_func = [&, this](this auto& self, MccCoordPairKind cp_kind) { if (cp_kind == MccCoordPairKind::COORDS_KIND_AZALT) { if constexpr (mccIsAppCoordPairKind) { // correct for refraction: alt -= dz_refr } azalt2hadec(az, alt, phi, ha, dec); cpair.setY(dec); if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { ra = lst + eo - ha; cpair.setX(ra); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { ra = lst + eo - ha; cpair.setX(ra); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { cpair.setX(ha); } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { cpair.setX(ha); // } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) { } else { // unsupported transformation!!! return; } } else if (cp_kind == MccCoordPairKind::COORDS_KIND_AZZD) { alt = half_pi - _y; self(MccCoordPairKind::COORDS_KIND_AZALT); } else if (cp_kind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { ra = lst + eo - ha; cpair.setX(ra); cpair.setY(dec); } else { hadec2azalt(ha, dec, phi, az, alt); if constexpr (mccIsAppCoordPairKind) { // RADEC_APP, HADEC_APP self(MccCoordPairKind::COORDS_KIND_AZALT); } else { // AZALT, AZZD cpair.setX(az); if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) { zd = half_pi - alt; cpair.setY(zd); } else { cpair.setY(alt); } } } } else if (cp_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { ra = lst + eo - ha; cpair.setX(ra); cpair.setY(dec); } else { hadec2azalt(ha, dec, phi, az, alt); if constexpr (mccIsObsCoordPairKind) { // RADEC_OBS, HADEC_OBS, AZALT, AZZD // correct for refraction: alt += dz_refr self(MccCoordPairKind::COORDS_KIND_AZALT); } } } else if (cp_kind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) { ha = lst + eo - ra; if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) { cpair.setX(ha); cpair.setY(dec); } else { self(MccCoordPairKind::COORDS_KIND_HADEC_OBS); } } else if (cp_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { ha = lst + eo - ra; if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { cpair.setX(ha); cpair.setY(dec); } else { self(MccCoordPairKind::COORDS_KIND_HADEC_APP); } } }; } }; } // end namespace mcc