#pragma once /* MOUNT CONTROL COMPONENTS LIBRARY */ /* ASTROMETRY ENGINE BASED ON ERFA-LIBRARY */ #include #include #include "mcc_astrom_iers.h" #include "mcc_mount_coord.h" #include "mcc_mount_astrom.h" namespace mcc::astrom::erfa { #include #include 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: 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 { 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 = 0.55; // observed wavelength in mkm MccAngle lat = "00:00:00"_dms; // site latitude MccAngle 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 MccAngle coord_t; typedef MccAngle prop_motion_t; typedef double parallax_t; typedef MccAngle gst_t; typedef MccAngle pa_t; typedef double eo_t; struct refract_result_t { double refa, refb; }; MccMountAstromEngineERFA() = default; MccMountAstromEngineERFA(engine_state_t state) : _currentState(std::move(state)) {} 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; } 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 = juldate.mjd, 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; } const auto arcsec2rad = std::numbers::pi / 180 / 3600; pol_pos->x *= arcsec2rad; pol_pos->y *= arcsec2rad; 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{}; mutable std::mutex _stateMutex; }; } // namespace mcc::astrom::erfa static_assert(mcc::traits::mcc_astrom_engine_c, "");