eddy/mountcontrol
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

...

Browse Source
This commit is contained in:
Timur A. Fatkhullin
2025-03-19 17:31:03 +03:00
parent cf37281e0f
commit 2a838830c8
4 changed files with 168 additions and 342 deletions
Download Patch File Download Diff File Expand all files Collapse all files

339
cxx/mount_astrom.h
View File

@@ -33,9 +33,10 @@ concept mcc_real_or_char_range =
} // namespace mcc::traits
namespace mcc::astro
namespace mcc::astrom
{
// 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,
@@ -64,11 +65,15 @@ static int mcc_julday(const std::chrono::system_clock::time_point& start_time,
return -2;
}
int64_t im = (unsigned)ymd.month();
int64_t id = (unsigned)ymd.day();
int64_t iy = (int)ymd.year();
// my = (im - 14) / 12;
// iypmy = (long) (iy + my);
int64_t my = -(14 - (unsigned)ymd.month()) / 12;
int64_t iypmy = (int)ymd.year() + my;
int64_t my = (im - 14LL) / 12LL;
int64_t iypmy = iy + my;
// (1461L * (iypmy + 4800L)) / 4L
// + (367L * (long) (im - 2 - 12 * my)) / 12L
@@ -76,17 +81,14 @@ static int mcc_julday(const std::chrono::system_clock::time_point& start_time,
// + (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;
int64_t mjd_int = (1461LL * (iypmy + 4800LL)) / 4LL + (367LL * (im - 2LL - 12LL * my)) / 12LL -
(3LL * ((iypmy + 4900LL) / 100LL)) / 4LL + id - 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;
using fd_t = std::chrono::duration<double, std::ratio<86400>>; // fraction of day
double mjd_float = static_cast<double>(mjd_int) + std::chrono::duration_cast<fd_t>(start_time - dd).count();
if constexpr (std::ranges::range<ResT>) {
double d_step = std::chrono::duration_cast<std::chrono::nanoseconds>(step).count() * nanosec;
double d_step = std::chrono::duration_cast<fd_t>(step).count();
size_t i = 0;
@@ -211,200 +213,20 @@ double mcc_time_to_alt_limit(const AT& alt_limit,
}
/* IERS bulletins
/* Classes to represent 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
namespace mcc::astrom
{
class MccLeapSeconds
class MccLeapSeconds final
{
public:
typedef std::chrono::system_clock::time_point time_point_t;
MccLeapSeconds()
{
// create default values
@@ -416,7 +238,7 @@ public:
~MccLeapSeconds() = default;
std::chrono::system_clock::time_point expireDate() const
time_point_t expireDate() const
{
return _expireDate;
}
@@ -494,7 +316,7 @@ public:
}
std::optional<double> operator[](const std::chrono::system_clock::time_point& tp) const
std::optional<double> operator[](const time_point_t& tp) const
{
if (tp > _expireDate) { // ???????!!!!!!!!!!!
return std::nullopt;
@@ -516,7 +338,7 @@ public:
{
double e_mjd;
astro::mcc_julday(_expireDate, e_mjd);
astrom::mcc_julday(_expireDate, e_mjd);
if (mjd > e_mjd) { // ???????!!!!!!!!!!!
return std::nullopt;
// return _db.back().tai_utc;
@@ -531,6 +353,16 @@ public:
return std::nullopt;
}
void dump(std::derived_from<std::basic_ostream<char>> auto& stream) const
{
stream << std::format("Leap seconds database expire date: {}", _expireDate) << '\n';
for (auto const& el : _db) {
stream << std::format("{} {} {}", el.mjd, el.ymd, el.tai_utc) << '\n';
}
}
private:
inline static const std::regex expr_date_rx{
"^ *# *File +expires +on +[0-8]{1,2} "
@@ -538,7 +370,7 @@ private:
inline static 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::chrono::system_clock::time_point _expireDate{};
time_point_t _expireDate{};
struct leapsecond_db_elem_t {
double mjd;
@@ -551,37 +383,74 @@ private:
class MccIersBulletinA
class MccIersBulletinA final
{
public:
typedef std::chrono::system_clock::time_point time_point_t;
struct pole_pos_t {
double x, y;
};
struct date_range_t {
std::chrono::year_month_day begin;
std::chrono::year_month_day end;
};
struct date_range_mjd_t {
double begin;
double end;
};
MccIersBulletinA()
{
// create pre-defined (default-state) database
std::istringstream ist(defaults::MCC_DEFAULT_IERS_BULLETIN_A_FILE);
load(ist);
}
~MccIersBulletinA() = default;
std::chrono::system_clock::time_point bulletinDate() const
{
return _date;
}
date_range_t dateRange() const
{
return {_db.front().ymd, _db.back().ymd};
}
date_range_mjd_t dateRangeMJD() const
{
return {_db.front().mjd, _db.back().mjd};
}
double TT_TAI() const
{
return _tt_tai;
}
// DUT1 = UT1 - UTC
std::optional<double> DUT1(const std::chrono::system_clock::time_point& tp) const
std::optional<double> DUT1(const time_point_t& tp) const
{
std::chrono::year_month_day ymd{std::chrono::floor<std::chrono::days>(tp)};
// use of the closest date
std::chrono::year_month_day ymd{std::chrono::round<std::chrono::days>(tp)};
for (auto const& el : _db | std::views::reverse) {
if (ymd >= el.ymd) {
if (ymd < _db.front().ymd) {
return std::nullopt;
}
if (ymd > _db.back().ymd) {
return std::nullopt;
}
auto el_prev = _db.front();
for (auto const& el : _db) {
if (ymd <= el.ymd) {
return el.dut1;
}
}
@@ -589,10 +458,20 @@ public:
return std::nullopt;
}
std::optional<double> DUT1(const double& mjd) const
std::optional<double> DUT1(double mjd) const
{
for (auto const& el : _db | std::views::reverse) {
if (mjd >= el.mjd) {
mjd = std::round(mjd); // round to closest integer MJD
if (mjd < _db.front().mjd) {
return std::nullopt;
}
if (mjd > _db.back().mjd) {
return std::nullopt;
}
for (auto const& el : _db) {
if (mjd <= el.mjd) {
return el.dut1;
}
}
@@ -600,12 +479,20 @@ public:
return std::nullopt;
}
std::optional<pole_pos_t> polePos(const std::chrono::system_clock::time_point& tp) const
std::optional<pole_pos_t> polarCoords(const time_point_t& tp) const
{
std::chrono::year_month_day ymd{std::chrono::floor<std::chrono::days>(tp)};
std::chrono::year_month_day ymd{std::chrono::round<std::chrono::days>(tp)};
for (auto const& el : _db | std::views::reverse) {
if (ymd >= el.ymd) {
if (ymd < _db.front().ymd) {
return std::nullopt;
}
if (ymd > _db.back().ymd) {
return std::nullopt;
}
for (auto const& el : _db) {
if (ymd <= el.ymd) {
return pole_pos_t{el.x, el.y};
}
}
@@ -613,10 +500,20 @@ public:
return std::nullopt;
}
std::optional<pole_pos_t> polePos(const double& mjd) const
std::optional<pole_pos_t> polarCoords(double mjd) const
{
for (auto const& el : _db | std::views::reverse) {
if (mjd >= el.mjd) {
mjd = std::round(mjd); // round to closest integer MJD
if (mjd < _db.front().mjd) {
return std::nullopt;
}
if (mjd > _db.back().mjd) {
return std::nullopt;
}
for (auto const& el : _db) {
if (mjd <= el.mjd) {
return pole_pos_t{el.x, el.y};
}
}
@@ -714,6 +611,16 @@ public:
return ok;
}
void dump(std::derived_from<std::basic_ostream<char>> auto& stream) const
{
stream << std::format("Bulletin A issue date: {}", _date) << '\n';
stream << std::format("TT-TAI: {}", _tt_tai) << '\n';
for (auto const& el : _db) {
stream << std::format("{} {} {:6.4f} {:6.4f} {:7.5f}", el.mjd, el.ymd, el.x, el.y, el.dut1) << '\n';
}
}
private:
inline static const std::regex bull_date_rx{
"^ *[0-9]{1,2} +(January|February|March|April|May|June|July|August|September|October|November|December) "
@@ -727,7 +634,7 @@ private:
"^ *[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]+ *$"};
std::chrono::system_clock::time_point _date;
time_point_t _date;
double _tt_tai;
struct earth_orient_db_elem_t {