#pragma once /* MOUNT CONTROL COMPONENTS LIBRARY */ /* ASTROMETRY ENGINE BASED ON ERFA-LIBRARY (THREAD-SAFE FOR ENGINE STATE MANIPULATIONS) */ #include #include #include "mcc_astrom_iers.h" #include "mcc_mount_concepts.h" #include "mcc_mount_coord.h" namespace mcc::astrom::erfa { enum class MccMountAstromEngineERFAErrorCode : int { ERROR_OK = 0, ERROR_INVALID_INPUT_ARG, ERROR_JULDATE_INVALID_YEAR, ERROR_JULDATE_INVALID_MONTH, 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, }; } namespace std { template <> class is_error_code_enum : public true_type { }; } // namespace std namespace mcc::astrom::erfa { #include #include /* error category definition */ // error category struct MccMountAstromEngineERFACategory : public std::error_category { MccMountAstromEngineERFACategory() : std::error_category() {} const char* name() const noexcept { return "ADC_GENERIC_DEVICE"; } std::string message(int ec) const { MccMountAstromEngineERFAErrorCode err = static_cast(ec); switch (err) { case MccMountAstromEngineERFAErrorCode::ERROR_OK: return "OK"; case MccMountAstromEngineERFAErrorCode::ERROR_INVALID_INPUT_ARG: return "invalid argument"; case MccMountAstromEngineERFAErrorCode::ERROR_JULDATE_INVALID_YEAR: return "invalid year number"; case MccMountAstromEngineERFAErrorCode::ERROR_JULDATE_INVALID_MONTH: return "invalid month number"; case MccMountAstromEngineERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR: return "unsupported coordinate pair"; case MccMountAstromEngineERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE: return "time point is out of range"; case MccMountAstromEngineERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE: return "time point is out of range"; case MccMountAstromEngineERFAErrorCode::ERROR_DUBIOUS_YEAR: return "dubious year"; case MccMountAstromEngineERFAErrorCode::ERROR_UNACCEPTABLE_DATE: return "unacceptable date"; case MccMountAstromEngineERFAErrorCode::ERROR_UPDATE_LEAPSECONDS: return "leap seconds update error"; case MccMountAstromEngineERFAErrorCode::ERROR_UPDATE_BULLETINA: return "bulletin A update error"; default: return "UNKNOWN"; } } static const MccMountAstromEngineERFACategory& get() { static const MccMountAstromEngineERFACategory constInst; return constInst; } }; inline std::error_code make_error_code(MccMountAstromEngineERFAErrorCode ec) { return std::error_code(static_cast(ec), MccMountAstromEngineERFACategory::get()); } class MccMountAstromEngineERFA { public: static constexpr double DEFAULT_WAVELENGTH = 0.55; // default observed wavelength in mkm typedef std::error_code error_t; /* use of the same type for representation of celestial and geodetic coordinates, celestial angles (e.g. P.A.), * and sideral time */ typedef double 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 = "00:00:00"_dms; // site latitude coord_t lon = "00:00:00"_dms; // site longitude double elev = 0.0; // site elevation (in meters) mcc::astrom::iers::MccLeapSeconds _leapSeconds{}; mcc::astrom::iers::MccIersBulletinA _bulletinA{}; }; typedef std::chrono::system_clock::time_point time_point_t; struct juldate_t { static constexpr double MJD0 = ERFA_DJM0; double mjd{51544.5}; // J2000.0 }; typedef coord_t sideral_time_t; typedef coord_t pa_t; typedef coord_t eo_t; struct refract_result_t { double refa, refb; }; MccMountAstromEngineERFA() = default; MccMountAstromEngineERFA(engine_state_t state) : _currentState(std::move(state)), _stateMutex(new std::mutex) {} MccMountAstromEngineERFA(const MccMountAstromEngineERFA&) = delete; MccMountAstromEngineERFA& operator=(const MccMountAstromEngineERFA&) = delete; MccMountAstromEngineERFA(MccMountAstromEngineERFA&& other) : _currentState(std::move(other._currentState)), _stateMutex(std::move(other._stateMutex)) {}; MccMountAstromEngineERFA& operator=(MccMountAstromEngineERFA&& other) { if (this != &other) { _currentState = std::move(other._currentState); _stateMutex = std::move(other._stateMutex); } return *this; } virtual ~MccMountAstromEngineERFA() = default; 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> auto& stream, char comment_sym = '#') { std::lock_guard lock{*_stateMutex}; if (!_currentState._leapSeconds.load(stream, comment_sym)) { return MccMountAstromEngineERFAErrorCode::ERROR_UPDATE_LEAPSECONDS; } return MccMountAstromEngineERFAErrorCode::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 MccMountAstromEngineERFAErrorCode::ERROR_UPDATE_LEAPSECONDS; } return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t updateBulletinA(std::derived_from> auto& stream, char comment_sym = '*') { std::lock_guard lock{*_stateMutex}; if (!_currentState._bulletinA.load(stream, comment_sym)) { return MccMountAstromEngineERFAErrorCode::ERROR_UPDATE_BULLETINA; } return MccMountAstromEngineERFAErrorCode::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 MccMountAstromEngineERFAErrorCode::ERROR_UPDATE_BULLETINA; } return MccMountAstromEngineERFAErrorCode::ERROR_OK; } /* time-related methods */ static time_point_t timePointNow() { return time_point_t::clock::now(); } // templated generic version template error_t greg2jul(TpT time_point, juldate_t& juldate) { using namespace std::literals::chrono_literals; auto dd = std::chrono::floor(time_point); std::chrono::year_month_day ymd{dd}; static constexpr std::chrono::year MIN_YEAR = -4799y; if (ymd.year() < MIN_YEAR) { return MccMountAstromEngineERFAErrorCode::ERROR_JULDATE_INVALID_YEAR; } if (!ymd.month().ok()) { return MccMountAstromEngineERFAErrorCode::ERROR_JULDATE_INVALID_MONTH; } int64_t im = (unsigned)ymd.month(); int64_t id = (unsigned)ymd.day(); int64_t iy = (int)ymd.year(); int64_t my = (im - 14LL) / 12LL; int64_t iypmy = iy + my; // integer part of result MJD int64_t mjd_int = (1461LL * (iypmy + 4800LL)) / 4LL + (367LL * (im - 2LL - 12LL * my)) / 12LL - (3LL * ((iypmy + 4900LL) / 100LL)) / 4LL + id - 2432076LL; using fd_t = std::chrono::duration>; // fraction of day juldate.mjd = static_cast(mjd_int) + std::chrono::duration_cast(time_point - dd).count(); return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t greg2jul(time_point_t time_point, juldate_t& juldate) { return greg2jul(time_point, juldate); } error_t terrestrialTime(juldate_t juldate, juldate_t& tt) { std::lock_guard lock{*_stateMutex}; using real_days_t = std::chrono::duration>; auto tai_utc = _currentState._leapSeconds[juldate.mjd]; if (tai_utc.has_value()) { tt.mjd += std::chrono::duration_cast(tai_utc.value()).count(); } else { return MccMountAstromEngineERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE; } auto tt_tai = _currentState._bulletinA.TT_TAI(); tt.mjd = juldate.mjd + std::chrono::duration_cast(tt_tai).count(); return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t apparentSiderTime(juldate_t juldate, sideral_time_t& gst, bool islocal = false) { // std::lock_guard lock{*_stateMutex}; using real_days_t = std::chrono::duration>; double ut1 = juldate.mjd; // double tt = juldate.mjd; { std::lock_guard lock{*_stateMutex}; auto dut1 = _currentState._bulletinA.DUT1(juldate.mjd); if (dut1.has_value()) { ut1 += std::chrono::duration_cast(dut1.value()).count(); } else { // out of range return MccMountAstromEngineERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; } } // auto tai_utc = _currentState._leapSeconds[juldate.mjd]; // if (tai_utc.has_value()) { // tt += std::chrono::duration_cast(tai_utc.value()).count(); // } else { // return ERROR_LEAPSECONDS_OUT_OF_RANGE; // } // auto tt_tai = _currentState._bulletinA.TT_TAI(); // tt += std::chrono::duration_cast(tt_tai).count(); // gst = eraGst06a(juldate.MJD0, ut1, juldate.MJD0, tt); juldate_t tt; auto err = terrestrialTime(juldate, tt); if (err != MccMountAstromEngineERFAErrorCode::ERROR_OK) { return err; } gst = eraGst06a(juldate.MJD0, ut1, juldate.MJD0, tt.mjd); if (islocal) { gst += _currentState.lon; } return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t eqOrigins(juldate_t juldate, eo_t& eo) { juldate_t tt; auto err = terrestrialTime(juldate, tt); if (err != MccMountAstromEngineERFAErrorCode::ERROR_OK) { return err; } eo = eraEo06a(tt.MJD0, tt.mjd); return MccMountAstromEngineERFAErrorCode::ERROR_OK; } /* atmospheric refraction-related methods */ error_t refraction(refract_result_t& refr) { std::lock_guard lock{*_stateMutex}; eraRefco(_currentState.meteo.pressure, _currentState.meteo.temperature, _currentState.meteo.humidity, _currentState.wavelength, &refr.refa, &refr.refb); return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t refractCorrection(const coord_t& alt, const refract_result_t& ref_params, coord_t& corr) { if (alt <= 0.0) { corr = 35.4 / 60.0 * std::numbers::pi / 180.0; // 35.4 arcminutes } else { auto tanALT = std::tan(alt); corr = ref_params.refa / tanALT + ref_params.refb / tanALT / tanALT / tanALT; } return MccMountAstromEngineERFAErrorCode::ERROR_OK; } /* coordinates conversional methods */ error_t icrs2obs(coord_t ra, coord_t dec, juldate_t juldate, coord_t& ra_app, coord_t& dec_app, coord_t& ha, coord_t& az, coord_t& alt, eo_t& eo) { std::lock_guard lock{*_stateMutex}; auto dut1 = _currentState._bulletinA.DUT1(juldate.mjd); if (!dut1.has_value()) { return MccMountAstromEngineERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; } auto pol_pos = _currentState._bulletinA.polarCoords(juldate.mjd); if (!pol_pos.has_value()) { return MccMountAstromEngineERFAErrorCode::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 ret = eraAtco13(ra, dec, 0.0, 0.0, 0.0, 0.0, juldate.MJD0, juldate.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 (ret == 1) { return MccMountAstromEngineERFAErrorCode::ERROR_DUBIOUS_YEAR; } else if (ret == -1) { return MccMountAstromEngineERFAErrorCode::ERROR_UNACCEPTABLE_DATE; } ra_app = ora; dec_app = odec; az = oaz; alt = std::numbers::pi / 2.0 - ozd; ha = oha; eo = eo_; return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t obs2icrs(MccCoordPairKind coord_kind, coord_t x, coord_t y, juldate_t juldate, coord_t ra, coord_t dec) { std::lock_guard lock{*_stateMutex}; auto dut1 = _currentState._bulletinA.DUT1(juldate.mjd); if (!dut1.has_value()) { return MccMountAstromEngineERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; } auto pol_pos = _currentState._bulletinA.polarCoords(juldate.mjd); if (!pol_pos.has_value()) { return MccMountAstromEngineERFAErrorCode::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 (coord_kind) { case mcc::MccCoordPairKind::COORDS_KIND_AZZD: type = "A"; break; case mcc::MccCoordPairKind::COORDS_KIND_AZALT: y = std::numbers::pi / 2.0 - 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 MccMountAstromEngineERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; } double ra_icrs, dec_icrs; int ret = eraAtoc13(type.c_str(), x, y, juldate.MJD0, juldate.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, &ra_icrs, &dec_icrs); if (ret == 1) { return MccMountAstromEngineERFAErrorCode::ERROR_DUBIOUS_YEAR; } else if (ret == -1) { return MccMountAstromEngineERFAErrorCode::ERROR_UNACCEPTABLE_DATE; } ra = ra_icrs; dec = dec_icrs; return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t hadec2azalt(coord_t ha, coord_t dec, coord_t& az, coord_t& alt) { std::lock_guard lock{*_stateMutex}; double r_az, r_alt; eraHd2ae(ha, dec, _currentState.lat, &r_az, &r_alt); az = r_az; alt = r_alt; return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t azalt2hadec(coord_t az, coord_t alt, coord_t& ha, coord_t& dec) { std::lock_guard lock{*_stateMutex}; double r_ha, r_dec; eraAe2hd(az, alt, _currentState.lat, &r_ha, &r_dec); ha = r_ha; dec = r_dec; return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t hadec2pa(coord_t ha, coord_t dec, pa_t& pa) { std::lock_guard lock{*_stateMutex}; pa = eraHd2pa(ha, dec, _currentState.lat); return MccMountAstromEngineERFAErrorCode::ERROR_OK; } error_t coord2coord(MccCoordPairKind coord_from_kind, coord_t x_from, coord_t y_from, time_point_t time_point_from, MccCoordPairKind coord_to_kind, coord_t& x_to, coord_t& y_to, time_point_t time_point_to) { error_t ret = MccMountAstromEngineERFAErrorCode::ERROR_OK; // no convertion if (time_point_from == time_point_to && coord_from_kind == coord_to_kind) { x_to = x_from; y_to = y_from; return ret; } juldate_t jd; eo_t eo; coord_t ra, dec, ha, az, alt; sideral_time_t lst; auto lst_eo = [&]() { ret = greg2jul(time_point_from, jd); if (!ret) { ret = apparentSiderTime(jd, lst, true); if (!ret) { ret = eqOrigins(jd, eo); } } }; // special case: to ICRS from apparent if (coord_to_kind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { if (coord_to_kind == MccCoordPairKind::COORDS_KIND_AZALT || coord_to_kind == MccCoordPairKind::COORDS_KIND_AZZD || coord_to_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP || coord_to_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { // ret = greg2jul(time_point_from, jd); if (!ret) { ret = obs2icrs(coord_from_kind, x_from, y_from, jd, x_to, y_to); } } else { ret = MccMountAstromEngineERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; } return ret; } // special case: from ICRS to apparent if (coord_from_kind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { if (coord_to_kind == MccCoordPairKind::COORDS_KIND_AZALT || coord_to_kind == MccCoordPairKind::COORDS_KIND_AZZD || coord_to_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP || coord_to_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { // ret = greg2jul(time_point_to, jd); if (!ret) { ret = icrs2obs(x_from, y_from, jd, ra, dec, ha, az, alt, eo); if (!ret) { if (coord_to_kind == MccCoordPairKind::COORDS_KIND_AZALT) { x_to = az; y_to = alt; } else if (coord_to_kind == MccCoordPairKind::COORDS_KIND_AZZD) { x_to = az; y_to = std::numbers::pi / 2.0 - alt; } else if (coord_to_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { x_to = ra; y_to = dec; } else if (coord_to_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { x_to = ha; y_to = dec; } } } } else { ret = MccMountAstromEngineERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; } return ret; } // from apparent to apparent if (time_point_from != time_point_to) { // first, convert 'from' coordinates to ICRS ret = greg2jul(time_point_from, jd); if (!ret) { ret = obs2icrs(coord_from_kind, x_from, y_from, jd, ra, dec); if (!ret) { // from ICRS to required return coord2coord(MccCoordPairKind::COORDS_KIND_RADEC_ICRS, ra, dec, time_point_from, coord_to_kind, x_to, y_to, time_point_to); } } } else { // same time points if (coord_from_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { // first, compute HA from CIO-based RA!!! lst_eo(); if (!ret) { ha = lst - x_from + eo; } else { return ret; } if (coord_to_kind == MccCoordPairKind::COORDS_KIND_AZALT) { ret = hadec2azalt(ha, y_from, x_to, y_to); } else if (coord_to_kind == MccCoordPairKind::COORDS_KIND_AZZD) { ret = hadec2azalt(ha, y_from, x_to, y_to); if (!ret) { y_to = std::numbers::pi / 2.0 - y_to; } } else if (coord_to_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { x_to = ha; y_to = y_from; } else { ret = MccMountAstromEngineERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; } } else if (coord_from_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { if (coord_to_kind == MccCoordPairKind::COORDS_KIND_AZALT) { ret = hadec2azalt(x_from, y_from, x_to, y_to); } else if (coord_to_kind == MccCoordPairKind::COORDS_KIND_AZZD) { ret = hadec2azalt(x_from, y_from, x_to, y_to); if (!ret) { y_to = std::numbers::pi / 2.0 - y_to; } } else if (coord_to_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { lst_eo(); if (!ret) { x_to = lst - x_from + eo; y_to = y_from; } } else { ret = MccMountAstromEngineERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; } } else if (coord_from_kind == MccCoordPairKind::COORDS_KIND_AZALT) { if (coord_to_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { ret = azalt2hadec(x_from, y_from, x_to, y_to); if (!ret) { lst_eo(); if (!ret) { x_to = lst - x_to + eo; } } } else if (coord_to_kind == MccCoordPairKind::COORDS_KIND_AZZD) { x_to = x_from; y_to = std::numbers::pi / 2.0 - y_from; } else if (coord_to_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) { ret = azalt2hadec(x_from, y_from, x_to, y_to); } else { ret = MccMountAstromEngineERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; } } else if (coord_from_kind == MccCoordPairKind::COORDS_KIND_AZZD) { return coord2coord(MccCoordPairKind::COORDS_KIND_AZALT, std::move(x_from), std::numbers::pi / 2.0 - y_to, time_point_from, coord_to_kind, x_to, y_to, time_point_to); } else { ret = MccMountAstromEngineERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR; } } 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 _stateMutex; }; } // namespace mcc::astrom::erfa static_assert(mcc::traits::mcc_astrom_engine_c, "");