diff --git a/cxx/CMakeLists.txt b/cxx/CMakeLists.txt index b969070..b442114 100644 --- a/cxx/CMakeLists.txt +++ b/cxx/CMakeLists.txt @@ -124,7 +124,7 @@ add_library(${CNTR_PROTO_LIB} STATIC ${CNTR_PROTO_LIB_SRC}) set(MCC_LIBRARY_SRC mcc_mount.h mcc_mount_coord.h mcc_mount_events_states.h mcc_finite_state_machine.h - mcc_mount_pec.h mcc_mount_pz.h mcc_traits.h mcc_mount_telemetry.h mcc_mount_config.h) + mcc_mount_pec.h mcc_mount_pz.h mcc_traits.h mcc_mount_telemetry.h mcc_mount_config.h mcc_mount_astro_erfa.h) set(MCC_LIBRARY mcc) add_library(${MCC_LIBRARY} INTERFACE ${MCC_LIBRARY_SRC}) target_compile_features(${MCC_LIBRARY} INTERFACE cxx_std_23) diff --git a/cxx/mcc_astrom_iers.h b/cxx/mcc_astrom_iers.h index 93fb29e..187fa8b 100644 --- a/cxx/mcc_astrom_iers.h +++ b/cxx/mcc_astrom_iers.h @@ -42,6 +42,11 @@ public: return _expireDate; } + auto expireMJD() const + { + return _expireMJD; + } + // load from stream bool load(std::derived_from> auto& stream, char comment_sym = '#') { @@ -118,7 +123,7 @@ public: int64_t mjd_int = (1461LL * (iypmy + 4800LL)) / 4LL + (367LL * (im - 2LL - 12LL * my)) / 12LL - (3LL * ((iypmy + 4900LL) / 100LL)) / 4LL + id - 2432076LL; - _expireJulDate = static_cast(mjd_int); + _expireMJD = static_cast(mjd_int); _db = std::move(db); @@ -169,7 +174,7 @@ public: { double e_mjd; - if (mjd > _expireJulDate) { // ???????!!!!!!!!!!! + if (mjd > _expireMJD) { // ???????!!!!!!!!!!! return std::nullopt; // return _db.back().tai_utc; } @@ -201,7 +206,7 @@ private: inline static const std::regex data_rx{"^ *[0-9]{5,}(\\.?[0-9]+) +[0-9]{1,2} +[0-9]{1,2} +[0-9]{4} +[0-9]{1,} *$"}; time_point_t _expireDate{}; - double _expireJulDate{}; + double _expireMJD{}; struct leapsecond_db_elem_t { double mjd; diff --git a/cxx/mcc_mount_astro_erfa.h b/cxx/mcc_mount_astro_erfa.h index ef7084f..f920d27 100644 --- a/cxx/mcc_mount_astro_erfa.h +++ b/cxx/mcc_mount_astro_erfa.h @@ -6,12 +6,14 @@ /* ASTROMETRY ENGINE BASED ON ERFA-LIBRARY */ #include -// #include +#include #include "mcc_astrom_iers.h" #include "mcc_mount_coord.h" +#include "mcc_mount_astrom.h" + namespace mcc::astrom::erfa { @@ -20,8 +22,28 @@ namespace mcc::astrom::erfa class MccMountAstromEngineERFA { +protected: + static constexpr std::array engineErrorString = {"OK", + "invalid argument", + "invalid year", + "invalid month", + "time point is out of range", + "time point is out of range", + "dubious year", + "unacceptable year"}; + public: - typedef int engine_err_t; + enum engine_err_t : size_t { + ERROR_OK = 0, + ERROR_INVALID_INPUT_ARG, + ERROR_JULDATE_INVALID_YEAR, + ERROR_JULDATE_INVALID_MONTH, + ERROR_BULLETINA_OUT_OF_RANGE, + ERROR_LEAPSECONDS_OUT_OF_RANGE, + ERROR_DUBIOUS_YEAR, + ERROR_UNACCEPTABLE_DATE + }; + // meteo parameters (to compute refraction) struct meteo_t { @@ -41,6 +63,7 @@ public: MccAngleLAT lat = "00:00:00"_dms; // site latitude MccAngleLON 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{}; @@ -54,6 +77,9 @@ public: }; typedef MccAngle coord_t; + typedef MccAngle prop_motion_t; + typedef double parallax_t; + typedef MccAngle gst_t; typedef MccAngle pa_t; typedef double eo_t; @@ -62,8 +88,220 @@ public: double refa, refb; }; + + MccMountAstromEngineERFA() = default; + + MccMountAstromEngineERFA(engine_state_t state) : _currentState(std::move(state)) {} + + virtual ~MccMountAstromEngineERFA() = default; + + void updateState(engine_state_t state) + { + std::lock_guard lock{_stateMutex}; + + _currentState = std::move(state); + } + + engine_state_t getState() const + { + std::scoped_lock lock{_stateMutex}; + + return _currentState; + } + + + std::string errorString(engine_err_t err) const + { + return engineErrorString[err]; + } + + /* time-related methods */ + + // templated generic version + template + engine_err_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 ERROR_JULDATE_INVALID_YEAR; + } + if (!ymd.month().ok()) { + return 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 ERROR_OK; + } + + engine_err_t greg2jul(time_point_t time_point, juldate_t& juldate) + { + return greg2jul(time_point, juldate); + } + + + engine_err_t apparentSiderTime(juldate_t juldate, gst_t& gst, bool islocal = false) + { + std::lock_guard lock{_stateMutex}; + + using real_days_t = std::chrono::duration>; + + double ut1, tt = juldate.mjd; + + auto dut1 = _currentState._bulletinA.DUT1(juldate.mjd); + + if (dut1.has_value()) { + ut1 += std::chrono::duration_cast(dut1.value()).count(); + } else { // out of range + return 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); + if (islocal) { + gst += _currentState.lon; + } + + return ERROR_OK; + } + + + /* atmospheric refraction-related methods */ + + engine_err_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 ERROR_OK; + } + + + /* coordinates conversional methods */ + + engine_err_t icrs2obs(coord_t ra, + coord_t dec, + prop_motion_t pm_ra, + prop_motion_t pm_dec, + parallax_t parallax, + 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 ERROR_BULLETINA_OUT_OF_RANGE; + } + + auto pol_pos = _currentState._bulletinA.polarCoords(juldate.mjd); + if (!pol_pos.has_value()) { + return ERROR_BULLETINA_OUT_OF_RANGE; + } + + double oaz, ozd, oha, odec, ora; + + int ret = eraAtco13(ra, dec, pm_ra, pm_dec, parallax, 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 ERROR_DUBIOUS_YEAR; + } else if (ret == -1) { + return ERROR_UNACCEPTABLE_DATE; + } + + ra_app = ora; + dec_app = odec; + az = oaz; + alt = std::numbers::pi / 2.0 - ozd; + ha = oha; + + return ERROR_OK; + } + + + engine_err_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 ERROR_OK; + } + + engine_err_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 ERROR_OK; + } + + + engine_err_t hadec2pa(coord_t ha, coord_t dec, pa_t& pa) + { + std::lock_guard lock{_stateMutex}; + + pa = eraHd2pa(ha, dec, _currentState.lat); + + return ERROR_OK; + } + protected: - engine_state_t _currentState; + engine_state_t _currentState{}; + + std::mutex _stateMutex; }; } // namespace mcc::astrom::erfa + + +static_assert(mcc::traits::mcc_astrom_engine_c, ""); diff --git a/cxx/mcc_mount_astrom.h b/cxx/mcc_mount_astrom.h index 388b793..db14a50 100644 --- a/cxx/mcc_mount_astrom.h +++ b/cxx/mcc_mount_astrom.h @@ -7,6 +7,7 @@ #include +#include // #include "mcc_traits.h" @@ -14,11 +15,14 @@ namespace mcc::traits { template -concept mcc_astrom_engine_c = requires(T t) { +concept mcc_astrom_engine_c = requires(T t, const T t_const) { typename T::engine_err_t; typename T::engine_state_t; + requires std::movable; typename T::coord_t; + typename T::prop_motion_t; + typename T::parallax_t; typename T::time_point_t; typename T::juldate_t; typename T::gst_t; @@ -30,15 +34,20 @@ concept mcc_astrom_engine_c = requires(T t) { { t.updateState(std::declval()) }; + { t_const.getState() } -> std::same_as; + + { t_const.errorString(std::declval()) } -> std::convertible_to; + /* coordinates conversional methods */ // ICRS RA and DEC to observed place: icrs2obs(ra, dec, jd, ra_app, dec_app, ha, az, alt, eo) { t.icrs2obs(std::declval(), std::declval(), - std::declval(), std::declval(), + std::declval(), std::declval(), + std::declval(), std::declval(), std::declval(), std::declval(), std::declval(), std::declval(), - std::declval()) + std::declval(), std::declval()) } -> std::same_as; // compute hour angle and declination from azimuth and altitude: hadec2azalt(ha, dec, az, alt) @@ -69,9 +78,10 @@ concept mcc_astrom_engine_c = requires(T t) { // requires mcc_systime_c>; // requires mcc_output_arg_c, typename T::juldate_t>; - // apparent sideral time: apparentSiderTime(jd, gst) + // apparent sideral time: apparentSiderTime(jd, gst, islocal) { - t.apparentSiderTime(std::declval(), std::declval()) + t.apparentSiderTime(std::declval(), std::declval(), + std::declval()) } -> std::same_as;