542 lines
17 KiB
C++
542 lines
17 KiB
C++
#pragma once
|
|
|
|
/* MOUNT CONTROL COMPONENTS LIBRARY */
|
|
|
|
|
|
/* CELESTIAL COORDINATES TRANSFORMATION ENGINE IMPLEMENTATION BASED ON ERFA LIBRARY */
|
|
|
|
#include <erfa.h>
|
|
#include <erfam.h>
|
|
#include <mutex>
|
|
|
|
|
|
#include "mcc_ccte_iers.h"
|
|
#include "mcc_defaults.h"
|
|
|
|
namespace mcc::ccte::erfa
|
|
{
|
|
|
|
|
|
enum class MccCCTE_ERFAErrorCode : int {
|
|
ERROR_OK = 0,
|
|
ERROR_INVALID_INPUT_ARG,
|
|
ERROR_julday_INVALID_YEAR,
|
|
ERROR_julday_INVALID_MONTH,
|
|
ERROR_julday_INVALID_DAY,
|
|
ERROR_UNSUPPORTED_COORD_PAIR,
|
|
ERROR_BULLETINA_OUT_OF_RANGE,
|
|
ERROR_LEAPSECONDS_OUT_OF_RANGE,
|
|
ERROR_DUBIOUS_YEAR,
|
|
ERROR_UNACCEPTABLE_DATE,
|
|
ERROR_UPDATE_LEAPSECONDS,
|
|
ERROR_UPDATE_BULLETINA,
|
|
ERROR_UNEXPECTED
|
|
};
|
|
|
|
} // namespace mcc::ccte::erfa
|
|
|
|
|
|
namespace std
|
|
{
|
|
|
|
template <>
|
|
class is_error_code_enum<mcc::ccte::erfa::MccCCTE_ERFAErrorCode> : public true_type
|
|
{
|
|
};
|
|
|
|
} // namespace std
|
|
|
|
|
|
|
|
namespace mcc::ccte::erfa
|
|
{
|
|
|
|
|
|
/* error category definition */
|
|
|
|
// error category
|
|
struct MccCCTE_ERFACategory : public std::error_category {
|
|
MccCCTE_ERFACategory() : std::error_category() {}
|
|
|
|
const char* name() const noexcept
|
|
{
|
|
return "CCTE-ERFA";
|
|
}
|
|
|
|
std::string message(int ec) const
|
|
{
|
|
MccCCTE_ERFAErrorCode err = static_cast<MccCCTE_ERFAErrorCode>(ec);
|
|
|
|
switch (err) {
|
|
case MccCCTE_ERFAErrorCode::ERROR_OK:
|
|
return "OK";
|
|
case MccCCTE_ERFAErrorCode::ERROR_INVALID_INPUT_ARG:
|
|
return "invalid argument";
|
|
case MccCCTE_ERFAErrorCode::ERROR_julday_INVALID_YEAR:
|
|
return "invalid year number";
|
|
case MccCCTE_ERFAErrorCode::ERROR_julday_INVALID_MONTH:
|
|
return "invalid month number";
|
|
case MccCCTE_ERFAErrorCode::ERROR_julday_INVALID_DAY:
|
|
return "invalid day number";
|
|
case MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR:
|
|
return "unsupported coordinate pair";
|
|
case MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE:
|
|
return "time point is out of range";
|
|
case MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE:
|
|
return "time point is out of range";
|
|
case MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR:
|
|
return "dubious year";
|
|
case MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE:
|
|
return "unacceptable date";
|
|
case MccCCTE_ERFAErrorCode::ERROR_UPDATE_LEAPSECONDS:
|
|
return "leap seconds update error";
|
|
case MccCCTE_ERFAErrorCode::ERROR_UPDATE_BULLETINA:
|
|
return "bulletin A update error";
|
|
case MccCCTE_ERFAErrorCode::ERROR_UNEXPECTED:
|
|
return "unexpected error value";
|
|
default:
|
|
return "UNKNOWN";
|
|
}
|
|
}
|
|
|
|
static const MccCCTE_ERFACategory& get()
|
|
{
|
|
static const MccCCTE_ERFACategory constInst;
|
|
return constInst;
|
|
}
|
|
};
|
|
|
|
|
|
inline std::error_code make_error_code(MccCCTE_ERFAErrorCode ec)
|
|
{
|
|
return std::error_code(static_cast<int>(ec), MccCCTE_ERFACategory::get());
|
|
}
|
|
|
|
|
|
|
|
class MccCCTE_ERFA
|
|
{
|
|
public:
|
|
static constexpr double DEFAULT_WAVELENGTH = 0.55; // default observed wavelength in mkm
|
|
|
|
typedef std::error_code error_t;
|
|
|
|
struct refract_model_t {
|
|
static constexpr std::string_view name()
|
|
{
|
|
return "ERFA";
|
|
}
|
|
|
|
double refa, refb;
|
|
};
|
|
|
|
/* use of the same type for representation of celestial and geodetic coordinates, celestial angles (e.g. P.A.),
|
|
* and sideral time */
|
|
typedef MccCelestialPoint::coord_t coord_t;
|
|
|
|
// meteo parameters (to compute refraction)
|
|
struct meteo_t {
|
|
typedef double temp_t;
|
|
typedef double humid_t;
|
|
typedef double press_t;
|
|
|
|
temp_t temperature; // Temperature in C
|
|
humid_t humidity; // humidity in % ([0.0, 1.0])
|
|
press_t pressure; // atmospheric presure in hPa=mB
|
|
};
|
|
|
|
struct engine_state_t {
|
|
meteo_t meteo{.temperature = 0.0, .humidity = 0.5, .pressure = 1010.0};
|
|
|
|
double wavelength = DEFAULT_WAVELENGTH; // observed wavelength in mkm
|
|
|
|
coord_t lat = 0.0; // site latitude
|
|
coord_t lon = 0.0; // site longitude
|
|
double elev = 0.0; // site elevation (in meters)
|
|
|
|
mcc::ccte::iers::MccLeapSeconds _leapSeconds{};
|
|
mcc::ccte::iers::MccIersBulletinA _bulletinA{};
|
|
};
|
|
|
|
MccCCTE_ERFA() : _stateMutex(new std::mutex) {}
|
|
|
|
virtual ~MccCCTE_ERFA() = default;
|
|
|
|
// engine state related methods
|
|
|
|
void setState(engine_state_t state)
|
|
{
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
_currentState = std::move(state);
|
|
}
|
|
|
|
engine_state_t getState() const
|
|
{
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
return _currentState;
|
|
}
|
|
|
|
void updateMeteo(meteo_t meteo)
|
|
{
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
_currentState.meteo = std::move(meteo);
|
|
}
|
|
|
|
error_t updateLeapSeconds(std::derived_from<std::basic_istream<char>> auto& stream, char comment_sym = '#')
|
|
{
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
if (!_currentState._leapSeconds.load(stream, comment_sym)) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_UPDATE_LEAPSECONDS;
|
|
}
|
|
|
|
return MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
}
|
|
|
|
|
|
error_t updateLeapSeconds(traits::mcc_input_char_range auto const& filename, char comment_sym = '#')
|
|
{
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
if (!_currentState._leapSeconds.load(filename, comment_sym)) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_UPDATE_LEAPSECONDS;
|
|
}
|
|
|
|
return MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
}
|
|
|
|
|
|
error_t updateBulletinA(std::derived_from<std::basic_istream<char>> auto& stream, char comment_sym = '*')
|
|
{
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
if (!_currentState._bulletinA.load(stream, comment_sym)) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_UPDATE_BULLETINA;
|
|
}
|
|
|
|
return MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
}
|
|
|
|
|
|
error_t updateBulletinA(traits::mcc_input_char_range auto const& filename, char comment_sym = '*')
|
|
{
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
if (!_currentState._bulletinA.load(filename, comment_sym)) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_UPDATE_BULLETINA;
|
|
}
|
|
|
|
return MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
}
|
|
|
|
|
|
// time-related methods
|
|
|
|
error_t timepointToJulday(mcc_time_point_c auto tp, mcc_julday_c auto* julday)
|
|
{
|
|
auto ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
if (julday == nullptr) {
|
|
return ret;
|
|
}
|
|
|
|
auto days = std::chrono::floor<std::chrono::days>(tp);
|
|
std::chrono::year_month_day ymd{days};
|
|
|
|
double mjd0;
|
|
|
|
int err = eraCal2jd(ymd.year(), (unsigned)ymd.month(), (unsigned)ymd.day(), &mjd0, &julday->mjd);
|
|
|
|
if (err != 0) {
|
|
ret = err == -1 ? MccCCTE_ERFAErrorCode::ERROR_julday_INVALID_YEAR
|
|
: err == -2 ? MccCCTE_ERFAErrorCode::ERROR_julday_INVALID_MONTH
|
|
: err == -3 ? MccCCTE_ERFAErrorCode::ERROR_julday_INVALID_DAY
|
|
: MccCCTE_ERFAErrorCode::ERROR_UNEXPECTED;
|
|
} else { // partial part of day
|
|
julday->mjd +=
|
|
std::chrono::duration_cast<std::chrono::duration<double, std::ratio<86400>>>(tp - days).count();
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
error_t timepointToAppSideral(mcc_time_point_c auto tp, mcc_angle_c auto* st, bool islocal = false)
|
|
{
|
|
auto ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
if (st == nullptr) {
|
|
return ret;
|
|
}
|
|
|
|
using real_days_t = std::chrono::duration<double, std::ratio<86400>>;
|
|
|
|
MccJulianDay julday;
|
|
|
|
ret = timepointToJulday(tp, &julday);
|
|
if (ret != MccCCTE_ERFAErrorCode::ERROR_OK) {
|
|
return ret;
|
|
}
|
|
|
|
double ut1 = julday.mjd;
|
|
double tt = julday.mjd;
|
|
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
auto dut1 = _currentState._bulletinA.DUT1(julday.mjd);
|
|
|
|
if (dut1.has_value()) {
|
|
ut1 += std::chrono::duration_cast<real_days_t>(dut1.value()).count();
|
|
} else { // out of range
|
|
return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
|
|
}
|
|
|
|
auto tai_utc = _currentState._leapSeconds[julday.mjd];
|
|
if (tai_utc.has_value()) {
|
|
tt += std::chrono::duration_cast<real_days_t>(tai_utc.value()).count();
|
|
} else {
|
|
return MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE;
|
|
}
|
|
|
|
|
|
auto tt_tai = _currentState._bulletinA.TT_TAI();
|
|
tt += std::chrono::duration_cast<real_days_t>(tt_tai).count();
|
|
|
|
*st = eraGst06a(julday.MJD0, ut1, julday.MJD0, tt);
|
|
|
|
if (islocal) {
|
|
*st += _currentState.lon;
|
|
}
|
|
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
// coordinates transformations
|
|
|
|
error_t transformCoordinates(mcc_celestial_point_c auto from_pt, mcc_celestial_point_c auto* to_pt)
|
|
{
|
|
error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
MccJulianDay jd;
|
|
|
|
if (to_pt == nullptr) {
|
|
return ret;
|
|
}
|
|
|
|
// no transformations
|
|
if (from_pt.time_point == to_pt->time_point && from_pt.pair_kind == to_pt->pair_kind) {
|
|
to_pt->X = from_pt.X;
|
|
to_pt->Y = from_pt.Y;
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
// special case: to ICRS from apparent
|
|
if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
|
|
if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_AZALT ||
|
|
from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_AZZD ||
|
|
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;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
// special case: from ICRS to apparent
|
|
if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
|
|
if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZALT ||
|
|
to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZZD ||
|
|
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;
|
|
}
|
|
|
|
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;
|
|
|
|
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;
|
|
}
|
|
} else {
|
|
ret = MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
}
|
|
|
|
|
|
// refraction
|
|
|
|
error_t refractionModel(refract_model_t* model)
|
|
{
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
eraRefco(_currentState.meteo.pressure, _currentState.meteo.temperature, _currentState.meteo.humidity,
|
|
_currentState.wavelength, &model->refa, &model->refb);
|
|
|
|
return MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
}
|
|
|
|
error_t refractionCorrection(mcc_celestial_point_c auto pt, mcc_angle_c auto* dZ)
|
|
{
|
|
error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
if (dZ == nullptr) {
|
|
return ret;
|
|
}
|
|
|
|
refract_model_t rmodel;
|
|
|
|
refractionModel(&rmodel);
|
|
|
|
if (pt.pair_kind == MccCoordPairKind::COORDS_KIND_AZZD) {
|
|
if (pt.Y >= std::numbers::pi / 2.0) {
|
|
*dZ = 35.4 / 60.0 * std::numbers::pi / 180.0; // 35.4 arcminutes
|
|
} else {
|
|
auto tanZ = tan(pt.Y);
|
|
*dZ = rmodel.refa * tanZ + rmodel.refb * tanZ * tanZ * tanZ;
|
|
}
|
|
} else {
|
|
MccCelestialPoint tr_pt{.pair_kind = MccCoordPairKind::COORDS_KIND_AZZD, .time_point = pt.time_point};
|
|
ret = transformCoordinates(std::move(pt), &tr_pt);
|
|
if (!ret) {
|
|
ret = refractionCorrection(std::move(tr_pt), dZ);
|
|
}
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
/* helper mathods */
|
|
|
|
auto leapSecondsExpireDate() const
|
|
{
|
|
return _currentState._leapSeconds.expireDate();
|
|
}
|
|
|
|
auto leapSecondsExpireMJD() const
|
|
{
|
|
return _currentState._leapSeconds.expireMJD();
|
|
}
|
|
|
|
|
|
auto bulletinADateRange() const
|
|
{
|
|
return _currentState._bulletinA.dateRange();
|
|
}
|
|
|
|
auto bulletinADateRangeMJD() const
|
|
{
|
|
return _currentState._bulletinA.dateRangeMJD();
|
|
}
|
|
|
|
|
|
protected:
|
|
engine_state_t _currentState{};
|
|
|
|
std::unique_ptr<std::mutex> _stateMutex;
|
|
};
|
|
|
|
} // namespace mcc::ccte::erfa
|