576 lines
18 KiB
C++
576 lines
18 KiB
C++
#pragma once
|
|
|
|
/* MOUNT CONTROL COMPONENTS LIBRARY */
|
|
|
|
|
|
/* ASTROMETRY ENGINE BASED ON ERFA-LIBRARY (THREAD-SAFE FOR ENGINE STATE MANIPULATIONS) */
|
|
|
|
#include <chrono>
|
|
#include <mutex>
|
|
|
|
|
|
#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<mcc::astrom::erfa::MccMountAstromEngineERFAErrorCode> : public true_type
|
|
{
|
|
};
|
|
|
|
} // namespace std
|
|
|
|
|
|
|
|
namespace mcc::astrom::erfa
|
|
{
|
|
|
|
#include <erfa.h>
|
|
#include <erfam.h>
|
|
|
|
|
|
/* 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<MccMountAstromEngineERFAErrorCode>(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<int>(ec), MccMountAstromEngineERFACategory::get());
|
|
}
|
|
|
|
|
|
|
|
/* A concept for ERFA-library compatible type to represent anglular quantities */
|
|
template <typename T>
|
|
concept mcc_erfa_angle_t = std::constructible_from<T, double> && std::convertible_to<T, double>;
|
|
|
|
template <mcc_erfa_angle_t AngleT = MccAngle>
|
|
class MccMountAstromEngineERFA
|
|
{
|
|
public:
|
|
static constexpr double DEFAULT_WAVELENGTH = 0.55; // default observed wavelength in mkm
|
|
|
|
typedef std::error_code error_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
|
|
|
|
AngleT lat = "00:00:00"_dms; // site latitude
|
|
AngleT 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 sideral_time_t;
|
|
// typedef MccAngle pa_t;
|
|
// typedef MccAngle eo_t;
|
|
|
|
|
|
/* use of the same type fro representation of celestial and geodetic coordinates, celestial angles (e.g. P.A.),
|
|
* sideral time */
|
|
typedef AngleT coord_t;
|
|
|
|
typedef AngleT sideral_time_t;
|
|
typedef AngleT pa_t;
|
|
typedef AngleT eo_t;
|
|
|
|
struct refract_result_t {
|
|
double refa, refb;
|
|
};
|
|
|
|
|
|
MccMountAstromEngineERFA() = default;
|
|
|
|
MccMountAstromEngineERFA(engine_state_t state) : _currentState(std::move(state)) {}
|
|
|
|
MccMountAstromEngineERFA(MccMountAstromEngineERFA&&) = default;
|
|
MccMountAstromEngineERFA& operator=(MccMountAstromEngineERFA&&) = default;
|
|
|
|
MccMountAstromEngineERFA(const MccMountAstromEngineERFA&) = delete;
|
|
MccMountAstromEngineERFA& operator=(const MccMountAstromEngineERFA&) = delete;
|
|
|
|
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<std::basic_istream<char>> 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<std::basic_istream<char>> 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 <mcc::traits::mcc_systime_c TpT>
|
|
error_t greg2jul(TpT time_point, juldate_t& juldate)
|
|
{
|
|
using namespace std::literals::chrono_literals;
|
|
|
|
auto dd = std::chrono::floor<std::chrono::days>(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<double, std::ratio<86400>>; // fraction of day
|
|
juldate.mjd = static_cast<double>(mjd_int) + std::chrono::duration_cast<fd_t>(time_point - dd).count();
|
|
|
|
return MccMountAstromEngineERFAErrorCode::ERROR_OK;
|
|
}
|
|
|
|
error_t greg2jul(time_point_t time_point, juldate_t& juldate)
|
|
{
|
|
return greg2jul<time_point_t>(time_point, juldate);
|
|
}
|
|
|
|
|
|
error_t terrestrialTime(juldate_t juldate, juldate_t& tt)
|
|
{
|
|
std::lock_guard lock{_stateMutex};
|
|
|
|
using real_days_t = std::chrono::duration<double, std::ratio<86400>>;
|
|
|
|
auto tai_utc = _currentState._leapSeconds[juldate.mjd];
|
|
if (tai_utc.has_value()) {
|
|
tt.mjd += std::chrono::duration_cast<real_days_t>(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<real_days_t>(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, std::ratio<86400>>;
|
|
|
|
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<real_days_t>(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<real_days_t>(tai_utc.value()).count();
|
|
// } else {
|
|
// return ERROR_LEAPSECONDS_OUT_OF_RANGE;
|
|
// }
|
|
|
|
|
|
// auto tt_tai = _currentState._bulletinA.TT_TAI();
|
|
// tt += std::chrono::duration_cast<real_days_t>(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;
|
|
}
|
|
|
|
|
|
/* 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{};
|
|
|
|
mutable std::mutex _stateMutex;
|
|
};
|
|
|
|
} // namespace mcc::astrom::erfa
|
|
|
|
|
|
static_assert(mcc::traits::mcc_astrom_engine_c<mcc::astrom::erfa::MccMountAstromEngineERFA<>>, "");
|