#pragma once #include "mcc_angle.h" #include "mcc_ccte_erfa_new.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_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_OBS // 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_AZALT // 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), _epoch(epoch) { } 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; } MccCelestialCoordEpoch epoch() const { return _epoch; } double MJD() const { return _epoch.MJD(); } // for something like: // auto [ra, dec, epoch] = coord_pair; operator std::tuple() const { return {_x, _y, _epoch}; } void setX(const CO_LON_T& x) { _x = x; } void setY(const CO_LAT_T& y) { _y = y; } void setEpoch(mcc_coord_epoch_c auto const& ep) { _epoch = ep; } protected: CO_LON_T _x; CO_LAT_T _y; MccCelestialCoordEpoch _epoch; }; template concept mcc_coord_pair_c = requires { requires mcc_angle_c; requires mcc_angle_c; requires std::derived_from>; }; /* predefined coordinate pairs */ template requires(std::derived_from && std::derived_from) 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()) : MccCoordPair(CO_LON_T{(double)x}, CO_LAT_T{(double)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; }; struct MccSkyRADEC_ICRS : MccNamedCoordPair { template requires(std::is_arithmetic_v && std::is_arithmetic_v) MccSkyRADEC_ICRS(CxT const& x, CyT const& y) : MccNamedCoordPair(x, y, MccCelestialCoordEpoch{}) { } 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&) { static_assert(false, "CANNOT SET EPOCH FOR ICRS-KIND COORDINATE PAIR!!!"); } }; struct MccSkyRADEC_APP : MccNamedCoordPair { using MccNamedCoordPair::MccNamedCoordPair; }; typedef MccNamedCoordPair MccSkyRADEC_OBS; // struct MccSkyRADEC_OBS : MccNamedCoordPair { // using MccNamedCoordPair::MccNamedCoordPair; // }; struct MccSkyHADEC_APP : MccNamedCoordPair { using MccNamedCoordPair::MccNamedCoordPair; }; struct MccSkyHADEC_OBS : MccNamedCoordPair { using MccNamedCoordPair::MccNamedCoordPair; }; struct MccSkyAZZD : MccNamedCoordPair { using MccNamedCoordPair::MccNamedCoordPair; }; struct MccSkyAZALT : MccNamedCoordPair { using MccNamedCoordPair::MccNamedCoordPair; }; struct MccGenXY : MccNamedCoordPair { using MccNamedCoordPair::MccNamedCoordPair; }; struct MccGeoLONLAT : MccNamedCoordPair { using MccNamedCoordPair::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 concept mcc_skypoint_c = std::derived_from && requires(T t) { typename T::meteo_t; { T::meteo(std::declval()) }; }; /* MCC-LIBRARY DEFAULT GENERIC SKY POINT CLASS IMPLEMENTATION */ template class MccGenericSkyPoint : public mcc_skypoint_interface_t { public: typedef CCTE_T ccte_t; static constexpr double MJD0 = 2400000.5; inline static CCTE_T cctEngine{}; // celestial coordinates transformation engine MccGenericSkyPoint() {} template MccGenericSkyPoint(const PT& coord_pair) : MccGenericSkyPoint() { auto self = from(coord_pair); } MccGenericSkyPoint(const MccGenericSkyPoint&) = default; MccGenericSkyPoint(MccGenericSkyPoint&&) = default; MccGenericSkyPoint& operator=(const MccGenericSkyPoint&) = default; MccGenericSkyPoint& operator=(MccGenericSkyPoint&&) = default; virtual ~MccGenericSkyPoint() = default; MccCelestialCoordEpoch epoch() const { return _epoch; } template MccGenericSkyPoint& 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()); } return *this; } MccGenericSkyPoint& operator=(mcc_coord_pair_c auto const& coord_pair) { return from(coord_pair); } template auto to(PT& cpair, PTs&... cpairs) const { toHelper(cpair); if constexpr (sizeof...(PTs)) { to(cpairs...); } } template operator PT() { PT res; to(res); return res; } protected: double _x{0.0}, _y{0.0}; MccCoordPairKind _pairKind{MccCoordPairKind::COORDS_KIND_RADEC_ICRS}; MccCelestialCoordEpoch _epoch{}; // J2000.0 template auto toHelper(PT& cpair) const { 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); }; typename CCTE_T::error_t ccte_err; double phi = cctEngine.getStateERFA().lat; 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!"); if (_pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS && PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { // from ICRS to ICRS - just copy and exit cpair = PT(typename PT::x_t(_x), typename PT::y_t(_y)); return; } // just copy coordinates and exit if (_pairKind == PT::pairKind && utils::isEqual(_epoch.MJD(), cpair.MJD())) { // cpair = PT(typename PT::x_t(_x), typename PT::y_t(_y), _epoch); cpair.setX(_x); cpair.setY(_y); return; } // if epochs are 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; if (mcc_is_obs_coordpair(_pairKind)) { ccte_err = cctEngine.obsToICRS(_pairKind, _epoch, _x, _y, &ra_icrs, &dec_icrs); } else if (mcc_is_app_coordpair(_pairKind)) { ccte_err = cctEngine.appToICRS(_pairKind, _epoch, _x, _y, &ra_icrs, &dec_icrs); } else { // unsupported transformation!!! return; } if (ccte_err) { 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(typename PT::x_t(ra_icrs), typename PT::y_t(dec_icrs)); return; } // here, the input coordinates and stored one are at the same epoch 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; } 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 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 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); cpair.setY(zd); } else { if constexpr (mccIsAppCoordPairKind) { // correct for refraction: alt -= dz_refr double dZ; ccte_err = cctEngine.refractionCorrection(half_pi - alt, &dZ); alt -= dZ; } 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 { // unsupported transformation!!! return; } } } else if (cp_kind == MccCoordPairKind::COORDS_KIND_AZZD) { alt = half_pi - zd; if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) { cpair.setX(az); cpair.setY(alt); } else { 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 double dZ; ccte_err = cctEngine.refractionReverseCorrection(half_pi - alt, &dZ); alt += dZ; 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); } } }; comp_func(pkind); // ran transformation } }; /* MCC-LIBRARY DEFAULT SKY POINT CLASS WITH ERFA-LIBRARY BASED ENGINE */ typedef MccGenericSkyPoint MccSkyPoint; } // end namespace mcc