...
This commit is contained in:
130
mcc/mcc_coord.h
130
mcc/mcc_coord.h
@@ -331,7 +331,58 @@ protected:
|
||||
template <mcc_coord_pair_c PT>
|
||||
auto toHelper(PT& cpair)
|
||||
{
|
||||
double ra_icrs, dec_icrs, ra_obs, dec_obs, ha_obs, az_obs, zd_obs, lst, eo;
|
||||
static constexpr double half_pi = std::numbers::pi / 2.0;
|
||||
|
||||
// HA, DEC to AZ, ALT (AZ from the South through the West)
|
||||
auto hadec2azalt = [](double ha, double dec, double phi, double& az, double& alt) {
|
||||
const auto cos_phi = std::cos(phi), sin_phi = std::sin(phi);
|
||||
const auto cos_dec = std::cos(dec), sin_dec = std::sin(dec);
|
||||
const auto cos_ha = std::cos(ha), sin_ha = std::sin(ha);
|
||||
|
||||
auto x = sin_phi * cos_dec * cos_ha - cos_phi * sin_dec;
|
||||
auto y = -cos_dec * sin_ha;
|
||||
auto z = cos_phi * cos_dec * cos_ha + sin_phi * sin_dec;
|
||||
|
||||
auto xx = x * x, yy = y * y;
|
||||
decltype(x) r;
|
||||
if (xx < yy) {
|
||||
r = yy * sqrt(1.0 + xx / yy);
|
||||
} else {
|
||||
r = xx * sqrt(1.0 + yy / xx);
|
||||
}
|
||||
|
||||
az = utils::isEqual(r, 0.0) ? 0.0 : std::atan2(y, x);
|
||||
if (az < 0.0) {
|
||||
az += std::numbers::pi * 2.0; // to range of [0, 2*PI]
|
||||
}
|
||||
|
||||
alt = std::atan2(z, r);
|
||||
};
|
||||
|
||||
// AZ, ALT to HA, DEC (AZ from the South through the West)
|
||||
auto azalt2hadec = [](double az, double alt, double phi, double& ha, double& dec) {
|
||||
const auto cos_phi = std::cos(phi), sin_phi = std::sin(phi);
|
||||
const auto cos_az = std::cos(az), sin_az = std::sin(az);
|
||||
const auto cos_alt = std::cos(alt), sin_alt = std::sin(alt);
|
||||
|
||||
auto x = sin_phi * cos_alt * cos_az + cos_phi * sin_alt;
|
||||
auto y = cos_alt * sin_az;
|
||||
auto z = -cos_phi * cos_alt * cos_az + sin_phi * sin_alt;
|
||||
|
||||
auto xx = x * x, yy = y * y;
|
||||
decltype(x) r;
|
||||
if (xx < yy) {
|
||||
r = yy * sqrt(1.0 + xx / yy);
|
||||
} else {
|
||||
r = xx * sqrt(1.0 + yy / xx);
|
||||
}
|
||||
|
||||
ha = utils::isEqual(r, 0.0) ? 0.0 : std::atan2(y, x);
|
||||
dec = std::atan2(z, r);
|
||||
};
|
||||
|
||||
double phi;
|
||||
double ra_icrs, dec_icrs, ra, dec, ha, az, zd, alt, lst, eo;
|
||||
|
||||
static_assert(PT::pairKind == MccCoordPairKind::COORDS_KIND_GENERIC, "UNSUPPORTED SKY POINT TRANSFORMATION!");
|
||||
static_assert(PT::pairKind == MccCoordPairKind::COORDS_KIND_UNKNOWN, "UNSUPPORTED SKY POINT TRANSFORMATION!");
|
||||
@@ -373,7 +424,7 @@ protected:
|
||||
}
|
||||
}
|
||||
|
||||
// here, APP or OBS to ICRS
|
||||
// here, from APP or OBS to ICRS
|
||||
if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS &&
|
||||
PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
|
||||
cpair = PT(PT::x_t(ra_icrs), PT::y_t(dec_icrs), MccCelestialCoordEpoch{});
|
||||
@@ -381,32 +432,105 @@ protected:
|
||||
}
|
||||
|
||||
if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
|
||||
if constexpr (mccIsObsCoordPairKind<PT::pairKind>) {
|
||||
// cct_engine.icrsToObs(...)
|
||||
} else if constexpr (mccIsAppCoordPairKind<PT::pairKind>) {
|
||||
// cct_engine.icrsToApp(...)
|
||||
} else {
|
||||
static_assert(false, "UNSUPPORTED SKY POINT TRANSFORMATION!");
|
||||
}
|
||||
|
||||
|
||||
if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP) {
|
||||
// cct_engine.icrsToApp(...)
|
||||
cpair.setX(ra);
|
||||
cpair.setY(dec);
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) {
|
||||
// cct_engine.icrsToObs(...)
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) {
|
||||
// cct_engine.icrsToApp(...)
|
||||
cpair.setX(ha);
|
||||
cpair.setY(dec);
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
|
||||
// cct_engine.icrsToObs(...)
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) {
|
||||
// cct_engine.icrsToObs(...)
|
||||
cpair.setX(az);
|
||||
cpair.setY(zd);
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) {
|
||||
// cct_engine.icrsToObs(...)
|
||||
cpair.setX(az);
|
||||
cpair.setY(half_pi - zd);
|
||||
} else {
|
||||
static_assert(false, "UNSUPPORTED SKY POINT TRANSFORMATION!");
|
||||
}
|
||||
} else if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_APP) {
|
||||
// compute EO, LST, refract model
|
||||
ha = lst + eo - cpair.x(); // from CIO based RA
|
||||
|
||||
if constexpr (PT::pairKind != MccCoordPairKind::COORDS_KIND_HADEC_APP) {
|
||||
hadec2azalt(ha, _y, phi, az, alt); // to app az, alt
|
||||
// to observed zenithal distance: alt += z_corr
|
||||
}
|
||||
|
||||
if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) {
|
||||
azalt2hadec(az, alt, phi, ha, dec); // to obs ha, dec
|
||||
|
||||
ra = lst + eo - ha; // CIO based RA
|
||||
cpair.setX(ra);
|
||||
cpair.setY(dec);
|
||||
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) {
|
||||
// ha = lst + eo - ra?
|
||||
cpair.setX(ha);
|
||||
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
|
||||
azalt2hadec(az, alt, phi, ha, dec); // to obs ha, dec
|
||||
cpair.setX(ha);
|
||||
cpair.setY(dec);
|
||||
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) {
|
||||
cpair.setX(az); // ????????????!!!!!!!!!
|
||||
cpair.setY(half_pi - alt);
|
||||
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) {
|
||||
cpair.setX(az); // ????????????!!!!!!!!!
|
||||
cpair.setY(alt);
|
||||
} else {
|
||||
static_assert(false, "UNSUPPORTED SKY POINT TRANSFORMATION!");
|
||||
}
|
||||
} else if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) {
|
||||
ha = lst + eo - cpair.x(); // from CIO based RA
|
||||
|
||||
if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP) {
|
||||
hadec2azalt(ha, _y, phi, az, alt);
|
||||
// to apparent zenithal distance: alt -= z_corr;
|
||||
azalt2hadec(az, alt, phi, ha, dec); // to app ha,dec
|
||||
|
||||
ra = lst + eo - ha; // CIO based RA
|
||||
cpair.setX(ra);
|
||||
cpair.setY(dec);
|
||||
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) {
|
||||
hadec2azalt(ha, _y, phi, az, alt);
|
||||
// to apparent zenithal distance: alt -= z_corr;
|
||||
azalt2hadec(az, alt, phi, ha, dec); // to app ha,dec
|
||||
|
||||
cpair.setX(ha);
|
||||
cpair.setY(dec);
|
||||
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
|
||||
cpair.setX(ha);
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) {
|
||||
hadec2azalt(ha, _y, phi, az, alt);
|
||||
cpair.setX(az);
|
||||
cpair.setY(half_pi - alt);
|
||||
} else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) {
|
||||
hadec2azalt(ha, _y, phi, az, alt);
|
||||
cpair.setX(az);
|
||||
cpair.setY(alt);
|
||||
} else {
|
||||
static_assert(false, "UNSUPPORTED SKY POINT TRANSFORMATION!");
|
||||
}
|
||||
} else if (pkind == MccCoordPairKind::COORDS_KIND_HADEC_APP) {
|
||||
} else if (pkind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
|
||||
} else if (pkind == MccCoordPairKind::COORDS_KIND_AZZD) {
|
||||
|
||||
Reference in New Issue
Block a user