874 lines
28 KiB
C++
874 lines
28 KiB
C++
#pragma once
|
|
|
|
/* MOUNT CONTROL COMPONENTS LIBRARY */
|
|
|
|
|
|
/* CELESTIAL COORDINATES TRANSFORMATION ENGINE IMPLEMENTATION BASED ON ERFA LIBRARY */
|
|
|
|
|
|
// NOTE: according to definition of astronomical azimuth it is counted from the South through the West, but
|
|
// in the ERFA the azimuth is counted from the North through the East!!!
|
|
//
|
|
// The implementation expects input azimuth as the "South"-defined one and transforms ERFA-outputs
|
|
// to the "South"-defined one respectively
|
|
|
|
|
|
|
|
#include <erfa.h>
|
|
#include <erfam.h>
|
|
#include <mutex>
|
|
|
|
|
|
#include "mcc_angle.h"
|
|
#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 mcc_CCTE_interface_t<std::error_code>
|
|
{
|
|
static constexpr double PI_2 = std::numbers::pi / 2.0;
|
|
|
|
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) {}
|
|
|
|
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;
|
|
|
|
|
|
std::string_view name() const
|
|
{
|
|
return "ERFA-CCTE-ENGINE";
|
|
}
|
|
|
|
// 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)
|
|
{
|
|
error_t 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((int)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)
|
|
{
|
|
error_t 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) {
|
|
//
|
|
} else {
|
|
ret = MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR;
|
|
}
|
|
|
|
ret = obs2icrs(from_pt, to_pt);
|
|
|
|
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) {
|
|
//
|
|
MccEqtHrzCoords eqhrz;
|
|
|
|
ret = icrs2obs(from_pt, &eqhrz);
|
|
|
|
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;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
// from apparent to apparent
|
|
if (from_pt.time_point != to_pt->time_point) { // first, convert 'from' coordinates to ICRS
|
|
MccCelestialPoint pt{.pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_ICRS};
|
|
pt.time_point = to_pt->time_point;
|
|
|
|
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) {
|
|
from_pt.X = ha;
|
|
hadec2azalt(from_pt, to_pt);
|
|
to_pt->Y = PI_2 - to_pt->Y;
|
|
} else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZALT) {
|
|
from_pt.X = ha;
|
|
hadec2azalt(from_pt, to_pt);
|
|
} 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) {
|
|
hadec2azalt(from_pt, to_pt);
|
|
to_pt->Y = PI_2 - to_pt->Y;
|
|
} else if (to_pt->pair_kind == MccCoordPairKind::COORDS_KIND_AZALT) {
|
|
hadec2azalt(from_pt, to_pt);
|
|
} 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) {
|
|
azalt2hadec(from_pt, to_pt);
|
|
} 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) {
|
|
azalt2hadec(from_pt, to_pt);
|
|
lst_eo();
|
|
if (!ret) {
|
|
to_pt->X = lst - to_pt->X + eo;
|
|
}
|
|
} 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;
|
|
}
|
|
|
|
|
|
// 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;
|
|
|
|
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;
|
|
}
|
|
|
|
|
|
// NOTE: according to definition of astronomical azimuth it is counted from the South through the West, but
|
|
// in the ERFA the azimuth is counted from the North through the East!!!
|
|
//
|
|
result->AZ = MccAngle(az + std::numbers::pi).normalize<MccAngle::NORM_KIND_0_360>();
|
|
result->ZD = zd;
|
|
result->HA = ha;
|
|
result->RA_APP = ra;
|
|
result->DEC_APP = dec;
|
|
result->ALT = MccCCTE_ERFA::PI_2 - zd;
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
error_t obs2icrs(mcc_celestial_point_c auto 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:
|
|
// NOTE: according to definition of astronomical azimuth it is counted from the South through the West,
|
|
// but in the ERFA the azimuth is counted from the North through the East!!!
|
|
//
|
|
from_pt.X += std::numbers::pi;
|
|
type = "A";
|
|
break;
|
|
case mcc::MccCoordPairKind::COORDS_KIND_AZALT:
|
|
// NOTE: according to definition of astronomical azimuth it is counted from the South through the West,
|
|
// but in the ERFA the azimuth is counted from the North through the East!!!
|
|
//
|
|
from_pt.X += std::numbers::pi;
|
|
from_pt.Y = MccCCTE_ERFA::PI_2 - 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) {
|
|
ret = MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR;
|
|
} else if (err == -1) {
|
|
ret = MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
void azalt2hadec(mcc_celestial_point_c auto const& from_pt, mcc_celestial_point_c auto* to_pt)
|
|
{
|
|
double ha, dec;
|
|
|
|
// NOTE: according to definition of astronomical azimuth it is counted from the South through the West, but
|
|
// in the ERFA the azimuth is counted from the North through the East!!!
|
|
//
|
|
eraAe2hd(from_pt.X + std::numbers::pi, from_pt.Y, _currentState.lat, &ha, &dec);
|
|
|
|
to_pt->X = ha;
|
|
to_pt->Y = dec;
|
|
}
|
|
|
|
|
|
void hadec2azalt(mcc_celestial_point_c auto const& from_pt, mcc_celestial_point_c auto* to_pt)
|
|
{
|
|
double az, alt;
|
|
|
|
eraHd2ae(from_pt.X, from_pt.Y, _currentState.lat, &az, &alt);
|
|
|
|
// NOTE: according to definition of astronomical azimuth it is counted from the South through the West, but
|
|
// in the ERFA the azimuth is counted from the North through the East!!!
|
|
//
|
|
// convert it to "from the South through the West"
|
|
to_pt->X = MccAngle(az - std::numbers::pi).normalize<MccAngle::NORM_KIND_0_360>();
|
|
to_pt->Y = alt;
|
|
}
|
|
|
|
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;
|
|
}
|
|
};
|
|
|
|
static_assert(mcc_ccte_c<MccCCTE_ERFA>, "");
|
|
|
|
} // namespace mcc::ccte::erfa
|