diff --git a/mcc/mcc_ccte_erfa.h b/mcc/mcc_ccte_erfa.h index 9c48d81..4155092 100644 --- a/mcc/mcc_ccte_erfa.h +++ b/mcc/mcc_ccte_erfa.h @@ -299,12 +299,21 @@ public: using real_days_t = std::chrono::duration>; - double ut1 = jd.MJD(); - double tt = jd.MJD(); + using jd_t = decltype(jd); + + double mjd; + if constexpr (std::floating_point) { + mjd = jd - ERFA_DJM0; + } else { + mjd = jd.MJD(); + } + + double ut1 = mjd; + double tt = mjd; std::lock_guard lock{*_stateMutex}; - auto dut1 = _currentState._bulletinA.DUT1(jd.MJD()); + auto dut1 = _currentState._bulletinA.DUT1(mjd); if (dut1.has_value()) { ut1 += std::chrono::duration_cast(dut1.value()).count(); @@ -312,7 +321,7 @@ public: return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE; } - auto tai_utc = _currentState._leapSeconds[jd.MJD()]; + auto tai_utc = _currentState._leapSeconds[mjd]; if (tai_utc.has_value()) { tt += std::chrono::duration_cast(tai_utc.value()).count(); } else { @@ -471,7 +480,8 @@ public: double eo, lst, ha, dec, az, alt; auto lst_eo = [&]() { - ret = eqOrigins(from_pt.time_point, &eo); + // ret = eqOrigins(from_pt.time_point, &eo); + ret = equationOrigins(from_pt.time_point, &eo); if (!ret) { ret = timepointToAppSideral(from_pt.time_point, &lst, true); } @@ -566,7 +576,7 @@ public: return icrs2obs(from_pt, to_pt); } - // from apparent: copy corresponded coordinates and compute other ones + // from apparent: copy corresponded coordinates and compute other ones (ignore to_pt->pair_kind) if (from_pt.pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) { to_pt->RA_APP = from_pt.X; to_pt->DEC_APP = from_pt.Y; @@ -647,6 +657,58 @@ public: return ret; } + // equation of the origins + + error_t equationOrigins(mcc_julday_c auto const& jd, mcc_angle_c auto* eo) + { + if (eo == nullptr) { + return MccCCTE_ERFAErrorCode::ERROR_OK; + } + + auto ret = MccCCTE_ERFAErrorCode::ERROR_OK; + + double tt, mjd; + + if constexpr (std::floating_point) { + mjd = jd - ERFA_DJM0; + } else { + mjd = jd.MJD(); + } + + std::lock_guard lock{*_stateMutex}; + + using real_days_t = std::chrono::duration>; + + auto tai_utc = _currentState._leapSeconds[mjd]; + if (tai_utc.has_value()) { + tt = jd.MJD() + 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; + } + + error_t equationOrigins(mcc_time_point_c auto const& tp, mcc_angle_c auto* eo) + { + if (eo == nullptr) { + return MccCCTE_ERFAErrorCode::ERROR_OK; + } + + MccJulianDay jd; + auto ret = timepointToJulday(tp, &jd); + if (ret) { + return ret; + } + + return equationOrigins(std::move(jd), eo); + } + // refraction diff --git a/mcc/mcc_generics.h b/mcc/mcc_generics.h index 77cb801..7ee0ada 100644 --- a/mcc/mcc_generics.h +++ b/mcc/mcc_generics.h @@ -85,11 +85,19 @@ concept mcc_fp_type_like_c = template concept mcc_angle_c = mcc_fp_type_like_c && requires(T v, double vd) { + // mandatory arithmetic operations { v + v } -> std::same_as; { v - v } -> std::same_as; { v += v } -> std::same_as; { v -= v } -> std::same_as; + { vd + v } -> std::same_as; + { vd - v } -> std::same_as; + { v + vd } -> std::same_as; + { v - vd } -> std::same_as; + { v += vd } -> std::same_as; + { v -= vd } -> std::same_as; + { v * vd } -> std::same_as; { v / vd } -> std::same_as; }; @@ -101,7 +109,8 @@ concept mcc_angle_c = mcc_fp_type_like_c && requires(T v, double vd) { * USE OF STL std::chrono::time_point */ template -concept mcc_time_point_c = requires(T t) { [](std::chrono::time_point) {}(t); }; +concept mcc_time_point_c = traits::mcc_systime_c; +// concept mcc_time_point_c = requires(T t) { [](std::chrono::time_point) {}(t); }; /* JULIAN DAY CLASS CONCEPT */ @@ -123,6 +132,16 @@ concept mcc_error_c = std::convertible_to || requires(const T t) { }; +template +static constexpr ErrT mcc_deduce_error(const DErrT& err, const ErrT& default_err) +{ + if constexpr (std::same_as) { + return err; + } else { + return default_err; + } +} + /* ATMOSPHERIC REFRACTION MODEL CLASS CONCEPT */ template @@ -144,6 +163,27 @@ concept mcc_celestial_point_c = requires(T t) { }; +static constexpr void mcc_copy_celestial_point(mcc_celestial_point_c auto const& from_pt, + mcc_celestial_point_c auto* to_pt) +{ + if (to_pt == nullptr) { + return; + } + + using from_pt_t = std::remove_cvref_t; + using to_pt_t = std::remove_cvref_t; + + if constexpr (std::derived_from && std::copyable) { + *to_pt = from_pt; + return; + } + + to_pt->pair_kind = from_pt.pair_kind; + to_pt->time_point = std::chrono::time_point_casttime_point)::duration>(from_pt); + to_pt->X = (double)from_pt.X; + to_pt->Y = (double)from_pt.Y; +} + /* CELESTIAL POINT WITH APPARENT EQUATORIAL AND HORIZONTAL CLASS CONCEPT */ @@ -159,6 +199,34 @@ concept mcc_eqt_hrz_coord_c = mcc_celestial_point_c && requires(T t) { }; +static constexpr void mcc_copy_eqt_hrz_coord(mcc_eqt_hrz_coord_c auto const& from_pt, mcc_eqt_hrz_coord_c auto* to_pt) +{ + if (to_pt == nullptr) { + return; + } + + using from_pt_t = std::remove_cvref_t; + using to_pt_t = std::remove_cvref_t; + + if constexpr (std::derived_from && std::copyable) { + *to_pt = from_pt; + return; + } + + to_pt->pair_kind = from_pt.pair_kind; + to_pt->time_point = std::chrono::time_point_casttime_point)::duration>(from_pt); + to_pt->X = (double)from_pt.X; + to_pt->Y = (double)from_pt.Y; + + to_pt->RA_APP = (double)from_pt.RA_APP; + to_pt->DEC_APP = (double)from_pt.DEC_APP; + + to_pt->HA = (double)from_pt.HA; + to_pt->AZ = (double)from_pt.AZ; + to_pt->ZD = (double)from_pt.ZD; + to_pt->ALT = (double)from_pt.ALT; +} + /* CELESTIAL COORDINATES TRANSFORMATION ENGINE */ @@ -180,6 +248,12 @@ struct mcc_CCTE_interface_t { return std::forward(self).timepointToAppSideral(std::move(tp), st, islocal); } + template SelfT> + RetT juldayToAppSideral(this SelfT&& self, mcc_julday_c auto jd, mcc_angle_c auto* st, bool islocal = false) + { + return std::forward(self).timepointToAppSideral(std::move(jd), st, islocal); + } + // NOTE: ASSUMING THE AZINUTH IS COUNTED FROM THE SOUTH THROUGH THE WEST!!! template SelfT> RetT transformCoordinates(this SelfT&& self, mcc_celestial_point_c auto from_pt, mcc_celestial_point_c auto* to_pt) @@ -230,8 +304,43 @@ concept mcc_pointing_target_coord_c = mcc_eqt_hrz_coord_c && requires(T t) { requires mcc_angle_c; // ICRS declination }; + +static constexpr void mcc_copy_pointing_target_coord(mcc_pointing_target_coord_c auto const& from_pt, + mcc_pointing_target_coord_c auto* to_pt) +{ + if (to_pt == nullptr) { + return; + } + + using from_pt_t = std::remove_cvref_t; + using to_pt_t = std::remove_cvref_t; + + if constexpr (std::derived_from && std::copyable) { + *to_pt = from_pt; + return; + } + + to_pt->pair_kind = from_pt.pair_kind; + to_pt->time_point = std::chrono::time_point_casttime_point)::duration>(from_pt); + to_pt->X = (double)from_pt.X; + to_pt->Y = (double)from_pt.Y; + + to_pt->RA_ICRS = (double)from_pt.RA_ICRS; + to_pt->DEC_ICRS = (double)from_pt.DEC_ICRS; + + to_pt->RA_APP = (double)from_pt.RA_APP; + to_pt->DEC_APP = (double)from_pt.DEC_APP; + + to_pt->HA = (double)from_pt.HA; + to_pt->AZ = (double)from_pt.AZ; + to_pt->ZD = (double)from_pt.ZD; + to_pt->ALT = (double)from_pt.ALT; +} + + + template -concept mcc_telemetry_data_c = mcc_eqt_hrz_coord_c && requires(T t) { +concept mcc_telemetry_data_c = mcc_eqt_hrz_coord_c && std::default_initializable && requires(T t) { // target target coordinates requires mcc_pointing_target_coord_c; @@ -249,6 +358,46 @@ concept mcc_telemetry_data_c = mcc_eqt_hrz_coord_c && requires(T t) { }; +static constexpr void mcc_copy_telemetry_data(mcc_telemetry_data_c auto const& from_pt, + mcc_telemetry_data_c auto* to_pt) +{ + if (to_pt == nullptr) { + return; + } + + using from_pt_t = std::remove_cvref_t; + using to_pt_t = std::remove_cvref_t; + + if constexpr (std::derived_from && std::copyable) { + *to_pt = from_pt; + return; + } + + to_pt->pair_kind = from_pt.pair_kind; + to_pt->time_point = std::chrono::time_point_casttime_point)::duration>(from_pt); + + to_pt->X = (double)from_pt.X; + to_pt->Y = (double)from_pt.Y; + to_pt->speedX = (double)from_pt.speedX; + to_pt->speedY = (double)from_pt.speedY; + + to_pt->RA_APP = (double)from_pt.RA_APP; + to_pt->DEC_APP = (double)from_pt.DEC_APP; + + to_pt->HA = (double)from_pt.HA; + to_pt->AZ = (double)from_pt.AZ; + to_pt->ZD = (double)from_pt.ZD; + to_pt->ALT = (double)from_pt.ALT; + + to_pt->pcmX = (double)from_pt.pcmX; + to_pt->pcmY = (double)from_pt.pcmY; + + to_pt->refCorr = (double)from_pt.refCorr; + + mcc_copy_pointing_target_coord(from_pt.target, &to_pt->target); +} + + /* MOUNT TELEMETRY MANAGER CLASS CONCEPT */ diff --git a/mcc/mcc_telemetry.h b/mcc/mcc_telemetry.h index 850ed2c..ce7377c 100644 --- a/mcc/mcc_telemetry.h +++ b/mcc/mcc_telemetry.h @@ -14,7 +14,13 @@ namespace mcc static constexpr double mcc_sideral_to_UT1_ratio = 1.002737909350795; // sideral/UT1 -enum MccTelemetryErrorCode : int { ERROR_OK, ERROR_COORD_TRANSFROM, ERROR_PCM_COMP }; +enum MccTelemetryErrorCode : int { + ERROR_OK, + ERROR_NULLPTR, + ERROR_COORD_TRANSFORM, + ERROR_PCM_COMP, + ERROR_HARDWARE_GETPOS +}; } // namespace mcc @@ -50,10 +56,14 @@ struct MccTelemetryCategory : public std::error_category { switch (err) { case MccTelemetryErrorCode::ERROR_OK: return "OK"; - case MccTelemetryErrorCode::ERROR_COORD_TRANSFROM: + case MccTelemetryErrorCode::ERROR_NULLPTR: + return "nullptr input argument"; + case MccTelemetryErrorCode::ERROR_COORD_TRANSFORM: return "coordinate transformation error"; case MccTelemetryErrorCode::ERROR_PCM_COMP: return "PCM computation error"; + case MccTelemetryErrorCode::ERROR_HARDWARE_GETPOS: + return "cannot get hardware position"; default: return "UNKNOWN"; } @@ -87,101 +97,212 @@ public: typedef std::error_code error_t; - MccTelemetry(mcc_ccte_c auto* ccte, mcc_PCM_c auto* pcm, mcc_hardware_c auto* hardware) + MccTelemetry(mcc_ccte_c auto* ccte, mcc_PCM_c auto* pcm, mcc_hardware_c auto* hardware) : _data() { + _data.target.pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_ICRS; + using ccte_t = std::remove_cvref_t; using pcm_t = std::remove_cvref_t; using hardware_t = std::remove_cvref_t; + + _updateTargetFunc = [ccte, this]() { + // + // ICRS coordinates of the taget must be already set + // + _data.target.pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_ICRS; + _data.target.X = _data.target.RA_ICRS; + _data.target.Y = _data.target.DEC_ICRS; + + // update apparent cordinates + auto ret = ccte->transformCoordinates(_data.target, &_data.target); + + return mcc_deduce_error(ret, MccTelemetryErrorCode::ERROR_COORD_TRANSFORM); + }; + _updateFunc = [ccte, pcm, hardware, this](MccTelemetryData* data) { + // first, update mount quantities typename hardware_t::axes_pos_t hw_pos; auto hw_err = hardware->getPos(&hw_pos); + if (hw_err) { + return mcc_deduce_error(hw_err, MccTelemetryErrorCode::ERROR_HARDWARE_GETPOS); + } - if (!hw_err) { - data->time_point = - std::chrono::time_point_casttime_point)::duration>(hw_pos.time_point); + double eo; - auto ccte_err = ccte->timepointToJulday(data->time_point, &data->JD); + data->time_point = + std::chrono::time_point_casttime_point)::duration>(hw_pos.time_point); + + auto ccte_err = ccte->timepointToJulday(data->time_point, &data->JD); + if (!ccte_err) { + ccte_err = ccte->juldayToAppSideral(data->JD, &data->LST, true); if (!ccte_err) { - ccte_err = ccte->timepointToAppSideral(data->time_point, &data->LST, true); - if (!ccte_err) { - data->speedX = (double)hw_pos.speedX; - data->speedY = (double)hw_pos.speedY; - - struct { - double dx, dy; - } pcm_res; - - auto pcm_err = pcm->computePCM(data, &pcm_res); - if (pcm_err) { - if constexpr (std::same_as) { - return pcm_err; - } else { - return MccTelemetryErrorCode::ERROR_PCM_COMP; - } - } - - MccCelestialPoint pt{.pair_kind = MccCoordPairKind::COORDS_KIND_AZALT, - .time_point = data->time_point}; - - if constexpr (mccIsEquatorialMount(pcm_t::mountType)) { - data->HA = (double)hw_pos.X + pcm_res.dx; - data->DEC_APP = (double)hw_pos.Y + pcm_res.dy; - - data->X = data->HA; - data->Y = data->DEC_APP; - - data->pair_kind = MccCoordPairKind::COORDS_KIND_HADEC_APP; - - ccte_err = ccte->transformCoordinates(*data, &pt); - if (!ccte_err) { - data->AZ = pt.X; - data->ALT = pt.Y; - data->ZD = std::numbers::pi / 2.0 - data->ALT; - } - } else if constexpr (mccIsAltAzMount(pcm_t::mountType)) { - data->AZ = (double)hw_pos.X + pcm_res.dx; - data->ALT = (double)hw_pos.Y + pcm_res.dy; - data->ZD = std::numbers::pi / 2.0 - data->ALT; - - data->X = data->AZ; - data->Y = data->ALT; - - data->pair_kind = MccCoordPairKind::COORDS_KIND_AZALT; - - pt.pair_kind = MccCoordPairKind::COORDS_KIND_HADEC_APP; - ccte_err = ccte->transformCoordinates(*data, &pt); - if (!ccte) { - data->HA = pt.X; - data->DEC_APP = pt.Y; - } - - } else { - static_assert(false, "UNKNOWN MOUNT TYPE!"); - } - } + ccte_err = ccte->equationOrigins(data->JD, &eo); } } + + if (ccte_err) { + return mcc_deduce_error(ccte_err, MccTelemetryErrorCode::ERROR_COORD_TRANSFORM); + } + + data->speedX = (double)hw_pos.speedX; + data->speedY = (double)hw_pos.speedY; + + struct { + double dx, dy; + } pcm_res; + + auto pcm_err = pcm->computePCM(data, &pcm_res); + if (pcm_err) { + return mcc_deduce_error(pcm_err, MccTelemetryErrorCode::ERROR_PCM_COMP); + } + + data->pcmX = pcm_res.dx; + data->pcmY = pcm_res.dy; + + MccCelestialPoint pt{.pair_kind = MccCoordPairKind::COORDS_KIND_AZALT, .time_point = data->time_point}; + + if constexpr (mccIsEquatorialMount(pcm_t::mountType)) { + data->HA = (double)hw_pos.X + pcm_res.dx; + data->DEC_APP = (double)hw_pos.Y + pcm_res.dy; + data->RA_APP = (double)data->LST - (double)data->HA + eo; + + data->X = data->HA; + data->Y = data->DEC_APP; + + data->pair_kind = MccCoordPairKind::COORDS_KIND_HADEC_APP; + + ccte_err = ccte->transformCoordinates(*data, &pt); + if (!ccte_err) { + data->AZ = pt.X; + data->ALT = pt.Y; + data->ZD = std::numbers::pi / 2.0 - data->ALT; + } + } else if constexpr (mccIsAltAzMount(pcm_t::mountType)) { + data->AZ = (double)hw_pos.X + pcm_res.dx; + data->ALT = (double)hw_pos.Y + pcm_res.dy; + data->ZD = std::numbers::pi / 2.0 - data->ALT; + + data->X = data->AZ; + data->Y = data->ALT; + + data->pair_kind = MccCoordPairKind::COORDS_KIND_AZALT; + + pt.pair_kind = MccCoordPairKind::COORDS_KIND_HADEC_APP; + ccte_err = ccte->transformCoordinates(*data, &pt); + if (!ccte) { + data->HA = pt.X; + data->DEC_APP = pt.Y; + data->RA_APP = (double)data->LST - (double)data->HA + eo; + } + + } else { + static_assert(false, "UNKNOWN MOUNT TYPE!"); + } + + if (!ccte_err) { + data->pair_kind = MccCoordPairKind::COORDS_KIND_AZZD; + data->X = data->AZ; + data->Y = data->ZD; + + ccte_err = ccte->refractionCorrection(data, &data->refCorr); + if (!ccte_err) { + // hardware encoders coordinates + data->X = (double)hw_pos.X; + data->Y = (double)hw_pos.Y; + + // update target (assuming target ICRS coordinates are already set) + + data->target.time_point = + std::chrono::time_point_casttarget.time_point)::duration>( + data->time_point); + + data->target.pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_ICRS; + data->target.X = data->target.RA_ICRS; + data->target.Y = data->target.DEC_ICRS; + + ccte_err = ccte->transformCoordinates(data->target, &data->target); + } + } + + if (ccte_err) { + return mcc_deduce_error(ccte_err, MccTelemetryErrorCode::ERROR_COORD_TRANSFORM); + } + + if constexpr (mccIsEquatorialMount(pcm_t::mountType)) { + data->pair_kind = MccCoordPairKind::COORDS_KIND_HADEC_APP; + } else if constexpr (mccIsAltAzMount(pcm_t::mountType)) { + data->pair_kind = MccCoordPairKind::COORDS_KIND_AZALT; + } else { + static_assert(false, "UNKNOWN MOUNT TYPE!"); + } + + + return MccTelemetryErrorCode::ERROR_OK; + }; + + + _setTargetFunc = [ccte, this](MccCelestialPoint const& pt) { + // in the case of apparent input coordinates + // one must ensure the same time points + + _data.target.pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_ICRS; + _data.target.time_point = + std::chrono::time_point_cast(pt.time_point); + + auto ret = ccte->transformCoordinates(pt, &_data.target); + + if (!ret) { + if (pt.pair_kind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) { + _data.target.RA_ICRS = _data.target.X; + _data.target.DEC_ICRS = _data.target.Y; + + // update apparent coordinates + ret = _updateTargetFunc(); + } else { // apparent coordinates were computed above + // compute ICRS coordinates + MccCelestialPoint cpt{.pair_kind = MccCoordPairKind::COORDS_KIND_RADEC_ICRS}; + ret = ccte->transformCoordinates(pt, &cpt); + + _data.target.RA_ICRS = cpt.X; + _data.target.DEC_ICRS = cpt.Y; + } + } + + return mcc_deduce_error(ret, MccTelemetryErrorCode::ERROR_COORD_TRANSFORM); }; } - error_t telemetryData(mcc_telemetry_data_c auto* data) + virtual ~MccTelemetry() = default; + + error_t telemetryData(mcc_telemetry_data_c auto* tdata) { - error_t ret = MccTelemetryErrorCode::ERROR_OK; - - // first, update mount coordinates + if (tdata == nullptr) { + return MccTelemetryErrorCode::ERROR_NULLPTR; + } + error_t ret = _updateFunc(&_data); + if (!ret) { + mcc_copy_telemetry_data(_data, tdata); + } return ret; } - error_t setPointingTarget(mcc_celestial_point_c auto pt) {} + error_t setPointingTarget(mcc_celestial_point_c auto pt) + { + // return + } protected: - std::function _updateFunc; + MccTelemetryData _data; + + std ::function _updateTargetFunc{}; + std::function _updateFunc{}; + std::function _setTargetFunc{}; }; } // namespace mcc diff --git a/mcc/tests/ccte_test.cpp b/mcc/tests/ccte_test.cpp index 9808773..b126af5 100644 --- a/mcc/tests/ccte_test.cpp +++ b/mcc/tests/ccte_test.cpp @@ -47,9 +47,11 @@ int main() std::cout << "JD = " << std::setprecision(12) << (double)jd << "\n"; - mcc::MccAngle lst; + mcc::MccAngle lst, eo; ret = erfa.timepointToAppSideral(now, &lst, true); std::cout << "LST = " << lst.sexagesimal(true) << "\n"; + ret = erfa.equationOrigins(now, &eo); + std::cout << "EO = " << eo.sexagesimal(true) << "\n"; std::cout << "\n\n";