diff --git a/mcc/mcc_ccte_erfa_new.h b/mcc/mcc_ccte_erfa_new.h index c8aab44..a73954a 100644 --- a/mcc/mcc_ccte_erfa_new.h +++ b/mcc/mcc_ccte_erfa_new.h @@ -251,35 +251,145 @@ public: // 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) + // 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 icrs2obs(MccSkyRADEC_ICRS const& radec_icrs, - MccSkyRADEC_OBS* radec_obs, - MccSkyAZZD* azzd, - mcc_angle_c auto* ha_obs) { + 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; + MccCelestialCoordEpoch ep; + ep.fromMJD((radec_obs->MJD())); + + auto err = icrsToObs(radec_icrs.x(), radec_icrs.y(), ep, &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->setMJD(radec_obs->MJD()); + azzd->setX(az); + azzd->setY(zd); + } + + if (ha_obs) { + *ha_obs = ha; + } + } + + return err; }; - 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) + + // 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 obs2icrs(mcc_coord_pair_c auto const& xy_obs, MccSkyRADEC_ICRS* radec_icrs) {} + 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; + + MccCelestialCoordEpoch ep; + ep.fromMJD((radec_app->MJD())); + + auto err = icrsToApp(radec_icrs.x(), radec_icrs.y(), ep, &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->setMJD(radec_app->MJD()); + 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; + + MccCelestialCoordEpoch ep; + ep.fromMJD((xy_obs->MJD())); + + auto err = obsToICRS(xy_obs.pair_kind, ep, 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 equationOrigins(const mcc_julday_c auto& mjd, mcc_angle_c auto* eo) @@ -415,6 +525,197 @@ protected: engine_state_t _currentState{}; std::unique_ptr _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 (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(); + } + + 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; + } };