306 lines
9.4 KiB
C++
306 lines
9.4 KiB
C++
#pragma once
|
|
|
|
#include "mcc_coord.h"
|
|
#include "mount_astrom.h"
|
|
|
|
namespace mcc
|
|
{
|
|
|
|
namespace traits
|
|
{
|
|
|
|
template <typename T>
|
|
concept mcc_prohibited_zone_c = requires(T t, const T const_t) {
|
|
typename T::pzcontext_t;
|
|
|
|
{ const_t.name() } -> std::same_as<std::string_view>;
|
|
{ const_t.desc() } -> std::same_as<std::string_view>;
|
|
|
|
// check if given coordinates are within the zone
|
|
{
|
|
t.inZone(std::declval<const MccAngle&>(), std::declval<const MccAngle&>(),
|
|
std::declval<typename T::pzcontext_t>())
|
|
} -> std::convertible_to<bool>;
|
|
|
|
// a time duration to reach the zone (0 - if already in the zone, chrono::duration<>::max() if never
|
|
// reach the zone)
|
|
{
|
|
t.timeTo(std::declval<const MccAngle&>(), std::declval<const MccAngle&>(),
|
|
std::declval<typename T::pzcontext_t>(), std::declval<const std::chrono::system_clock::time_point&>())
|
|
} -> mcc_time_duration_c;
|
|
|
|
// a time duration to exit the zone (0 - if already out of the zone, chrono::duration<>::max() if
|
|
// never exit the zone)
|
|
{
|
|
t.timeFrom(std::declval<const MccAngle&>(), std::declval<const MccAngle&>(),
|
|
std::declval<typename T::pzcontext_t>(),
|
|
std::declval<const std::chrono::system_clock::time_point&>())
|
|
} -> mcc_time_duration_c;
|
|
};
|
|
|
|
} // namespace traits
|
|
|
|
|
|
class MccProhibitedZone
|
|
{
|
|
public:
|
|
enum pzcoords_kind_t { COORDS_KIND_RADEC_IRCS, COORDS_KIND_RADEC_APP, COORDS_KIND_HADEC_APP, COORDS_KIND_AZALT };
|
|
|
|
struct pzcontext_t {
|
|
typedef std::chrono::duration<double> real_secs_t; // seconds duration in double
|
|
|
|
pzcoords_kind_t coords_kind;
|
|
|
|
real_secs_t dut1; // UT1-UTC
|
|
real_secs_t tt_tai; // TT-TAI
|
|
real_secs_t tai_utc; // TAI-UTC (leap seconds)
|
|
|
|
MccAngle lat, lon; // site geographic coordinates
|
|
};
|
|
|
|
virtual ~MccProhibitedZone() = default;
|
|
};
|
|
|
|
class MccMinAltPZ : public MccProhibitedZone
|
|
{
|
|
public:
|
|
MccMinAltPZ(const MccAngle& min_alt) : _minAlt(min_alt) {}
|
|
|
|
MccAngle minAlt() const
|
|
{
|
|
return _minAlt;
|
|
}
|
|
|
|
private:
|
|
double _minAlt;
|
|
|
|
bool inZoneImpl(const MccAngle& alt, const MccAngle&)
|
|
{
|
|
return alt <= _minAlt;
|
|
}
|
|
};
|
|
|
|
|
|
class MccMaxAltPZ
|
|
{
|
|
public:
|
|
MccMaxAltPZ(double max_alt) : _maxAlt(max_alt) {}
|
|
|
|
double maxAlt() const
|
|
{
|
|
return _maxAlt;
|
|
}
|
|
|
|
private:
|
|
double _maxAlt;
|
|
|
|
bool inZoneImpl(const MccAngle& alt, const MccAngle&)
|
|
{
|
|
return alt >= _maxAlt;
|
|
}
|
|
};
|
|
|
|
|
|
enum class MccAltLimitKind { MIN_ALT_LIMIT, MAX_ALT_LIMIT };
|
|
|
|
template <MccAltLimitKind KIND = MccAltLimitKind::MIN_ALT_LIMIT>
|
|
class MccAltLimitPZ
|
|
{
|
|
public:
|
|
static constexpr MccAltLimitKind altLimitKind = KIND;
|
|
|
|
MccAltLimitPZ(const MccAngle& alt_limit) : _altLimit(alt_limit)
|
|
{
|
|
_altLimit.normalize<MccAngle::NORM_KIND_90_90>();
|
|
}
|
|
|
|
std::string_view name() const
|
|
{
|
|
return KIND == MccAltLimitKind::MIN_ALT_LIMIT ? "MINALT-ZONE"
|
|
: KIND == MccAltLimitKind::MAX_ALT_LIMIT ? "MAXALT-ZONE"
|
|
: "ALTLIMIT-UNKNOWN";
|
|
}
|
|
|
|
std::string_view desc() const
|
|
{
|
|
return KIND == MccAltLimitKind::MIN_ALT_LIMIT ? "Minimal altitude prohibited zone"
|
|
: KIND == MccAltLimitKind::MAX_ALT_LIMIT ? "Maximal altitude prohibited zone"
|
|
: "Unknown altitude prohibited zone";
|
|
}
|
|
|
|
bool inZone(const MccAngle& x, const MccAngle& y, const MccProhibitedZone::pzcontext_t& context)
|
|
{
|
|
if (context.coords_kind == MccProhibitedZone::COORDS_KIND_AZALT) { // trivial case
|
|
if constexpr (KIND == MccAltLimitKind::MIN_ALT_LIMIT) {
|
|
return y <= _altLimit;
|
|
} else if constexpr (KIND == MccAltLimitKind::MAX_ALT_LIMIT) {
|
|
return y >= _altLimit;
|
|
}
|
|
} else if (context.coords_kind == MccProhibitedZone::COORDS_KIND_RADEC_APP) {
|
|
auto dd =
|
|
astrom::mcc_time_to_alt(_altLimit, x, y, context.lat, context.lon, std::chrono::system_clock::now(),
|
|
context.dut1, context.tt_tai, context.tai_utc);
|
|
}
|
|
}
|
|
|
|
|
|
auto timeTo(const MccAngle& x,
|
|
const MccAngle& y,
|
|
const MccProhibitedZone::pzcontext_t& context,
|
|
traits::mcc_systime_c auto const& utc = std::chrono::system_clock::now())
|
|
{
|
|
if (context.coords_kind == MccProhibitedZone::COORDS_KIND_RADEC_APP) {
|
|
auto dd = astrom::mcc_time_to_alt(_altLimit, x, y, context.lat, context.lon, utc, context.dut1,
|
|
context.tt_tai, context.tai_utc);
|
|
if (std::isnan(dd.first)) { // error!
|
|
throw std::system_error(std::make_error_code(std::errc::invalid_argument));
|
|
}
|
|
|
|
if (std::isinf(dd.first)) { // never reach zone
|
|
return dd.first;
|
|
}
|
|
|
|
|
|
if constexpr (KIND == MccAltLimitKind::MIN_ALT_LIMIT) {
|
|
if (dd.first <= dd.second) { // in zone
|
|
return decltype(dd.first)(0.0);
|
|
}
|
|
} else if constexpr (KIND == MccAltLimitKind::MAX_ALT_LIMIT) {
|
|
if (dd.first >= dd.second) { // in zone
|
|
return decltype(dd.first)(0.0);
|
|
}
|
|
}
|
|
|
|
return dd.second;
|
|
} else {
|
|
throw std::system_error(std::make_error_code(std::errc::operation_not_supported));
|
|
}
|
|
}
|
|
|
|
auto timeFrom(const MccAngle& x,
|
|
const MccAngle& y,
|
|
const MccProhibitedZone::pzcontext_t& context,
|
|
traits::mcc_systime_c auto const& utc = std::chrono::system_clock::now())
|
|
{
|
|
if (context.coords_kind == MccProhibitedZone::COORDS_KIND_RADEC_APP) {
|
|
auto dd = astrom::mcc_time_to_alt(_altLimit, x, y, context.lat, context.lon, utc, context.dut1,
|
|
context.tt_tai, context.tai_utc);
|
|
if (std::isnan(dd.first)) { // error!
|
|
throw std::system_error(std::make_error_code(std::errc::invalid_argument));
|
|
}
|
|
|
|
if (std::isinf(dd.first)) { // never reach zone
|
|
return dd.first;
|
|
}
|
|
|
|
if constexpr (KIND == MccAltLimitKind::MIN_ALT_LIMIT) {
|
|
if (dd.first > dd.second) { // not in zone
|
|
return decltype(dd.first)(0.0);
|
|
}
|
|
} else if constexpr (KIND == MccAltLimitKind::MAX_ALT_LIMIT) {
|
|
if (dd.first < dd.second) { // not in zone
|
|
return decltype(dd.first)(0.0);
|
|
}
|
|
}
|
|
|
|
return dd.first;
|
|
} else {
|
|
throw std::system_error(std::make_error_code(std::errc::operation_not_supported));
|
|
}
|
|
}
|
|
|
|
private:
|
|
MccAngle _altLimit;
|
|
};
|
|
|
|
|
|
/*
|
|
* a general planar ellipse equation:
|
|
* A*(x-xc)^2 + B*(x-xc)(y-yc) + C*(y-yc)^2 = 1
|
|
*
|
|
* A = cos(t)^2/a^2 + sin(t)^2/b^2
|
|
* B = sin(2*t)/a^2 - sin(2*t)/b^2
|
|
* C = cos(t)^2/b^2 + sin(t)^2/a^2
|
|
*
|
|
* t - angle between X-axis and big semi-axis (a)
|
|
*
|
|
*
|
|
* Ellipse on unit sphere (Az, Alt):
|
|
* x^2/a^2 + y^2/b^2 = 1,
|
|
* where a = tan(alpha), b = tan(beta),
|
|
* a and b - big and small spherical semi-axis
|
|
* also:
|
|
* let delta is a spherical distance between ellipse focuses, then
|
|
* d = tan(delta) and b^2 = (a^2 - d^2)/(1 + d^2)
|
|
*
|
|
* tangent coordinates:
|
|
* x = tan(Az)
|
|
* y = tan(Alt)*sqrt(1+tan(Az)^2)
|
|
*
|
|
*
|
|
* --------------------------------------------
|
|
*
|
|
* let P - point with (Az_P, Zd_P),
|
|
* (Az_C, Zd_C) - center of ellipse, a and b big and small semi-axis,
|
|
* vec_a and vec_b - unit vectors along a an b (it lie in tangent surface to point C!!!), then
|
|
*
|
|
* ((vec_P-vec_C)*vec_b)^2/a^2 + ((vec_P-vec_C)*vec_b)^2/b^2 <= 1.0
|
|
* {((vec_P-vec_C)*vec_b)^2/tan(a)^2 + ((vec_P-vec_C)*vec_b)^2/tan(b)^2 <= 1.0}
|
|
* * - dot-product
|
|
*
|
|
* vec_P = (sin(Zd_P)*cos(Az_P), sin(Zd_P)*sin(Az_P), cos(Zd_P))
|
|
* vec_C = (sin(Zd_C)*cos(Az_C), sin(Zd_C)*sin(Az_C), cos(Zd_C))
|
|
*
|
|
*/
|
|
class MccEllipsePZ
|
|
{
|
|
public:
|
|
MccEllipsePZ(const MccAngle& xc,
|
|
const MccAngle& yc,
|
|
const MccAngle& a,
|
|
const MccAngle& b,
|
|
const MccAngle& theta = 0.0)
|
|
: _xc(xc), _yc(yc), _a(a), _b(b), _theta(theta), _tanXc(std::tan(xc)), _tanYc(std::tan(yc))
|
|
{
|
|
_tanYc *= std::sqrt(1.0 + _tanXc * _tanXc);
|
|
|
|
auto _tan2A = tan(a);
|
|
auto _tan2B = tan(b);
|
|
|
|
_tan2A *= _tan2A;
|
|
_tan2B *= _tan2B;
|
|
|
|
auto ct = cos(theta);
|
|
auto ct2 = ct * ct;
|
|
auto st = sin(theta);
|
|
auto st2 = st * st;
|
|
auto sin2t = sin(2.0 * theta);
|
|
|
|
|
|
cxx = ct2 / _tan2A + st2 / _tan2B;
|
|
cyy = st2 / _tan2A + ct2 / _tan2B;
|
|
|
|
cxy = sin2t / _tan2A - sin2t / _tan2B;
|
|
}
|
|
|
|
// circle
|
|
MccEllipsePZ(const MccAngle& xc, const MccAngle& yc, const MccAngle& a) : mcc::MccEllipsePZ(xc, yc, a, a) {}
|
|
|
|
|
|
private:
|
|
double _xc, _yc, _a, _b, _theta;
|
|
double _tanXc, _tanYc, cxx, cxy, cyy;
|
|
|
|
bool inZoneImpl(const MccAngle& x, const MccAngle& y)
|
|
{
|
|
auto tanX = tan(x);
|
|
auto tanY = tan(y) * sqrt(1.0 + tanX * tanX);
|
|
auto xr = tanX - _tanXc;
|
|
auto yr = tanY - _tanYc;
|
|
|
|
return (cxx * xr * xr + cxy * xr * yr + cyy * yr * yr) <= 1.0;
|
|
}
|
|
};
|
|
|
|
} // namespace mcc
|