837 lines
26 KiB
C++
837 lines
26 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;
|
|
}
|
|
|
|
|
|
// apparent sideral time (Greenwitch or local)
|
|
error_t apparentSideralTime(mcc_coord_epoch_c auto const& epoch, mcc_angle_c auto* st, bool islocal = false)
|
|
{
|
|
error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
if (st == nullptr) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
|
|
}
|
|
|
|
using real_days_t = std::chrono::duration<double, std::ratio<86400>>;
|
|
|
|
double ut1 = epoch.MJD();
|
|
double tt = epoch.MJD();
|
|
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
auto dut1 = _currentState._bulletinA.DUT1(epoch.MJD());
|
|
|
|
if (dut1.has_value()) {
|
|
ut1 += std::chrono::duration_cast<real_days_t>(dut1.value()).count();
|
|
} else { // out of range
|
|
return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
|
|
}
|
|
|
|
auto tai_utc = _currentState._leapSeconds[epoch.MJD()];
|
|
if (tai_utc.has_value()) {
|
|
tt += std::chrono::duration_cast<real_days_t>(tai_utc.value()).count();
|
|
} else {
|
|
return MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE;
|
|
}
|
|
|
|
|
|
auto tt_tai = _currentState._bulletinA.TT_TAI();
|
|
tt += std::chrono::duration_cast<real_days_t>(tt_tai).count();
|
|
|
|
*st = eraGst06a(ERFA_DJM0, ut1, ERFA_DJM0, tt);
|
|
|
|
if (islocal) {
|
|
*st = eraAnp(*st + _currentState.lon);
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
// ICRS to observed
|
|
// returned azimuth is counted from the South through the West
|
|
error_t icrsToObs(mcc_angle_c auto const& ra_icrs,
|
|
mcc_angle_c auto const& dec_icrs,
|
|
mcc_coord_epoch_c auto const& epoch,
|
|
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)
|
|
{
|
|
return icrsTo(true, ra_icrs, dec_icrs, epoch, ra_obs, dec_obs, ha_obs, az, zd, obj_params);
|
|
}
|
|
|
|
error_t icrsToObs(MccSkyRADEC_ICRS const& radec_icrs,
|
|
MccSkyRADEC_OBS* radec_obs,
|
|
MccSkyAZZD* azzd,
|
|
mcc_angle_c auto* ha_obs,
|
|
obj_pars_t* obj_params = nullptr)
|
|
{
|
|
double ra_obs, dec_obs, az, zd, ha;
|
|
|
|
auto err =
|
|
icrsToObs(radec_icrs.x(), radec_icrs.y(), radec_icrs.epoch(), &ra_obs, &dec_obs, &ha, &az, &zd, obj_params);
|
|
|
|
if (!err) {
|
|
if (radec_obs) {
|
|
radec_obs->setX(ra_obs);
|
|
radec_obs->setY(dec_obs);
|
|
}
|
|
|
|
if (azzd) {
|
|
azzd->setEpoch(radec_obs->epoch());
|
|
azzd->setX(az);
|
|
azzd->setY(zd);
|
|
}
|
|
|
|
if (ha_obs) {
|
|
*ha_obs = ha;
|
|
}
|
|
}
|
|
|
|
return err;
|
|
};
|
|
|
|
|
|
// ICRS to apparent (in vacuo)
|
|
// returned azimuth is counted from the South through the West
|
|
error_t icrsToApp(mcc_angle_c auto const& ra_icrs,
|
|
mcc_angle_c auto const& dec_icrs,
|
|
mcc_coord_epoch_c auto const& epoch,
|
|
mcc_angle_c auto* ra_app,
|
|
mcc_angle_c auto* dec_app,
|
|
mcc_angle_c auto* ha_app,
|
|
mcc_angle_c auto* az,
|
|
mcc_angle_c auto* zd, // should be interpretated as zenithal distance corrected for refraction
|
|
obj_pars_t* obj_params = nullptr)
|
|
{
|
|
return icrsTo(false, ra_icrs, dec_icrs, epoch, ra_app, dec_app, ha_app, az, zd, obj_params);
|
|
}
|
|
|
|
|
|
error_t icrsToApp(MccSkyRADEC_ICRS const& radec_icrs,
|
|
MccSkyRADEC_OBS* radec_app,
|
|
MccSkyAZZD* azzd,
|
|
mcc_angle_c auto* ha_app,
|
|
obj_pars_t* obj_params = nullptr)
|
|
{
|
|
double ra_app, dec_app, az, zd, ha;
|
|
|
|
auto err =
|
|
icrsToApp(radec_icrs.x(), radec_icrs.y(), radec_icrs.epoch(), &ra_app, &dec_app, &ha, &az, &zd, obj_params);
|
|
|
|
if (!err) {
|
|
if (radec_app) {
|
|
radec_app->setX(ra_app);
|
|
radec_app->setY(dec_app);
|
|
}
|
|
|
|
if (azzd) {
|
|
azzd->setEpoch(radec_app->epoch());
|
|
azzd->setX(az);
|
|
azzd->setY(zd);
|
|
}
|
|
|
|
if (ha_app) {
|
|
*ha_app = ha;
|
|
}
|
|
}
|
|
|
|
return err;
|
|
}
|
|
|
|
|
|
error_t obsToICRS(MccCoordPairKind obs_type,
|
|
mcc_coord_epoch_c auto const& epoch,
|
|
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)
|
|
{
|
|
return toICRS(true, obs_type, epoch, co_lon, co_lat, ra_icrs, dec_icrs);
|
|
}
|
|
|
|
|
|
error_t obsToICRS(mcc_coord_pair_c auto const& xy_obs, MccSkyRADEC_ICRS* radec_icrs)
|
|
{
|
|
double ra, dec;
|
|
|
|
auto err = obsToICRS(xy_obs.pair_kind, xy_obs.epoch(), xy_obs.x(), xy_obs.y(), &ra, &dec);
|
|
if (err) {
|
|
return err;
|
|
}
|
|
|
|
if (radec_icrs) {
|
|
radec_icrs->setX(ra);
|
|
radec_icrs->setY(dec);
|
|
}
|
|
|
|
return err;
|
|
}
|
|
|
|
|
|
error_t appToICRS(MccCoordPairKind app_type,
|
|
mcc_coord_epoch_c auto const& epoch,
|
|
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)
|
|
{
|
|
return toICRS(false, app_type, epoch, co_lon, co_lat, ra_icrs, dec_icrs);
|
|
}
|
|
|
|
|
|
error_t appToICRS(mcc_coord_pair_c auto const& xy_app, MccSkyRADEC_ICRS* radec_icrs)
|
|
{
|
|
double ra, dec;
|
|
|
|
auto err = appToICRS(xy_app.pair_kind, xy_app.epoch(), xy_app.x(), xy_app.y(), &ra, &dec);
|
|
if (!err) {
|
|
if (radec_icrs) {
|
|
radec_icrs->setX(ra);
|
|
radec_icrs->setY(dec);
|
|
}
|
|
}
|
|
|
|
return err;
|
|
}
|
|
|
|
|
|
error_t equationOrigins(mcc_coord_epoch_c auto const& epoch, mcc_angle_c auto* eo)
|
|
{
|
|
if (eo == nullptr) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
|
|
}
|
|
|
|
error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
using real_days_t = std::chrono::duration<double, std::ratio<86400>>;
|
|
|
|
double mjd = epoch.MJD();
|
|
|
|
auto tai_utc = _currentState._leapSeconds[mjd];
|
|
if (tai_utc.has_value()) {
|
|
double tt = mjd;
|
|
tt += std::chrono::duration_cast<real_days_t>(tai_utc.value()).count();
|
|
|
|
auto tt_tai = _currentState._bulletinA.TT_TAI();
|
|
tt += +std::chrono::duration_cast<real_days_t>(tt_tai).count();
|
|
|
|
*eo = eraEo06a(ERFA_DJM0, tt);
|
|
} else {
|
|
ret = MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
// 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)
|
|
template <typename ZAPP_T = std::nullptr_t>
|
|
error_t refractionCorrection(mcc_angle_c auto Zobs, mcc_angle_c auto* dZ, ZAPP_T Zapp = nullptr)
|
|
requires(std::is_null_pointer_v<ZAPP_T> ||
|
|
(std::is_pointer_v<ZAPP_T> && mcc_angle_c<std::remove_pointer_t<ZAPP_T>>))
|
|
{
|
|
error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
if (dZ == nullptr) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
|
|
}
|
|
|
|
refract_model_t rmodel;
|
|
ret = refractionModel(&rmodel);
|
|
|
|
if (!ret) {
|
|
ret = refractionCorrection(rmodel, Zobs, dZ, Zapp);
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
// Zobs must be observed zenithal distance (Zapp = Zobs + dZ -- corrected (in vacuo) zenithal distance)
|
|
template <typename ZAPP_T = std::nullptr_t>
|
|
error_t refractionCorrection(const refract_model_t& rmodel,
|
|
mcc_angle_c auto Zobs,
|
|
mcc_angle_c auto* dZ,
|
|
ZAPP_T Zapp = nullptr)
|
|
requires(std::is_null_pointer_v<ZAPP_T> ||
|
|
(std::is_pointer_v<ZAPP_T> && mcc_angle_c<std::remove_pointer_t<ZAPP_T>>))
|
|
{
|
|
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 constexpr (!std::is_null_pointer_v<ZAPP_T>) {
|
|
*Zapp = Zobs + *dZ;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
// Zapp must be topocentric (in vacuo) zenithal distance (Zobs = Zapp - dZ -- observed, i.e. affected by refraction,
|
|
// zenithal distance)
|
|
template <typename ZOBS_T = std::nullptr_t>
|
|
error_t refractionReverseCorrection(mcc_angle_c auto Zapp, mcc_angle_c auto* dZ, ZOBS_T Zobs = nullptr)
|
|
requires(std::is_null_pointer_v<ZOBS_T> ||
|
|
(std::is_pointer_v<ZOBS_T> && mcc_angle_c<std::remove_pointer_t<ZOBS_T>>))
|
|
{
|
|
error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
if (dZ == nullptr) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
|
|
}
|
|
|
|
refract_model_t rmodel;
|
|
ret = refractionModel(&rmodel);
|
|
|
|
if (!ret) {
|
|
ret = refractionReverseCorrection(rmodel, Zapp, dZ, Zobs);
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
|
|
// Zapp must be topocentric (in vacuo) zenithal distance (Zobs = Zapp - dZ -- observed, i.e. affected by refraction,
|
|
// zenithal distance)
|
|
template <typename ZOBS_T = std::nullptr_t>
|
|
error_t refractionReverseCorrection(const refract_model_t& rmodel,
|
|
mcc_angle_c auto Zapp,
|
|
mcc_angle_c auto* dZ,
|
|
ZOBS_T Zobs = nullptr)
|
|
requires(std::is_null_pointer_v<ZOBS_T> ||
|
|
(std::is_pointer_v<ZOBS_T> && mcc_angle_c<std::remove_pointer_t<ZOBS_T>>))
|
|
{
|
|
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 constexpr (!std::is_null_pointer_v<ZOBS_T>) {
|
|
*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;
|
|
|
|
|
|
error_t icrsTo(bool observed, // true - observed, false - apparent
|
|
mcc_angle_c auto const& ra_icrs,
|
|
mcc_angle_c auto const& dec_icrs,
|
|
mcc_coord_epoch_c auto const& epoch,
|
|
mcc_angle_c auto* ra,
|
|
mcc_angle_c auto* dec,
|
|
mcc_angle_c auto* ha,
|
|
mcc_angle_c auto* az,
|
|
mcc_angle_c auto* zd,
|
|
obj_pars_t* obj_params = nullptr)
|
|
{
|
|
int err;
|
|
double r, d, h, a, z, eo;
|
|
double pressure = 0.0; // 0 for apparent coordinates type (see ERFA's refco.c: if pressure is zero then
|
|
// refraction is also zero)
|
|
|
|
|
|
error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
if (observed) {
|
|
pressure = _currentState.meteo.pressure;
|
|
}
|
|
|
|
auto dut1 = _currentState._bulletinA.DUT1(epoch.MJD());
|
|
|
|
if (!dut1.has_value()) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
|
|
}
|
|
|
|
auto pol_pos = _currentState._bulletinA.polarCoords(epoch.MJD());
|
|
if (!pol_pos.has_value()) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
|
|
}
|
|
|
|
// const auto arcsec2rad = std::numbers::pi / 180 / 3600;
|
|
const auto arcsec2rad = 1.0_arcsecs;
|
|
pol_pos->x *= arcsec2rad;
|
|
pol_pos->y *= arcsec2rad;
|
|
|
|
if (obj_params) {
|
|
err = eraAtco13(ra_icrs, dec_icrs, obj_params->pm_RA, obj_params->pm_DEC, obj_params->parallax,
|
|
obj_params->radvel, ERFA_DJM0, epoch.MJD(), dut1->count(), _currentState.lon,
|
|
_currentState.lat, _currentState.elev, pol_pos->x, pol_pos->y, pressure,
|
|
_currentState.meteo.temperature, _currentState.meteo.humidity, _currentState.wavelength, &a,
|
|
&z, &h, &d, &r, &eo);
|
|
} else {
|
|
err = eraAtco13(ra_icrs, dec_icrs, 0.0, 0.0, 0.0, 0.0, ERFA_DJM0, epoch.MJD(), dut1->count(),
|
|
_currentState.lon, _currentState.lat, _currentState.elev, pol_pos->x, pol_pos->y, pressure,
|
|
_currentState.meteo.temperature, _currentState.meteo.humidity, _currentState.wavelength, &a,
|
|
&z, &h, &d, &r, &eo);
|
|
}
|
|
|
|
if (err == 1) {
|
|
ret = MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR;
|
|
} else if (err == -1) {
|
|
ret = MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE;
|
|
}
|
|
|
|
|
|
if (ra) {
|
|
*ra = r;
|
|
}
|
|
|
|
if (dec) {
|
|
*dec = d;
|
|
}
|
|
|
|
if (ha) {
|
|
*ha = h;
|
|
}
|
|
|
|
if (az) {
|
|
// NOTE: according to definition of astronomical azimuth it is counted from the South through the West, but
|
|
// in the ERFA the azimuth is counted from the North through the East!!!
|
|
//
|
|
*az = MccAngle(a - std::numbers::pi).normalize<MccAngle::NORM_KIND_0_360>();
|
|
// *az = MccAngle(a + std::numbers::pi).normalize<MccAngle::NORM_KIND_0_360>();
|
|
}
|
|
|
|
if (zd) {
|
|
*zd = z;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
error_t toICRS(bool observed, // true - observed, false - apparent
|
|
MccCoordPairKind pair_type,
|
|
mcc_coord_epoch_c auto const& epoch,
|
|
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 ret = MccCCTE_ERFAErrorCode::ERROR_OK;
|
|
|
|
// check coordinate pair consistency
|
|
if (mcc_is_app_coordpair(pair_type) && observed) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR;
|
|
}
|
|
if (mcc_is_obs_coordpair(pair_type) && !observed) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR;
|
|
}
|
|
|
|
std::lock_guard lock{*_stateMutex};
|
|
|
|
auto dut1 = _currentState._bulletinA.DUT1(epoch.MJD());
|
|
|
|
if (!dut1.has_value()) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
|
|
}
|
|
|
|
auto pol_pos = _currentState._bulletinA.polarCoords(epoch.MJD());
|
|
if (!pol_pos.has_value()) {
|
|
return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
|
|
}
|
|
|
|
// const auto arcsec2rad = std::numbers::pi / 180 / 3600;
|
|
const auto arcsec2rad = 1.0_arcsecs;
|
|
pol_pos->x *= arcsec2rad;
|
|
pol_pos->y *= arcsec2rad;
|
|
|
|
std::string type;
|
|
double x, y, ra, dec;
|
|
double pressure = 0.0;
|
|
|
|
if (observed) {
|
|
pressure = _currentState.meteo.pressure;
|
|
}
|
|
|
|
switch (pair_type) {
|
|
case mcc::MccCoordPairKind::COORDS_KIND_AZZD:
|
|
// NOTE: according to definition of astronomical azimuth it is counted from the South through the West,
|
|
// but in the ERFA the azimuth is counted from the North through the East!!!
|
|
//
|
|
x = co_lon + std::numbers::pi;
|
|
y = co_lat;
|
|
type = "A";
|
|
break;
|
|
case mcc::MccCoordPairKind::COORDS_KIND_AZALT:
|
|
// NOTE: according to definition of astronomical azimuth it is counted from the South through the West,
|
|
// but in the ERFA the azimuth is counted from the North through the East!!!
|
|
//
|
|
x = co_lon + std::numbers::pi;
|
|
y = MccCCTE_ERFA::PI_2 - co_lat; // altitude to zenithal distance
|
|
type = "A";
|
|
break;
|
|
case mcc::MccCoordPairKind::COORDS_KIND_HADEC_OBS:
|
|
type = "H";
|
|
x = co_lon;
|
|
y = co_lat;
|
|
break;
|
|
case mcc::MccCoordPairKind::COORDS_KIND_RADEC_OBS:
|
|
type = "R";
|
|
x = co_lon;
|
|
y = co_lat;
|
|
break;
|
|
case mcc::MccCoordPairKind::COORDS_KIND_HADEC_APP:
|
|
type = "H";
|
|
x = co_lon;
|
|
y = co_lat;
|
|
break;
|
|
case mcc::MccCoordPairKind::COORDS_KIND_RADEC_APP:
|
|
type = "R";
|
|
x = co_lon;
|
|
y = co_lat;
|
|
break;
|
|
default:
|
|
return MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR;
|
|
};
|
|
|
|
int err =
|
|
eraAtoc13(type.c_str(), x, y, ERFA_DJM0, epoch.MJD(), dut1->count(), _currentState.lon, _currentState.lat,
|
|
_currentState.elev, pol_pos->x, pol_pos->y, pressure, _currentState.meteo.temperature,
|
|
_currentState.meteo.humidity, _currentState.wavelength, &ra, &dec);
|
|
|
|
if (err == 1) {
|
|
ret = MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR;
|
|
} else if (err == -1) {
|
|
ret = MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE;
|
|
}
|
|
|
|
if (ra) {
|
|
*ra_icrs = ra;
|
|
}
|
|
|
|
if (dec) {
|
|
*dec_icrs = dec;
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
};
|
|
|
|
|
|
} // namespace mcc::ccte::erfa
|