102 lines
2.9 KiB
C++
102 lines
2.9 KiB
C++
#pragma once
|
|
|
|
/*********************************
|
|
* MOUNT CONTROL COMPONENTS *
|
|
* *
|
|
* astrometry functions *
|
|
*********************************/
|
|
|
|
#include <chrono>
|
|
#include <xsimd/xsimd.hpp>
|
|
|
|
|
|
namespace mcc::traits
|
|
{
|
|
|
|
template <typename T>
|
|
concept mcc_scalar_or_simd_c = xsimd::is_batch<T>::value || std::is_arithmetic_v<T>;
|
|
|
|
} // namespace mcc::traits
|
|
|
|
|
|
namespace mcc::astro
|
|
{
|
|
|
|
// modified Julian date (based on ERFA eraCal2jd)
|
|
// template <traits::mcc_scalar_or_simd_c T>
|
|
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)
|
|
{
|
|
if (mjd.empty()) {
|
|
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;
|
|
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;
|
|
|
|
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>{});
|
|
|
|
// vectorized part
|
|
size_t i = 0;
|
|
for (; i < vec_size; i += vec_size) {
|
|
res_reg += step_reg;
|
|
|
|
res_reg.store_aligned(mjd.data() + i);
|
|
}
|
|
|
|
// scalar part
|
|
for (size_t j = i; j < mjd.size(); ++j) {
|
|
mjd[j] = mjd_float + j * d_step;
|
|
}
|
|
|
|
|
|
return 0;
|
|
}
|
|
|
|
template <traits::mcc_scalar_or_simd_c T, typename CT, typename DT>
|
|
T mcc_time_to_alt_limit(const T& alt_limit, const T& RA, const T& DEC, const std::chrono::time_point<CT, DT>& now)
|
|
{
|
|
}
|
|
|
|
} // namespace mcc::astro
|