eddy/mountcontrol
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

...

Browse Source
This commit is contained in:
Timur A. Fatkhullin
2025-08-18 12:22:45 +03:00
parent 6c10c6b6ff
commit ca74d1dd73
5 changed files with 1007 additions and 116 deletions
Download Patch File Download Diff File Expand all files Collapse all files

468
mcc/mcc_ccte_erfa.h
View File

@@ -10,6 +10,7 @@
#include <mutex>
#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<decltype(*to_pt)>;
MccGenericCelestialPoint<typename tp_pt_t::coord_t> 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<std::mutex> _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<double, std::ratio<86400>>;
auto tai_utc = _currentState._leapSeconds[jd.mjd];
if (tai_utc.has_value()) {
tt.mjd = jd.mjd + std::chrono::duration_cast<real_days_t>(tai_utc.value()).count();
auto tt_tai = _currentState._bulletinA.TT_TAI();
tt.mjd = jd.mjd + std::chrono::duration_cast<real_days_t>(tt_tai).count();
*eo = eraEo06a(tt.MJD0, tt.mjd);
} else {
ret = MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE;
}
return ret;
}
};
} // namespace mcc::ccte::erfa