401 lines
11 KiB
C++
401 lines
11 KiB
C++
#pragma once
|
|
|
|
#include <erfa.h>
|
|
#include <erfam.h>
|
|
#include <mutex>
|
|
|
|
#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<mcc::ccte::erfa::MccCCTE_ERFAErrorCode> : 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<MccCCTE_ERFAErrorCode>(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<int>(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<std::basic_istream<char>> 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<std::basic_istream<char>> 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;
|
|
}
|
|
}
|
|
|
|
// 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<std::mutex> _stateMutex;
|
|
};
|
|
|
|
|
|
} // namespace mcc::ccte::erfa
|