diff --git a/cxx/mount_astrom.h b/cxx/mount_astrom.h index 8e9eac5..f3098e3 100644 --- a/cxx/mount_astrom.h +++ b/cxx/mount_astrom.h @@ -55,6 +55,7 @@ template >) { size_t mjd_size = 0; if constexpr (std::ranges::range) { @@ -101,7 +102,7 @@ static int mcc_julday(traits::mcc_systime_c auto const& start_time, constexpr size_t reg_size = xsimd::batch::size; size_t vec_size = mjd_size - mjd_size % reg_size; - xsimd::batch res_reg{mjd_float}; + xsimd::batch res_reg(mjd_float); xsimd::batch step_reg = [d_step](std::index_sequence) { return xsimd::batch{(Is * d_step)...}; }(std::make_index_sequence{}); @@ -111,7 +112,7 @@ static int mcc_julday(traits::mcc_systime_c auto const& start_time, auto ptr = mjd.begin(); // vectorized part - for (; i < vec_size; i += vec_size) { + for (; i < vec_size; i += reg_size) { res_reg += step_reg; if constexpr (std::ranges::contiguous_range) { @@ -143,6 +144,18 @@ static int mcc_julday(traits::mcc_systime_c auto const& start_time, return 0; } +template +static int mcc_julday(traits::mcc_systime_c auto const& start_time, + ResT& mjd, + const DT& step = std::chrono::milliseconds(100)) + requires std::is_pointer_v> +{ + auto sp = std::span(mjd); + int ret = mcc_julday(start_time, sp, step); + return ret; +} + + /* * angles are in degrees or sexagimal string form * RA and DEC are apparent! diff --git a/cxx/tests/astrom_test.cpp b/cxx/tests/astrom_test.cpp index 53b0c37..ea87971 100644 --- a/cxx/tests/astrom_test.cpp +++ b/cxx/tests/astrom_test.cpp @@ -68,12 +68,39 @@ int main(int argc, char* argv[]) // const double rr = 180.0 / std::numbers::pi * 60.0; // std::cout << "A(arcmin) = " << A * rr << "; B(arcmin) = " << B * rr << "\n"; + float alt_lim = 10.0; + std::string_view ra_str = "06:30:00.0", dec_str = "00:00:00.0"; + std::cout << "\n\nTimes to object (RA = " << ra_str << ", DEC = " << dec_str << ") sets to given altitude (" + << alt_lim << " degrees):\n"; using namespace std::chrono_literals; - double stm = mcc::astrom::mcc_time_to_alt_limit(10.0, "06:30:00.0", "00:00:00.0", 43.646711, 41.440732, + double stm = mcc::astrom::mcc_time_to_alt_limit(alt_lim, ra_str, dec_str, 43.646711, 41.440732, std::chrono::system_clock::now(), 0.041s, 32.184s, 37.0s); // auto stm_d = mcc::astrom::mcc_chrono_radians{stm}; std::cout << "STM: " << stm * 12.0 / std::numbers::pi * 60.0 << " minutes\n"; + + std::cout << "\n\n\n"; + const size_t N = 1000; + double mjds[N]; + now = std::chrono::system_clock::now(); + st = mcc::astrom::mcc_julday(now, mjds, std::chrono::milliseconds(100)); + std::cout << std::format("comp time for {} MJDs = {}\n", N, std::chrono::system_clock::now() - now); + std::cout << std::format("MIN MJD: {:.16f}; \nMAX MJD: {:.16f}\n", mjds[0], mjds[N - 1]); + + double elong = 43.646711 * std::numbers::pi / 180.0; + double phi = 41.440732 * std::numbers::pi / 180.0; + double aob, zob, hob, dob, rob, eo; + + now = std::chrono::system_clock::now(); + + for (auto& el : mjds) { + st = erfa::eraAtco13(1.343523, 0.32352345, 0.0, 0.0, 0.0, 0.0, ERFA_DJM0, el, 0.041, elong, phi, 2100.0, 0.0, + 0.0, 700.0, 10.0, 0.8, 0.5, &aob, &zob, &hob, &dob, &rob, &eo); + } + std::cout << std::format( + "comp time for {} coordinate transf = {}\n", N, + std::chrono::duration_cast(std::chrono::system_clock::now() - now)); + return ecode; }