#pragma once #include #include #include #include "mcc_ccte_iers.h" #include "mcc_coord.h" #include "mcc_generics.h" namespace mcc::ccte::erfa { enum class MccCCTE_ERFAErrorCode : int { ERROR_OK = 0, ERROR_NULLPTR, 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 : 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(ec); switch (err) { case MccCCTE_ERFAErrorCode::ERROR_OK: return "OK"; case MccCCTE_ERFAErrorCode::ERROR_NULLPTR: return "input argument is the nullptr"; 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(ec), MccCCTE_ERFACategory::get()); } 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 typedef std::error_code error_t; struct refract_model_t { static constexpr std::string_view name() { return "ERFA"; } double refa, refb; }; // 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 }; // celestial object addition parameters struct obj_pars_t { double pm_RA = 0.0; // rads/year double pm_DEC = 0.0; // rads/year double parallax; // in arcsecs double radvel; // radial velocity (signed, km/s) }; struct engine_state_t { meteo_t meteo{.temperature = 0.0, .humidity = 0.5, .pressure = 1010.0}; double wavelength = DEFAULT_WAVELENGTH; // observed wavelength in mkm double lat = 0.0; // site latitude double 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 nameCCTE() const { return "ERFA-CCTE-ENGINE"; } // engine state related methods void setStateERFA(engine_state_t state) { std::lock_guard lock{*_stateMutex}; _currentState = std::move(state); } engine_state_t getStateERFA() const { std::lock_guard lock{*_stateMutex}; return _currentState; } void updateMeteoERFA(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 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> 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; } // ICRS to observed // azimuth must be counted from the South through the West error_t icrs2obs(mcc_angle_c auto ra_icrs, mcc_angle_c auto dec_icrs, mcc_angle_c auto* ra_obs, mcc_angle_c auto* dec_obs, mcc_angle_c auto* ha_obs, mcc_angle_c auto* az, mcc_angle_c auto* zd, obj_pars_t* obj_params = nullptr) { } error_t icrs2obs(MccSkyRADEC_ICRS const& radec_icrs, MccSkyRADEC_OBS* radec_obs, MccSkyAZZD* azzd, mcc_angle_c auto* ha_obs) { }; error_t obs2icrs(MccCoordPairKind obs_type, mcc_angle_c auto const& co_lon, mcc_angle_c auto const& co_lat, mcc_angle_c auto* ra_icrs, mcc_angle_c auto* dec_icrs) { } error_t obs2icrs(mcc_coord_pair_c auto const& xy_obs, MccSkyRADEC_ICRS* radec_icrs) {} error_t equationOrigins(const mcc_julday_c auto& mjd, mcc_angle_c auto* eo) { if (eo == nullptr) { return MccCCTE_ERFAErrorCode::ERROR_NULLPTR; } error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK; std::lock_guard lock{*_stateMutex}; using real_days_t = std::chrono::duration>; auto tai_utc = _currentState._leapSeconds[mjd]; if (tai_utc.has_value()) { double tt = mjd; tt += std::chrono::duration_cast(tai_utc.value()).count(); auto tt_tai = _currentState._bulletinA.TT_TAI(); tt += +std::chrono::duration_cast(tt_tai).count(); *eo = eraEo06a(ERFA_DJM0, tt); } else { ret = MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE; } return ret; } // refraction error_t refractionModel(refract_model_t* model) { if (model == nullptr) { return MccCCTE_ERFAErrorCode::ERROR_NULLPTR; } 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; } // Zobs must be observed zenithal distance (Zapp = Zobs + dZ -- corrected (in vacuo) zenithal distance) error_t refractionCorrection(const refract_model_t& rmodel, mcc_angle_c auto Zobs, mcc_angle_c auto* dZ, mcc_angle_c auto* Zapp = nullptr) { error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK; if (dZ == nullptr) { return MccCCTE_ERFAErrorCode::ERROR_NULLPTR; } if (Zobs >= std::numbers::pi / 2.0) { *dZ = 35.4 / 60.0 * std::numbers::pi / 180.0; // 35.4 arcminutes } else { auto tanZ = tan(Zobs); *dZ = rmodel.refa * tanZ + rmodel.refb * tanZ * tanZ * tanZ; } if (Zapp != nullptr) { *Zapp = Zobs + *dZ; } return ret; } // Zapp must be topocentric (in vacuo) zenithal distance (Zobs = Zapp - dZ -- observed, i.e. affected by refraction, // zenithal distance) error_t refractionReverseCorrection(const refract_model_t& rmodel, mcc_angle_c auto Zapp, mcc_angle_c auto* dZ, mcc_angle_c auto* Zobs = nullptr) { error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK; if (dZ == nullptr) { return MccCCTE_ERFAErrorCode::ERROR_NULLPTR; } if (Zapp >= std::numbers::pi / 2.0) { *dZ = 35.4 / 60.0 * std::numbers::pi / 180.0; // 35.4 arcminutes } else { auto tanZ = tan(Zapp); auto tanZ2 = tanZ * tanZ; auto b3 = 3.0 * rmodel.refb; // with Newton-Raphson correction *dZ = (rmodel.refa * tanZ + rmodel.refb * tanZ * tanZ2) / (1.0 + rmodel.refa + tanZ2 * (rmodel.refa + b3) + b3 * tanZ2 * tanZ2); } if (Zobs != nullptr) { *Zobs = Zapp - *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 _stateMutex; }; } // namespace mcc::ccte::erfa