From 0182844af6e17159c1b0d50fe85a19f8b60a22b1 Mon Sep 17 00:00:00 2001 From: "Timur A. Fatkhullin" Date: Tue, 24 Mar 2026 02:43:17 +0300 Subject: [PATCH] ... --- CMakeLists.txt | 3 + include/mcc/mcc_angle.h | 6 +- include/mcc/mcc_pcm.h | 2 +- include/mcc/mcc_pcm_construct.h | 24 ++++--- tests/mcc_pcm_test.cpp | 113 ++++++++++++++++++++++++++++++++ 5 files changed, 132 insertions(+), 16 deletions(-) create mode 100644 tests/mcc_pcm_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index ecc8063..9bdf70e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -335,6 +335,9 @@ if(BUILD_TESTS) add_executable(mcc_fitpack_test tests/mcc_fitpack_test.cpp) target_link_libraries(mcc_fitpack_test PRIVATE ${PROJECT_NAME}) + + add_executable(mcc_pcm_test tests/mcc_pcm_test.cpp) + target_link_libraries(mcc_pcm_test PRIVATE ${PROJECT_NAME}) else() # This is just a stub to allow access to the path and library settings for the ${PROJECT_NAME} target during development add_executable(just_stub EXCLUDE_FROM_ALL main.cpp) diff --git a/include/mcc/mcc_angle.h b/include/mcc/mcc_angle.h index 0fb2cf8..a15ac1f 100644 --- a/include/mcc/mcc_angle.h +++ b/include/mcc/mcc_angle.h @@ -179,7 +179,8 @@ public: template MccAngle& normalize() { - _angleInRads = std::fmod(_angleInRads, std::numbers::pi * 2.0); + // _angleInRads = std::fmod(_angleInRads, std::numbers::pi * 2.0); + _angleInRads = std::fmod(_angleInRads, MCC_TWO_PI); if constexpr (KIND == NORM_KIND_0_360) { if (_angleInRads < 0.0) { @@ -195,7 +196,8 @@ public: } else if constexpr (KIND == NORM_KIND_180_180) { if (_angleInRads > std::numbers::pi) { // _angleInRads = 2.0 * std::numbers::pi - _angleInRads; - _angleInRads = MCC_TWO_PI - _angleInRads; + // _angleInRads = MCC_TWO_PI - _angleInRads; + _angleInRads = _angleInRads - MCC_TWO_PI; } else if (_angleInRads < -std::numbers::pi) { // _angleInRads += 2.0 * std::numbers::pi; _angleInRads += MCC_TWO_PI; diff --git a/include/mcc/mcc_pcm.h b/include/mcc/mcc_pcm.h index 8ee6ef1..27400a3 100644 --- a/include/mcc/mcc_pcm.h +++ b/include/mcc/mcc_pcm.h @@ -148,8 +148,8 @@ public: // alt-azimuthal coeff_t tubeFlexure; - coeff_t forkFlexure; coeff_t DECaxisFlexure; // declination axis flexure + coeff_t forkFlexure; }; #ifdef USE_BSPLINE_PCM diff --git a/include/mcc/mcc_pcm_construct.h b/include/mcc/mcc_pcm_construct.h index e70d756..7eaef56 100644 --- a/include/mcc/mcc_pcm_construct.h +++ b/include/mcc/mcc_pcm_construct.h @@ -138,7 +138,7 @@ public: }; MccError addPoint(mcc_skypoint_c auto target_coords, mcc_coord_pair_c auto const& hw_counts) - requires(decltype(hw_counts)::pairKind == MccCoordPairKind::COORDS_KIND_XY) + requires(std::decay_t::pairKind == MccCoordPairKind::COORDS_KIND_XY) { auto err = target_coords.to(ref_coordpair_t::pairKind, hw_counts.epoch()); if (err) { @@ -232,7 +232,7 @@ public: min_data_size = 8; } - if (_table.size() < min_data_size) { + if (numberOfPoints() < min_data_size) { result.error = MccDefaultPCMConstructorErrorCode::ERROR_NOT_ENOUGH_DATA; return result; } @@ -240,7 +240,7 @@ public: // robust linear regression with Tukey's loss function result = robustLinearRegress(pcm_data, comp_params); } else { - if (_table.size() < min_data_size) { + if (numberOfPoints() < min_data_size) { result.error = MccDefaultPCMConstructorErrorCode::ERROR_NOT_ENOUGH_DATA; return result; } @@ -256,8 +256,8 @@ public: std::vector xres = _table.colon_res, yres = _table.colat_res; for (size_t i = 0; i < numberOfPoints(); ++i) { - _table.colon_res = _table.target_colon[i] - result.model_colon; - _table.colat_res = _table.target_colat[i] - result.model_colat; + _table.colon_res[i] = _table.target_colon[i] - result.model_colon[i]; + _table.colat_res[i] = _table.target_colat[i] - result.model_colat[i]; } auto r = bsplineFitting(pcm_data); @@ -266,11 +266,11 @@ public: if (!r.error) { for (size_t i = 0; i < numberOfPoints(); ++i) { - result.model_colon += r.model_colon; - result.model_colat += r.model_colat; + result.model_colon[i] += r.model_colon[i]; + result.model_colat[i] += r.model_colat[i]; - result.colon_res = _table.target_colon[i] - result.model_colon; - result.colat_res = _table.target_colat[i] - result.model_colat; + result.colon_res[i] = _table.target_colon[i] - result.model_colon[i]; + result.colat_res[i] = _table.target_colat[i] - result.model_colat[i]; } } @@ -424,11 +424,9 @@ protected: 1.0E-8; }; - compute_result_t result{.pcm_type = pcm_data.type, - .error = MccDefaultPCMConstructorErrorCode::ERROR_OK, - .max_iter = comp_params.max_iter}; + compute_result_t result{.pcm_type = pcm_data.type, .error = MccDefaultPCMConstructorErrorCode::ERROR_OK}; - const size_t n_coeffs = pcm_data.type == MOUNT_TYPE == MccMountType::FORK_TYPE ? 9 : 8; + const size_t n_coeffs = MOUNT_TYPE == MccMountType::FORK_TYPE ? 9 : 8; Eigen::MatrixXd X(2 * numberOfPoints(), n_coeffs); Eigen::VectorXd y(2 * numberOfPoints()), beta(n_coeffs), new_beta(n_coeffs); diff --git a/tests/mcc_pcm_test.cpp b/tests/mcc_pcm_test.cpp new file mode 100644 index 0000000..1bb2f15 --- /dev/null +++ b/tests/mcc_pcm_test.cpp @@ -0,0 +1,113 @@ +#include +#include + + +#include + +struct pcm_res_t { + double pcmX, pcmY; +}; + +static constexpr mcc::MccMountType MOUNT_TYPE{mcc::MccMountType::FORK_TYPE}; + +int main() +{ + using pcm_t = mcc::impl::MccDefaultPCM; + using pcm_const_t = mcc::impl::MccPCMConstructor; + + pcm_t pcm; + pcm_const_t pcm_const; + + pcm_t::pcm_data_t pcm_data{.type = mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY, + .siteLatitude = 43.6466666667_degs, + .geomCoefficients = {.zeroPointX = 2.434534_degs, + .zeroPointY = -3.382549_degs, + .collimationErr = 10.534_arcsecs, + .nonperpendErr = 21.86834_arcsecs, + .misalignErr1 = 7.345768_arcsecs, + .misalignErr2 = -12.87934_arcsecs, + .tubeFlexure = 5.9083_arcsecs, + .DECaxisFlexure = 3.48379812_arcsecs, + .forkFlexure = 9.342354_arcsecs}}; + + + size_t N = 50; + + std::vector encX(N), encY(N), ha(N), dec(N), pcmX(N), pcmY(N); + std::vector resX(N), resY(N); + pcm_res_t pcm_res; + + + + auto comp_hadec = [&]() { + mcc::impl::MccSkyPoint hadec; + for (size_t i = 0; i < N; ++i) { + auto err = pcm.computePCM(mcc::impl::MccGenXY{(double)encX[i], (double)encY[i]}, &pcm_res, &hadec); + if (err) { + std::println("Error in PCM computing: {}", err.message()); + return 1; + } + + pcmX[i] = pcm_res.pcmX; + pcmY[i] = pcm_res.pcmY; + + ha[i] = hadec.co_lon(); + dec[i] = hadec.co_lat(); + } + return 0; + }; + + pcm.setPCMData(pcm_data); + + double start_encX = 1.2423_degs, stop_encX = 354.896124_degs; + double start_encY = -32.2423_degs, stop_encY = 89.896124_degs; + double sigma = 24.3248792_arcsecs; + + double x_step = (stop_encX - start_encX) / (N - 1), y_step = (stop_encY - start_encY) / (N - 1); + + for (size_t i = 0; i < N; ++i) { + encX[i] = start_encX + i * x_step; + encY[i] = start_encY + i * y_step; + } + + comp_hadec(); + + std::println("X Y HA DEC (pcmX pcmY):"); + for (size_t i = 0; i < N; ++i) { + std::println("{:12.7f} {:12.7f} {:12.7f} {:12.7f} ({:10.7f} {:10.7f})", encX[i].degrees(), encY[i].degrees(), + ha[i].degrees(), dec[i].degrees(), pcmX[i].degrees(), pcmY[i].degrees()); + } + + + // add 'noise' + std::random_device rd{}; + std::mt19937 gen{rd()}; + std::normal_distribution dist{0.0, sigma}; + + mcc::impl::MccSkyPoint hadec; + + for (size_t i = 0; i < N; ++i) { + ha[i] += dist(gen); + dec[i] += dist(gen); + + resX[i] = ha[i] - encX[i]; + resY[i] = dec[i] - encY[i]; + + hadec.from(mcc::impl::MccSkyHADEC_OBS{(double)ha[i], (double)dec[i]}); + + pcm_const.addPoint(hadec, mcc::impl::MccGenXY{(double)encX[i], (double)encY[i]}); + } + + pcm_t::pcm_data_t fit_pcm_data{.type = pcm_data.type, .siteLatitude = pcm_data.siteLatitude}; + + auto r = pcm_const.computeModel(fit_pcm_data); + if (r.error) { + std::println("error: {}", r.error.message()); + return 1; + } + + std::println("FITTED COEFFS:"); + std::println("X0 = {}, Y0 = ", fit_pcm_data.geomCoefficients.zeroPointX, fit_pcm_data.geomCoefficients.zeroPointY); + + return 0; +} \ No newline at end of file