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-12 17:43:17 +03:00
parent a5bdf87a9b
commit 9e70ace4b7
4 changed files with 156 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files

131
cxx/mount_astrom.h
View File

@@ -19,6 +19,10 @@ 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> ||
@@ -31,13 +35,17 @@ namespace mcc::astro
{
// modified Julian date (based on ERFA eraCal2jd)
// template <traits::mcc_scalar_or_simd_c T>
template <traits::mcc_real_scalar_or_real_range_c ResT>
static int mcc_julday(const std::chrono::system_clock::time_point& start_time,
const std::chrono::system_clock::duration& step,
std::vector<double, xsimd::default_allocator<double>>& mjd)
ResT& mjd,
const std::chrono::system_clock::duration& step = std::chrono::milliseconds(100))
{
if (mjd.empty()) {
return -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;
@@ -74,27 +82,52 @@ static int mcc_julday(const std::chrono::system_clock::time_point& start_time,
double mjd_float = static_cast<double>(mjd_int) +
std::chrono::duration_cast<std::chrono::nanoseconds>(start_time - dd).count() * nanosec;
double d_step = std::chrono::duration_cast<std::chrono::nanoseconds>(step).count() * nanosec;
constexpr size_t reg_size = xsimd::batch<double>::size;
size_t vec_size = mjd.size() - mjd.size() % reg_size;
if constexpr (std::ranges::range<ResT>) {
double d_step = std::chrono::duration_cast<std::chrono::nanoseconds>(step).count() * nanosec;
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>{});
size_t i = 0;
// vectorized part
size_t i = 0;
for (; i < vec_size; i += vec_size) {
res_reg += step_reg;
#ifdef VEC_XSIMD
constexpr size_t reg_size = xsimd::batch<double>::size;
size_t vec_size = mjd_size - mjd_size % reg_size;
res_reg.store_aligned(mjd.data() + i);
}
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>{});
// scalar part
for (size_t j = i; j < mjd.size(); ++j) {
mjd[j] = mjd_float + j * d_step;
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;
}
@@ -174,7 +207,29 @@ struct leapsecond_db_elem_t {
double tai_utc; // TAI-UTC in seconds
};
typedef std::vector<leapsecond_db_elem_t> leapsecond_db_t;
// typedef std::vector<leapsecond_db_elem_t> leapsecond_db_t;
struct leapsecond_db_t {
std::chrono::system_clock::time_point expireDate;
std::vector<leapsecond_db_elem_t> db;
};
// init to some known state
static leapsecond_db_t CURRENT_LEAPSECONDS_DB;
// = {
// {41317.0, 1, 1, 1972, 10}, {41499.0, 1, 7, 1972, 11}, {41683.0, 1, 1, 1973, 12},
// {42048.0, 1, 1, 1974, 13}, {42413.0, 1, 1, 1975, 14}, {42778.0, 1, 1, 1976, 15},
// {43144.0, 1, 1, 1977, 16}, {43509.0, 1, 1, 1978, 17}, {43874.0, 1, 1, 1979, 18},
// {44239.0, 1, 1, 1980, 19}, {44786.0, 1, 7, 1981, 20}, {45151.0, 1, 7, 1982, 21},
// {45516.0, 1, 7, 1983, 22}, {46247.0, 1, 7, 1985, 23}, {47161.0, 1, 1, 1988, 24},
// {47892.0, 1, 1, 1990, 25}, {48257.0, 1, 1, 1991, 26}, {48804.0, 1, 7, 1992, 27},
// {49169.0, 1, 7, 1993, 28}, {49534.0, 1, 7, 1994, 29}, {50083.0, 1, 1, 1996, 30},
// {50630.0, 1, 7, 1997, 31}, {51179.0, 1, 1, 1999, 32}, {53736.0, 1, 1, 2006, 33},
// {54832.0, 1, 1, 2009, 34}, {56109.0, 1, 7, 2012, 35}, {57204.0, 1, 7, 2015, 36},
// {57754.0, 1, 1, 2017, 37}
// };
struct earth_orient_db_elem_t {
@@ -186,27 +241,11 @@ struct earth_orient_db_elem_t {
};
// typedef std::vector<earth_orient_db_elem_t> earth_orient_db_t;
struct earth_orient_db_t {
std::chrono::system_clock::time_point bulletinDate{};
std::vector<earth_orient_db_elem_t> db{};
};
// init to some known state
static leapsecond_db_t CURRENT_LEAPSECONDS_DB = {
{41317.0, 1, 1, 1972, 10}, {41499.0, 1, 7, 1972, 11}, {41683.0, 1, 1, 1973, 12},
{42048.0, 1, 1, 1974, 13}, {42413.0, 1, 1, 1975, 14}, {42778.0, 1, 1, 1976, 15},
{43144.0, 1, 1, 1977, 16}, {43509.0, 1, 1, 1978, 17}, {43874.0, 1, 1, 1979, 18},
{44239.0, 1, 1, 1980, 19}, {44786.0, 1, 7, 1981, 20}, {45151.0, 1, 7, 1982, 21},
{45516.0, 1, 7, 1983, 22}, {46247.0, 1, 7, 1985, 23}, {47161.0, 1, 1, 1988, 24},
{47892.0, 1, 1, 1990, 25}, {48257.0, 1, 1, 1991, 26}, {48804.0, 1, 7, 1992, 27},
{49169.0, 1, 7, 1993, 28}, {49534.0, 1, 7, 1994, 29}, {50083.0, 1, 1, 1996, 30},
{50630.0, 1, 7, 1997, 31}, {51179.0, 1, 1, 1999, 32}, {53736.0, 1, 1, 2006, 33},
{54832.0, 1, 1, 2009, 34}, {56109.0, 1, 7, 2012, 35}, {57204.0, 1, 7, 2015, 36},
{57754.0, 1, 1, 2017, 37}
};
static earth_orient_db_t CURRENT_EARTH_ORIENT_DB;
@@ -300,13 +339,20 @@ static bool mcc_parse_leapsecs(std::derived_from<std::basic_istream<char>> auto&
leapsecond_db_t db;
// leapsecond_db_elem_t db_elem;
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 {
@@ -314,14 +360,15 @@ static bool mcc_parse_leapsecs(std::derived_from<std::basic_istream<char>> auto&
}
if (std::regex_match(line, data_rx)) {
is.str(std::move(line));
is.str(line);
is >> mjd >> day >> month >> year >> tai_utc;
db.emplace_back(mjd, day, month, year, tai_utc);
db.db.emplace_back(mjd, day, month, year, tai_utc);
is.clear();
continue;
}
}
if (db.empty()) {
if (db.db.empty()) {
return false;
}