mountcontrol/mcc/mcc_coord.h
2026-01-16 12:23:58 +03:00

620 lines
22 KiB
C++

#pragma once
#include "mcc_angle.h"
#include "mcc_ccte_erfa_new.h"
#include "mcc_defaults.h"
#include "mcc_generics.h"
#include <erfa.h>
namespace mcc
{
template <mcc_angle_c CO_LON_T, mcc_angle_c CO_LAT_T>
class MccCoordPair
{
public:
typedef CO_LON_T x_t;
typedef CO_LAT_T y_t;
static constexpr MccCoordPairKind pairKind =
!(std::derived_from<CO_LON_T, MccAngle> ||
std::derived_from<CO_LAT_T, MccAngle>) // unknown type (possibly just double or float)
? MccCoordPairKind::COORDS_KIND_GENERIC
: (std::same_as<CO_LON_T, MccAngle> || std::same_as<CO_LAT_T, MccAngle>) // one of the types is MccAngle
? MccCoordPairKind::COORDS_KIND_GENERIC
// ICRS RA and DEC
: (std::same_as<CO_LON_T, MccAngleRA_ICRS> && std::same_as<CO_LAT_T, MccAngleDEC_ICRS>)
? MccCoordPairKind::COORDS_KIND_RADEC_ICRS
// apparent RA and DEC
: (std::same_as<CO_LON_T, MccAngleRA_APP> && std::same_as<CO_LAT_T, MccAngleDEC_APP>)
? MccCoordPairKind::COORDS_KIND_RADEC_APP
// observed RA and DEC
: (std::same_as<CO_LON_T, MccAngleRA_OBS> && std::same_as<CO_LAT_T, MccAngleDEC_OBS>)
? MccCoordPairKind::COORDS_KIND_RADEC_OBS
// apparent HA and DEC
: (std::same_as<CO_LON_T, MccAngleHA_APP> && std::same_as<CO_LAT_T, MccAngleDEC_APP>)
? MccCoordPairKind::COORDS_KIND_HADEC_APP
// observed HA and DEC
: (std::same_as<CO_LON_T, MccAngleHA_OBS> && std::same_as<CO_LAT_T, MccAngleDEC_OBS>)
? MccCoordPairKind::COORDS_KIND_HADEC_OBS
// apparent AZ and ZD
: (std::same_as<CO_LON_T, MccAngleAZ> && std::same_as<CO_LAT_T, MccAngleZD>)
? MccCoordPairKind::COORDS_KIND_AZZD
// apparent AZ and ALT
: (std::same_as<CO_LON_T, MccAngleAZ> && std::same_as<CO_LAT_T, MccAngleALT>)
? MccCoordPairKind::COORDS_KIND_AZALT
// general purpose X and Y
: (std::same_as<CO_LON_T, MccAngleX> && std::same_as<CO_LAT_T, MccAngleY>)
? MccCoordPairKind::COORDS_KIND_XY
// geographical longitude and latitude
: (std::same_as<CO_LON_T, MccAngleLON> && std::same_as<CO_LAT_T, MccAngleLAT>)
? MccCoordPairKind::COORDS_KIND_LONLAT
: MccCoordPairKind::COORDS_KIND_UNKNOWN;
template <mcc_coord_epoch_c EpT = MccCelestialCoordEpoch>
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<CO_LON_T, CO_LAT_T, MccCelestialCoordEpoch>() 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 <typename T>
concept mcc_coord_pair_c = requires {
requires mcc_angle_c<typename T::x_t>;
requires mcc_angle_c<typename T::y_t>;
requires std::derived_from<T, MccCoordPair<typename T::x_t, typename T::y_t>>;
};
/* predefined coordinate pairs */
template <mcc_angle_c CO_LON_T, mcc_angle_c CO_LAT_T>
requires(std::derived_from<CO_LON_T, MccAngle> && std::derived_from<CO_LAT_T, MccAngle>)
class MccNamedCoordPair : public MccCoordPair<CO_LON_T, CO_LAT_T>
{
public:
MccNamedCoordPair() : MccCoordPair<CO_LON_T, CO_LAT_T>(CO_LON_T{0.0}, CO_LAT_T{0.0}, MccCelestialCoordEpoch::now())
{
}
template <typename CxT, typename CyT, mcc_coord_epoch_c EpT = MccCelestialCoordEpoch>
requires(std::is_arithmetic_v<CxT> && std::is_arithmetic_v<CyT>)
MccNamedCoordPair(CxT const& x, CyT const& y, EpT const& epoch = EpT::now())
: MccCoordPair<CO_LON_T, CO_LAT_T>(CO_LON_T{(double)x}, CO_LAT_T{(double)y}, epoch)
{
}
template <mcc_coord_epoch_c EpT = MccCelestialCoordEpoch>
MccNamedCoordPair(MccAngle const& x, MccAngle const& y, EpT const& epoch = EpT::now())
: MccCoordPair<CO_LON_T, CO_LAT_T>(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<MccAngleRA_ICRS, MccAngleDEC_ICRS> {
template <typename CxT, typename CyT>
requires(std::is_arithmetic_v<CxT> && std::is_arithmetic_v<CyT>)
MccSkyRADEC_ICRS(CxT const& x, CyT const& y)
: MccNamedCoordPair<MccAngleRA_ICRS, MccAngleDEC_ICRS>(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<MccAngleRA_APP, MccAngleDEC_APP> {
using MccNamedCoordPair<MccAngleRA_APP, MccAngleDEC_APP>::MccNamedCoordPair;
};
typedef MccNamedCoordPair<MccAngleRA_OBS, MccAngleDEC_OBS> MccSkyRADEC_OBS;
// struct MccSkyRADEC_OBS : MccNamedCoordPair<MccAngleRA_OBS, MccAngleDEC_OBS> {
// using MccNamedCoordPair<MccAngleRA_OBS, MccAngleDEC_OBS>::MccNamedCoordPair;
// };
struct MccSkyHADEC_APP : MccNamedCoordPair<MccAngleHA_APP, MccAngleDEC_APP> {
using MccNamedCoordPair<MccAngleHA_APP, MccAngleDEC_APP>::MccNamedCoordPair;
};
struct MccSkyHADEC_OBS : MccNamedCoordPair<MccAngleHA_OBS, MccAngleDEC_OBS> {
using MccNamedCoordPair<MccAngleHA_OBS, MccAngleDEC_OBS>::MccNamedCoordPair;
};
struct MccSkyAZZD : MccNamedCoordPair<MccAngleAZ, MccAngleZD> {
using MccNamedCoordPair<MccAngleAZ, MccAngleZD>::MccNamedCoordPair;
};
struct MccSkyAZALT : MccNamedCoordPair<MccAngleAZ, MccAngleALT> {
using MccNamedCoordPair<MccAngleAZ, MccAngleALT>::MccNamedCoordPair;
};
struct MccGenXY : MccNamedCoordPair<MccAngleX, MccAngleY> {
using MccNamedCoordPair<MccAngleX, MccAngleY>::MccNamedCoordPair;
};
struct MccGeoLONLAT : MccNamedCoordPair<MccAngleLON, MccAngleLAT> {
using MccNamedCoordPair<MccAngleLON, MccAngleLAT>::MccNamedCoordPair;
};
struct mcc_skypoint_interface_t {
virtual ~mcc_skypoint_interface_t() = default;
// template <std::derived_from<mcc_skypoint_interface_t> SelfT, mcc_angle_c XT, mcc_angle_c YT>
// SelfT& from(this SelfT&& self, XT&& x, YT&& y)
// {
// return std::forward<SelfT>(self).from(std::forward<XT>(x), std::forward<YT>(y));
// }
template <std::derived_from<mcc_skypoint_interface_t> SelfT, mcc_coord_pair_c PT>
auto from(this SelfT&& self, PT&& cpair)
{
return std::forward<SelfT>(self).from(std::forward<PT>(cpair));
}
template <std::derived_from<mcc_skypoint_interface_t> SelfT, mcc_coord_pair_c PT>
auto operator=(this SelfT&& self, PT&& cpair)
{
return std::forward<SelfT>(self).operator=(std::forward<PT>(cpair));
}
template <std::derived_from<mcc_skypoint_interface_t> SelfT, mcc_coord_pair_c PT, mcc_coord_pair_c... PTs>
auto to(this SelfT&& self, PT& cpair, PTs&... cpairs)
{
return std::forward<SelfT>(self).to(cpair, cpairs...);
}
template <std::derived_from<mcc_skypoint_interface_t> SelfT, mcc_coord_pair_c PT>
operator PT(this SelfT&& self)
{
return std::forward<SelfT>(self).operator PT();
}
};
template <typename T>
concept mcc_skypoint_c = std::derived_from<T, mcc_skypoint_interface_t> && requires(T t) {
typename T::meteo_t;
{ T::meteo(std::declval<typename T::meteo_t const&>()) };
};
/* MCC-LIBRARY DEFAULT GENERIC SKY POINT CLASS IMPLEMENTATION */
template <typename CCTE_T>
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 <mcc_coord_pair_c PT>
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 <mcc_coord_pair_c PT>
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 <mcc_coord_pair_c PT, mcc_coord_pair_c... PTs>
auto to(PT& cpair, PTs&... cpairs) const
{
toHelper(cpair);
if constexpr (sizeof...(PTs)) {
to(cpairs...);
}
}
template <mcc_coord_pair_c PT>
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 <mcc_coord_pair_c PT>
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<PT::pairKind>) {
ccte_err = cctEngine.icrsToApp(ra_icrs, dec_icrs, cpair.epoch(), &ra, &dec, &ha, &az, &zd);
} else if constexpr (mccIsObsCoordPairKind<PT::pairKind>) {
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<PT::pairKind>) {
// 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<PT::pairKind>) { // 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<PT::pairKind>) { // 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<ccte::erfa::MccCCTE_ERFA> MccSkyPoint;
} // end namespace mcc