396 lines
12 KiB
C++
396 lines
12 KiB
C++
#pragma once
|
|
|
|
/*********************************
|
|
* MOUNT CONTROL COMPONENTS *
|
|
* *
|
|
* astrometry functions *
|
|
*********************************/
|
|
|
|
#include <chrono>
|
|
#include <xsimd/xsimd.hpp>
|
|
#include "utils.h"
|
|
|
|
|
|
namespace mcc::traits
|
|
{
|
|
|
|
template <typename T>
|
|
concept mcc_scalar_or_simd_c = xsimd::is_batch<T>::value || std::is_arithmetic_v<T>;
|
|
|
|
|
|
template <typename T>
|
|
concept mcc_real_scalar_or_real_range_c =
|
|
std::floating_point<T> || std::ranges::output_range<T, float> || std::ranges::output_range<T, double>;
|
|
|
|
template <typename T>
|
|
concept mcc_real_or_char_range =
|
|
std::floating_point<T> ||
|
|
(std::ranges::contiguous_range<T> && std::same_as<char, std::remove_cvref_t<std::ranges::range_value_t<T>>>);
|
|
|
|
} // namespace mcc::traits
|
|
|
|
|
|
namespace mcc::astro
|
|
{
|
|
|
|
// modified Julian date (based on ERFA eraCal2jd)
|
|
template <traits::mcc_real_scalar_or_real_range_c ResT>
|
|
static int mcc_julday(const std::chrono::system_clock::time_point& start_time,
|
|
ResT& mjd,
|
|
const std::chrono::system_clock::duration& step = std::chrono::milliseconds(100))
|
|
{
|
|
size_t mjd_size = 0;
|
|
if constexpr (std::ranges::range<ResT>) {
|
|
mjd_size = std::ranges::distance(mjd.begin(), mjd.end());
|
|
if (!mjd_size) {
|
|
return -100;
|
|
}
|
|
}
|
|
|
|
using namespace std::literals::chrono_literals;
|
|
|
|
auto dd = std::chrono::floor<std::chrono::days>(start_time);
|
|
std::chrono::year_month_day ymd{dd};
|
|
|
|
static constexpr std::chrono::year MIN_YEAR = -4799y;
|
|
|
|
if (ymd.year() < MIN_YEAR) {
|
|
return -1;
|
|
}
|
|
if (!ymd.month().ok()) {
|
|
return -2;
|
|
}
|
|
|
|
// my = (im - 14) / 12;
|
|
// iypmy = (long) (iy + my);
|
|
|
|
int64_t my = -(14 - (unsigned)ymd.month()) / 12;
|
|
int64_t iypmy = (int)ymd.year() + my;
|
|
|
|
// (1461L * (iypmy + 4800L)) / 4L
|
|
// + (367L * (long) (im - 2 - 12 * my)) / 12L
|
|
// - (3L * ((iypmy + 4900L) / 100L)) / 4L
|
|
// + (long) id - 2432076L
|
|
|
|
// integer part of result MJD
|
|
int64_t mjd_int = 1461LL * (iypmy + 480LL) / 4LL +
|
|
(367LL * ((int64_t)(unsigned)ymd.month() - 2LL - 12LL * my)) / 12LL -
|
|
(3LL * (iypmy + 4900LL) / 100LL) / 4LL + (int64_t)(unsigned)ymd.day() - 2432076LL;
|
|
|
|
constexpr double nanosec = 1.0 / 24.0 / 3600.0 / 1.0E-9; // 1 nanosecond in days
|
|
|
|
double mjd_float = static_cast<double>(mjd_int) +
|
|
std::chrono::duration_cast<std::chrono::nanoseconds>(start_time - dd).count() * nanosec;
|
|
|
|
if constexpr (std::ranges::range<ResT>) {
|
|
double d_step = std::chrono::duration_cast<std::chrono::nanoseconds>(step).count() * nanosec;
|
|
|
|
size_t i = 0;
|
|
|
|
#ifdef VEC_XSIMD
|
|
constexpr size_t reg_size = xsimd::batch<double>::size;
|
|
size_t vec_size = mjd_size - mjd_size % reg_size;
|
|
|
|
xsimd::batch<double> res_reg{mjd_float};
|
|
xsimd::batch<double> step_reg = [d_step]<size_t... Is>(std::index_sequence<Is...>) {
|
|
return xsimd::batch<double>{(Is * d_step)...};
|
|
}(std::make_index_sequence<reg_size>{});
|
|
|
|
alignas(xsimd::batch<double>::arch_type::alignment()) double arr[reg_size];
|
|
|
|
auto ptr = mjd.begin();
|
|
|
|
// vectorized part
|
|
for (; i < vec_size; i += vec_size) {
|
|
res_reg += step_reg;
|
|
|
|
if constexpr (std::ranges::contiguous_range<ResT>) {
|
|
res_reg.store_unaligned(mjd.data() + i);
|
|
// res_reg.store_aligned(mjd.data() + i);
|
|
} else {
|
|
res_reg.store_aligned(arr);
|
|
std::ranges::copy(arr, ptr);
|
|
|
|
if constexpr (std::ranges::random_access_range<ResT>) {
|
|
ptr += reg_size;
|
|
} else {
|
|
for (int k = 0; k < reg_size; ++k) {
|
|
++ptr;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
// scalar part
|
|
for (size_t j = i; j < mjd_size; ++j, ++ptr) {
|
|
*ptr = mjd_float + j * d_step;
|
|
}
|
|
} else { // result is scalar
|
|
mjd = mjd_float;
|
|
}
|
|
|
|
|
|
return 0;
|
|
}
|
|
|
|
/*
|
|
* angles are in degrees or sexagimal string form
|
|
*
|
|
* returns
|
|
* NaN if object is non-rising or "alt_limit" < 0, Inf is circumpolar
|
|
*/
|
|
template <traits::mcc_real_or_char_range AT,
|
|
traits::mcc_real_or_char_range RT,
|
|
traits::mcc_real_or_char_range DT,
|
|
traits::mcc_real_or_char_range LT>
|
|
double mcc_time_to_alt_limit(const AT& alt_limit,
|
|
const RT& RA,
|
|
const DT& DEC,
|
|
const LT& LAT,
|
|
const std::chrono::system_clock::time_point& now)
|
|
{
|
|
// sin(alt) = sin(DEC)*sin(phi) + cos(DEC)*cos(phi)*cos(HA)
|
|
// HA = LST - RA
|
|
|
|
double ra, dec, lat, alt;
|
|
|
|
if constexpr (std::floating_point<AT>) {
|
|
alt = alt_limit * utils::deg2radCoeff;
|
|
} else {
|
|
alt = utils::parsAngleString(alt_limit);
|
|
alt *= utils::deg2radCoeff;
|
|
}
|
|
|
|
if (alt < 0.0) {
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
}
|
|
|
|
if constexpr (std::floating_point<RT>) {
|
|
ra = RA * utils::deg2radCoeff;
|
|
} else {
|
|
ra = utils::parsAngleString(RA, true);
|
|
ra *= utils::deg2radCoeff;
|
|
}
|
|
|
|
if constexpr (std::floating_point<DT>) {
|
|
dec = DEC * utils::deg2radCoeff;
|
|
} else {
|
|
dec = utils::parsAngleString(DEC);
|
|
dec *= utils::deg2radCoeff;
|
|
}
|
|
|
|
if constexpr (std::floating_point<LT>) {
|
|
lat = LAT * utils::deg2radCoeff;
|
|
} else {
|
|
lat = utils::parsAngleString(LAT);
|
|
lat *= utils::deg2radCoeff;
|
|
}
|
|
|
|
if (lat >= 0.0) { // north hemisphere
|
|
if (dec < (lat - std::numbers::pi / 2.0)) { // never rises above horizont
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
}
|
|
} else { // south hemisphere
|
|
if (dec > (lat + std::numbers::pi / 2.0)) { // never rises above horizont
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
}
|
|
}
|
|
|
|
double cos_ha = (std::sin(alt) - std::sin(dec) * std::sin(lat)) / std::cos(dec) / std::cos(lat);
|
|
if (std::abs(cos_ha) > 1.0) { // circumpolar (it never sets below horizon)
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
double lst = std::acos(cos_ha) + ra;
|
|
|
|
|
|
return 0.0;
|
|
}
|
|
|
|
|
|
/* IERS bulletins
|
|
*
|
|
* BULLETIN A: https://datacenter.iers.org/data/latestVersion/bulletinA.txt
|
|
* leapseconds: https://hpiers.obspm.fr/iers/bul/bulc/Leap_Second.dat
|
|
*
|
|
*/
|
|
|
|
enum iers_db_state_t { IERS_DB_STATE_UNINITIALIZED, IERS_DB_STATE_OK, IERS_DB_STATE_PARSE_ERROR };
|
|
|
|
struct leapsecond_db_elem_t {
|
|
double mjd;
|
|
unsigned day, month;
|
|
int year;
|
|
// std::chrono::year_month_day ymd;
|
|
double tai_utc; // TAI-UTC in seconds
|
|
};
|
|
|
|
struct leapsecond_db_t {
|
|
iers_db_state_t state = IERS_DB_STATE_UNINITIALIZED;
|
|
std::chrono::system_clock::time_point expireDate{};
|
|
std::vector<leapsecond_db_elem_t> db{};
|
|
};
|
|
|
|
|
|
|
|
struct earth_orient_db_elem_t {
|
|
int year;
|
|
unsigned month, day;
|
|
double mjd;
|
|
double x, y; // Polar coordinates in arcsecs
|
|
double dut1; // UT1-UTC in seconds
|
|
};
|
|
|
|
|
|
struct earth_orient_db_t {
|
|
iers_db_state_t state = IERS_DB_STATE_UNINITIALIZED;
|
|
std::chrono::system_clock::time_point bulletinDate{};
|
|
double tt_tai = 0.0; // TT-TAI
|
|
std::vector<earth_orient_db_elem_t> db{};
|
|
};
|
|
|
|
|
|
|
|
static earth_orient_db_t mcc_parse_bulletinA(std::derived_from<std::basic_istream<char>> auto& stream,
|
|
char comment_sym = '*')
|
|
{
|
|
const std::regex bull_date_rx{
|
|
"^ *[0-9]{1,2} +(January|February|March|April|May|June|July|August|September|October|November|December) "
|
|
"+[0-9]{4,} +Vol\\. +[XMLCDVI]+ +No\\. +[0-9]+ *$"};
|
|
|
|
const std::regex bull_tt_tai_rx{"^ *TT += +TAI +\\+ +[0-9]+\\.[0-9]+ +seconds *$"};
|
|
|
|
const std::regex bull_tab_title_rx{"^ *MJD +x\\(arcsec\\) +y\\(arcsec\\) +UT1-UTC\\(sec\\) *$"};
|
|
|
|
// 2025 3 7 60741 0.0663 0.3341 0.04348
|
|
const std::regex bull_tab_vals_rx{
|
|
"^ *[0-9]{4,} +[0-9]{1,2} +[0-9]{1,2} +[0-9]{5,} +[0-9]+\\.[0-9]+ +[0-9]+\\.[0-9]+ +[0-9]+\\.[0-9]+ *$"};
|
|
|
|
earth_orient_db_t db;
|
|
enum { TAB_STATE_SEEK, TAB_STATE_START };
|
|
int tab_state = TAB_STATE_SEEK;
|
|
|
|
int year;
|
|
unsigned month, day;
|
|
double mjd, x, y, dut1;
|
|
std::istringstream is;
|
|
|
|
for (std::string line; std::getline(stream, line);) {
|
|
if (line.empty()) {
|
|
continue;
|
|
}
|
|
|
|
auto sv = utils::trimSpaces(line, utils::TrimType::TRIM_LEFT);
|
|
|
|
if (sv.size()) {
|
|
if (sv[0] == comment_sym) { // comment string
|
|
continue;
|
|
}
|
|
|
|
if (tab_state == TAB_STATE_START) {
|
|
if (std::regex_match(sv.begin(), sv.end(), bull_tab_vals_rx)) {
|
|
// is.str({sv.begin(), sv.end()});
|
|
is.str(line);
|
|
is >> year >> month >> day >> mjd >> x >> y >> dut1;
|
|
db.db.emplace_back(year, month, day, mjd, x, y, dut1);
|
|
is.clear();
|
|
} else { // end of the table - just stop parsing
|
|
break;
|
|
}
|
|
|
|
continue;
|
|
}
|
|
|
|
if (std::regex_match(sv.begin(), sv.end(), bull_date_rx)) {
|
|
is.str({sv.begin(), sv.end()});
|
|
is >> std::chrono::parse("%d %B %Y", db.bulletinDate);
|
|
continue;
|
|
}
|
|
|
|
if (std::regex_match(sv.begin(), sv.end(), bull_tt_tai_rx)) {
|
|
is.str({sv.begin(), sv.end()});
|
|
std::string dummy;
|
|
is >> dummy >> dummy >> dummy >> dummy >> db.tt_tai;
|
|
continue;
|
|
}
|
|
|
|
if (std::regex_match(sv.begin(), sv.end(), bull_tab_title_rx)) {
|
|
tab_state = TAB_STATE_START;
|
|
continue;
|
|
}
|
|
|
|
} else { // empty string (only spaces)
|
|
continue;
|
|
}
|
|
}
|
|
|
|
if (db.db.empty()) {
|
|
db.state = IERS_DB_STATE_PARSE_ERROR;
|
|
} else {
|
|
db.state = IERS_DB_STATE_OK;
|
|
}
|
|
|
|
return db;
|
|
}
|
|
|
|
|
|
static leapsecond_db_t mcc_parse_leapsecs(std::derived_from<std::basic_istream<char>> auto& stream,
|
|
char comment_sym = '#')
|
|
{
|
|
// # File expires on 28 December 2025
|
|
const std::regex expr_date_rx{
|
|
"^ *# *File +expires +on +[0-8]{1,2} "
|
|
"+(January|February|March|April|May|June|July|August|September|October|November|December) +[0-9]{4} *$"};
|
|
|
|
const std::regex data_rx{"^ *[0-9]{5,}(\\.?[0-9]+) +[0-9]{1,2} +[0-9]{1,2} +[0-9]{4} +[0-9]{1,} *$"};
|
|
|
|
std::istringstream is;
|
|
double mjd;
|
|
unsigned day, month;
|
|
int year;
|
|
double tai_utc;
|
|
|
|
|
|
leapsecond_db_t db;
|
|
|
|
for (std::string line; std::getline(stream, line);) {
|
|
auto sv = utils::trimSpaces(line, utils::TrimType::TRIM_LEFT);
|
|
|
|
if (sv.size()) {
|
|
if (sv[0] == comment_sym) { // comment string
|
|
if (std::regex_match(line, expr_date_rx)) {
|
|
auto pos = line.find("on");
|
|
sv = utils::trimSpaces(std::string_view{line.begin() + pos + 2, line.end()},
|
|
utils::TrimType::TRIM_LEFT);
|
|
is.str({sv.begin(), sv.end()});
|
|
is >> std::chrono::parse("%d %B %Y", db.expireDate);
|
|
is.clear();
|
|
}
|
|
continue;
|
|
}
|
|
} else {
|
|
continue;
|
|
}
|
|
|
|
if (std::regex_match(line, data_rx)) {
|
|
is.str(line);
|
|
is >> mjd >> day >> month >> year >> tai_utc;
|
|
db.db.emplace_back(mjd, day, month, year, tai_utc);
|
|
is.clear();
|
|
continue;
|
|
}
|
|
}
|
|
|
|
if (db.db.empty()) {
|
|
db.state = IERS_DB_STATE_PARSE_ERROR;
|
|
} else {
|
|
db.state = IERS_DB_STATE_OK;
|
|
}
|
|
|
|
return db;
|
|
}
|
|
|
|
|
|
|
|
} // namespace mcc::astro
|