This commit is contained in:
2026-01-15 19:11:16 +03:00
parent 09cf5f9c19
commit 01d5657b1b
8 changed files with 867 additions and 637 deletions

View File

@@ -1,5 +1,6 @@
#pragma once
#include "mcc_angle.h"
#include "mcc_ccte_iers.h"
#include "mcc_defaults.h"
#include "mcc_generics.h"
@@ -42,7 +43,7 @@ public:
? 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_AZZD
? 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
@@ -52,7 +53,7 @@ public:
: 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), _mjd(epoch.MJD)
MccCoordPair(CO_LON_T const& x, CO_LAT_T const& y, EpT const& epoch = EpT::now()) : _x(x), _y(y), _epoch(epoch)
{
}
@@ -74,16 +75,22 @@ public:
return _y;
}
MccCelestialCoordEpoch epoch() const
{
return _epoch;
}
double MJD() const
{
return _mjd;
return _epoch.MJD();
}
// for something like:
// auto [ra, dec, mjd] = coord_pair;
operator std::tuple<CO_LON_T, CO_LAT_T, double>() const
// auto [ra, dec, epoch] = coord_pair;
operator std::tuple<CO_LON_T, CO_LAT_T, MccCelestialCoordEpoch>() const
{
return {_x, _y, _mjd};
return {_x, _y, _epoch};
}
void setX(const CO_LON_T& x)
@@ -96,16 +103,16 @@ public:
_y = y;
}
void setMJD(double mjd)
void setEpoch(mcc_coord_epoch_c auto const& ep)
{
_mjd = mjd;
_epoch = ep;
}
protected:
CO_LON_T _x;
CO_LAT_T _y;
double _mjd;
MccCelestialCoordEpoch _epoch;
};
@@ -128,7 +135,7 @@ public:
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{x}, CO_LAT_T{y}, epoch)
: MccCoordPair<CO_LON_T, CO_LAT_T>(CO_LON_T{(double)x}, CO_LAT_T{(double)y}, epoch)
{
}
@@ -148,15 +155,62 @@ public:
virtual ~MccNamedCoordPair() = default;
};
using MccSkyRADEC_ICRS = MccNamedCoordPair<MccAngleRA_ICRS, MccAngleDEC_ICRS>;
using MccSkyRADEC_APP = MccNamedCoordPair<MccAngleRA_APP, MccAngleDEC_APP>;
using MccSkyRADEC_OBS = MccNamedCoordPair<MccAngleRA_OBS, MccAngleDEC_OBS>;
using MccSkyHADEC_APP = MccNamedCoordPair<MccAngleHA_APP, MccAngleDEC_APP>;
using MccSkyHADEC_OBS = MccNamedCoordPair<MccAngleHA_OBS, MccAngleDEC_OBS>;
using MccSkyAZZD = MccNamedCoordPair<MccAngleAZ, MccAngleZD>;
using MccSkyAZALT = MccNamedCoordPair<MccAngleAZ, MccAngleALT>;
using MccGenXY = MccNamedCoordPair<MccAngleX, MccAngleY>;
using MccGeoLONLAT = MccNamedCoordPair<MccAngleLON, MccAngleLAT>;
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& ep)
{
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;
};
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 {
@@ -211,69 +265,31 @@ concept mcc_skypoint_c = std::derived_from<T, mcc_skypoint_interface_t> && requi
/* MCC-LIBRARY DEFAULT SKY POINT CLASS IMPLEMENTATION BASED ON ERFA-LIBRARY */
class MccSkyPoint : public mcc_skypoint_interface_t
template <typename CCTE_T>
class MccGenericSkyPoint : public mcc_skypoint_interface_t
{
public:
typedef CCTE_T ccte_t;
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
};
inline static CCTE_T cctEngine{}; // celestial coordinates transformation engine
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() {}
MccGenericSkyPoint() {}
template <mcc_coord_pair_c PT>
MccSkyPoint(const PT& coord_pair) : MccSkyPoint()
MccGenericSkyPoint(const PT& coord_pair) : MccGenericSkyPoint()
{
auto self = from(coord_pair);
}
MccSkyPoint(const MccSkyPoint&) = default;
MccSkyPoint(MccSkyPoint&&) = default;
MccGenericSkyPoint(const MccGenericSkyPoint&) = default;
MccGenericSkyPoint(MccGenericSkyPoint&&) = default;
MccSkyPoint& operator=(const MccSkyPoint&) = default;
MccSkyPoint& operator=(MccSkyPoint&&) = default;
MccGenericSkyPoint& operator=(const MccGenericSkyPoint&) = default;
MccGenericSkyPoint& operator=(MccGenericSkyPoint&&) = default;
virtual ~MccSkyPoint() = default;
virtual ~MccGenericSkyPoint() = default;
MccCelestialCoordEpoch epoch() const
{
@@ -281,7 +297,7 @@ public:
}
template <mcc_coord_pair_c PT>
MccSkyPoint& from(const PT& coord_pair)
MccGenericSkyPoint& from(const PT& coord_pair)
{
_x = coord_pair.x();
_y = coord_pair.y();
@@ -293,9 +309,11 @@ public:
} else {
_epoch.fromMJD(coord_pair.MJD());
}
return *this;
}
MccSkyPoint& operator=(mcc_coord_pair_c auto const& coord_pair)
MccGenericSkyPoint& operator=(mcc_coord_pair_c auto const& coord_pair)
{
return from(coord_pair);
}
@@ -312,18 +330,6 @@ public:
}
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
@@ -331,10 +337,6 @@ protected:
template <mcc_coord_pair_c PT>
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)
@@ -385,27 +387,62 @@ protected:
dec = std::atan2(z, r);
};
double phi;
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!");
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;
if (_pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
if constexpr (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;
} else { // from ICRS to apparent or observed
if constexpr (mccIsAppCoordPairKind<PT::pairKind>) {
ccte_err = cctEngine.icrsToApp(_x, _y, cpair.epoch(), &ra, &dec, &ha, &az, &zd);
} else if constexpr (mccIsObsCoordPairKind<PT::pairKind>) {
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!");
}
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);
// cpair = PT(typename PT::x_t(_x), typename PT::y_t(_y), _epoch);
cpair.setX(_x);
cpair.setY(_y);
return;
}
// if epochs is not the same then
// 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;
@@ -413,16 +450,17 @@ protected:
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(...)
ccte_err = cctEngine.obsToICRS(_pairKind, _epoch, _x, _y, &ra_icrs, &dec_icrs);
} else if (mcc_is_app_coordpair(_pairKind)) {
// cct_engine.appToICRS(...)
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;
@@ -432,119 +470,14 @@ protected:
// 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{});
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
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<PT::pairKind>) {
// cct_engine.icrsToObs(...)
} else if constexpr (mccIsAppCoordPairKind<PT::pairKind>) {
// 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!!!
ccte_err = cctEngine.equationOrigins(cpair.MJD(), &eo);
if (ccte_err) {
return;
}
@@ -561,11 +494,15 @@ 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;
}
auto comp_func = [&, this](this auto& self, MccCoordPairKind cp_kind) {
// coordinate transformation lambda (possibly recursive!!!)
auto comp_func = [&, this](this auto&& self, MccCoordPairKind cp_kind) -> void {
if (cp_kind == MccCoordPairKind::COORDS_KIND_AZALT) {
if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) {
zd = half_pi - alt;
@@ -574,6 +511,9 @@ protected:
} 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);
@@ -629,6 +569,9 @@ protected:
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);
}
}
@@ -650,6 +593,8 @@ protected:
}
}
};
comp_func(pkind); // ran transformation
}
};