diff --git a/include/mcc/mcc_pcm.h b/include/mcc/mcc_pcm.h index 671251a..a13606c 100644 --- a/include/mcc/mcc_pcm.h +++ b/include/mcc/mcc_pcm.h @@ -485,6 +485,17 @@ private: } } + + if constexpr (mcc_is_equatorial_mount) { // HA - DEC system (needs normalization) + if (!inverse) { + x = MccAngleHA_OBS(x); + y = MccAngleDEC_OBS(y); + } + + // in the case of inverse PCM it is assumed that angles are already normalized (inputs are celestial + // coordiinates) + } + pcm_geom_coeffs_t* geom_coeffs = &_pcmData.geomCoefficients; #ifdef USE_BSPLINE_PCM diff --git a/include/mcc/mcc_pcm_fit.h b/include/mcc/mcc_pcm_fit.h index b5a894d..073374d 100644 --- a/include/mcc/mcc_pcm_fit.h +++ b/include/mcc/mcc_pcm_fit.h @@ -112,9 +112,9 @@ public: size_t final_iter{}; // number of final iteration - std::vector model_colon{}, model_colat{}; // fitting model values - std::vector colon_res{}, colat_res{}; // target - model - std::vector colon_weight{}, colat_weight; // Tukey's weights + std::vector model_colonRES{}, model_colatRES{}; // fitting model values: fit(target-hw) + std::vector colon_err{}, colat_err{}; // fitting residuals: (target-hw) - fit + std::vector colon_weight{}, colat_weight; // Tukey's weights #ifdef USE_BSPLINE_PCM std::array bspline_fit_err{ @@ -122,8 +122,8 @@ public: // 0 - pcmX(X,Y), 1 - pcmY(X,Y), 2 - revPcmX(COLON, COLAT), 3 - revPcmY(COLON, COLAT) // quantities below are computed only fo pcm_type == MccDefaultPCMType::PCM_TYPE_BSPLINE - std::vector inv_model_colon{}, inv_model_colat{}; // fitted inverse model values - std::vector inv_colon_res{}, inv_colat_res{}; // encoder - model + std::vector inv_model_colonRES{}, inv_model_colatRES{}; // fitted inverse model values + std::vector inv_colon_err{}, inv_colat_err{}; // encoder - model #endif }; @@ -134,6 +134,15 @@ public: }; + struct pcm_table_elem_t { + ref_coordpair_t target{}; // celestial + MccGenXY hw{}; // hardware (encoder) + MccGenXY res{}; // residuals: celestial - hardware + }; + + typedef std::vector pcm_table_t; + + size_t numberOfPoints() const { return _targetCOLON.size(); @@ -206,21 +215,31 @@ public: } - // something like: - // auto [tag, hw, err] = getPoint(7); + template + MccError addPoint(COLON_T const& colon, + COLAT_T const& colat, + HW_X_T const& x, + HW_Y_T const& y, + mcc_coord_epoch_c auto const& epoch) + requires(std::is_arithmetic_v && std::is_arithmetic_v && std::is_arithmetic_v && + std::is_arithmetic_v) + { + ref_coordpair_t cp{typename ref_coordpair_t::x_t{colon}, typename ref_coordpair_t::y_t{colat}, epoch}; + MccGenXY hw{x, y, epoch}; + + return addPoint(cp, hw); + } + + template - std::tuple getPoint(size_t idx) + MccError getPoint(size_t idx, std::tuple& point) requires(std::same_as && std::same_as && HW_T::pairKind == MccCoordPairKind::COORDS_KIND_XY) { - std::tuple point; - if (idx > numberOfPoints()) { - std::get<2>(point) = MccPCMFitterErrorCode::ERROR_INVALID_INDEX; + return MccPCMFitterErrorCode::ERROR_INVALID_INDEX; } else { - std::get<2>(point) = MccPCMFitterErrorCode::ERROR_OK; - std::get<0>(point).setX(_targetCOLON[idx]); std::get<0>(point).setY(_targetCOLAT[idx]); std::get<0>(point).setEpoch(_epoch[idx]); @@ -230,7 +249,103 @@ public: std::get<1>(point).setEpoch(_epoch[idx]); } - return point; + return MccPCMFitterErrorCode::ERROR_OK; + } + + + pcm_table_t getPCMTable() const + { + pcm_table_t tab; + + if (numberOfPoints()) { + for (size_t i = 0; i < numberOfPoints(); ++i) { + tab.emplace_back(pcm_table_elem_t{.target{typename ref_coordpair_t::x_t{_targetCOLON[i]}, + typename ref_coordpair_t::y_t{_targetCOLAT[i]}, _epoch[i]}, + .hw{MccAngleX{_hwX[i]}, MccAngleY{_hwY[i]}, _epoch[i]}, + .res{MccAngleX{_colonRES[i]}, MccAngleY{_colatRES[i]}, _epoch[i]}}); + } + } + + return tab; + } + + + // + // for B-splines interior knots along axes must be given in 'pcm_data' + // NOTE: the size of the interior knots array must be at least 2 as + // it are interpretated as border knots and final full knots set is: + // knots = [input_knots[0], input_knots[0], input_knots[0], input_knots[0], input_knots[1], input_knots[2], + // ..., input_knots[N-1], input_knots[N-1], input_knots[N-1], input_knots[N-1]], where N = input_knots.size() + // + // WARNING: the input knots for inverse B-spline are ignored so the direct and inverse B-spline coefficients are + // calculated on the same mesh! + compute_result_t computeModel(MccDefaultPCM::pcm_data_t& pcm_data, + compute_params_t const& comp_params = {}) + { + compute_result_t result{.pcm_type = pcm_data.type, .error = MccPCMFitterErrorCode::ERROR_OK}; + + size_t min_data_size = 2; // 2 is for BSPLINE + + if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY +#ifdef USE_BSPLINE_PCM + || pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE +#endif + ) { + if constexpr (MOUNT_TYPE == MccMountType::FORK_TYPE) { + min_data_size = 9; + } else { + min_data_size = 8; + } + + if (numberOfPoints() < min_data_size) { + result.error = MccPCMFitterErrorCode::ERROR_NOT_ENOUGH_DATA; + return result; + } + + // robust linear regression with Tukey's loss function + result = robustLinearRegress(pcm_data, comp_params); + } else { + if (numberOfPoints() < min_data_size) { + result.error = MccPCMFitterErrorCode::ERROR_NOT_ENOUGH_DATA; + return result; + } + } + + +#ifdef USE_BSPLINE_PCM + if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) { + return bsplineFitting(pcm_data); + } else if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE) { + // the fitting for geometrical coefficients is already done above so + // one must fit residuals by bivariate B-splines + + std::vector xres = _colonRES, yres = _colatRES; + for (size_t i = 0; i < numberOfPoints(); ++i) { + _colonRES[i] = _targetCOLON[i] - result.model_colonRES[i]; + _colatRES[i] = _targetCOLAT[i] - result.model_colatRES[i]; + } + + auto r = bsplineFitting(pcm_data); + result.error = r.error; + result.bspline_fit_err = r.bspline_fit_err; + + if (!r.error) { + for (size_t i = 0; i < numberOfPoints(); ++i) { + result.model_colonRES[i] += r.model_colonRES[i]; + result.model_colatRES[i] += r.model_colatRES[i]; + + result.colon_err[i] = _targetCOLON[i] - result.model_colonRES[i]; + result.colat_err[i] = _targetCOLAT[i] - result.model_colatRES[i]; + } + } + + // restore original residuals + _colonRES = xres; + _colatRES = yres; + } +#endif + + return result; } protected: @@ -242,21 +357,21 @@ protected: void computeModelResi(compute_result_t& result) { - result.colon_res.resize(numberOfPoints()); - result.colat_res.resize(numberOfPoints()); + result.colon_err.resize(numberOfPoints()); + result.colat_err.resize(numberOfPoints()); for (size_t i = 0; i < numberOfPoints(); ++i) { - result.colon_res[i] = _colonRES[i] - result.model_colon[i]; // = target - model - result.colat_res[i] = _colatRES[i] - result.model_colat[i]; // = target - model + result.colon_err[i] = _colonRES[i] - result.model_colonRES[i]; // = target - model + result.colat_err[i] = _colatRES[i] - result.model_colatRES[i]; // = target - model } #ifdef USE_BSPLINE_PCM if (result.pcm_type == MccDefaultPCMType::PCM_TYPE_BSPLINE) { - result.inv_colon_res.resize(numberOfPoints()); - result.inv_colat_res.resize(numberOfPoints()); + result.inv_colon_err.resize(numberOfPoints()); + result.inv_colat_err.resize(numberOfPoints()); for (size_t i = 0; i < numberOfPoints(); ++i) { - result.inv_colon_res[i] = _colonRES[i] - result.inv_model_colon[i]; // = hw - model - result.inv_colat_res[i] = _colatRES[i] - result.inv_model_colat[i]; // = hw - model + result.inv_colon_err[i] = _colonRES[i] - result.inv_model_colonRES[i]; // = hw - model + result.inv_colat_err[i] = _colatRES[i] - result.inv_model_colatRES[i]; // = hw - model } } #endif @@ -279,9 +394,25 @@ protected: // full knot vectors: [t_min, t_min, t_min, t_min, ..., t_1, t_2, ..., t_max, t_max, t_max, t_max] std::vector tx(pcm_data.bspline.knotsX.size() + 6), ty(pcm_data.bspline.knotsY.size() + 6); - std::vector&theta_hw = _hwY, theta_tag = _targetCOLAT, phi_hw = _hwX, phi_tag = _targetCOLON; + std::vector&theta_hw = _hwY, &theta_tag = _targetCOLAT, &phi_hw = _hwX, &phi_tag = _targetCOLON; std::vector theta, phi; + // full set of B-spline knots + tx[0] = tx[1] = tx[2] = pcm_data.bspline.knotsX[0]; + tx[tx.size() - 1] = tx[tx.size() - 2] = tx[tx.size() - 3] = + pcm_data.bspline.knotsX[pcm_data.bspline.knotsX.size() - 1]; + + ty[0] = ty[1] = ty[2] = pcm_data.bspline.knotsY[0]; + ty[ty.size() - 1] = ty[ty.size() - 2] = ty[ty.size() - 3] = + pcm_data.bspline.knotsY[pcm_data.bspline.knotsY.size() - 1]; + + for (size_t i = 0; i < pcm_data.bspline.knotsX.size(); ++i) { + tx[3 + i] = pcm_data.bspline.knotsX[i]; + } + for (size_t i = 0; i < pcm_data.bspline.knotsY.size(); ++i) { + ty[3 + i] = pcm_data.bspline.knotsY[i]; + } + /* WARNING: FITPACK B-spline on sphere: in the fitting routines the first angle argument is THETA - co-latitude @@ -303,21 +434,14 @@ protected: theta_hw[i] = _hwY[i] + MCC_HALF_PI; phi_hw[i] = _hwX[i] + std::numbers::pi; } - } - tx[0] = tx[1] = tx[2] = pcm_data.bspline.knotsX[0]; - tx[tx.size() - 1] = tx[tx.size() - 2] = tx[tx.size() - 3] = - pcm_data.bspline.knotsX[pcm_data.bspline.knotsX.size() - 1]; + for (size_t i = 0; i < tx.size(); ++i) { + tx[i] += std::numbers::pi; + } - ty[0] = ty[1] = ty[2] = pcm_data.bspline.knotsY[0] + MCC_HALF_PI; - ty[ty.size() - 1] = ty[ty.size() - 2] = ty[ty.size() - 3] = - pcm_data.bspline.knotsY[pcm_data.bspline.knotsY.size() - 1] + MCC_HALF_PI; - - for (size_t i = 0; i < pcm_data.bspline.knotsX.size(); ++i) { - tx[3 + i] = pcm_data.bspline.knotsX[i]; - } - for (size_t i = 0; i < pcm_data.bspline.knotsY.size(); ++i) { - ty[3 + i] = pcm_data.bspline.knotsY[i] + MCC_HALF_PI; + for (size_t i = 0; i < ty.size(); ++i) { + ty[i] += MCC_HALF_PI; + } } size_t Ncoeffs = (tx.size() - 4) * (ty.size() - 4); @@ -343,10 +467,10 @@ protected: bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.coeffsX, theta_hw, phi_hw, - result.model_colon); // get fitted residuals!!! + result.model_colonRES); // get fitted residuals!!! bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.coeffsY, theta_hw, phi_hw, - result.model_colat); // get fitted residuals!!! + result.model_colatRES); // get fitted residuals!!! if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) { @@ -358,16 +482,16 @@ protected: // inverse (encoder = celestial + pcm) - std::vector colon_res = _colonRES; - std::vector colat_res = _colatRES; + std::vector colon_err = _colonRES; + std::vector colat_err = _colatRES; if constexpr (std::same_as) { theta_tag = theta; phi_tag = phi; } - for (size_t i = 0; i < colat_res.size(); ++i) { - colon_res[i] = -colon_res[i]; - colat_res[i] = -colat_res[i]; + for (size_t i = 0; i < colat_err.size(); ++i) { + colon_err[i] = -colon_err[i]; + colat_err[i] = -colat_err[i]; if constexpr (std::same_as) { theta_tag[i] = _targetCOLAT[i] + MCC_HALF_PI; @@ -375,14 +499,14 @@ protected: } } - result.bspline_fit_err[2] = bsplines::fitpack_sphere_fit(theta_tag, phi_tag, colon_res, 1.0, ty, tx, + result.bspline_fit_err[2] = bsplines::fitpack_sphere_fit(theta_tag, phi_tag, colon_err, 1.0, ty, tx, pcm_data.bspline.inverseCoeffsX, resi2x); if (result.bspline_fit_err[2] > 0) { result.error = MccPCMFitterErrorCode::ERROR_BSPLINE_FIT; return result; } - result.bspline_fit_err[3] = bsplines::fitpack_sphere_fit(theta_tag, phi_tag, colat_res, 1.0, ty, tx, + result.bspline_fit_err[3] = bsplines::fitpack_sphere_fit(theta_tag, phi_tag, colat_err, 1.0, ty, tx, pcm_data.bspline.inverseCoeffsY, resi2y); if (result.bspline_fit_err[3] > 0) { result.error = MccPCMFitterErrorCode::ERROR_BSPLINE_FIT; @@ -390,10 +514,10 @@ protected: } bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.inverseCoeffsX, theta_tag, phi_tag, - result.inv_model_colon); // get fitted residuals!!! + result.inv_model_colonRES); // get fitted residuals!!! bsplines::fitpack_eval_spl2d(ty, tx, pcm_data.bspline.inverseCoeffsY, theta_tag, phi_tag, - result.inv_model_colat); // get fitted residuals!!! + result.inv_model_colatRES); // get fitted residuals!!! } computeModelResi(result); @@ -546,10 +670,10 @@ protected: pcm_data.geomCoefficients.forkFlexure = new_beta(8); } - result.model_colon = {model.begin(), model.begin() + numberOfPoints()}; - result.model_colat = {model.begin() + numberOfPoints(), model.end()}; - result.colon_res.resize(numberOfPoints()); - result.colat_res.resize(numberOfPoints()); + result.model_colonRES = {model.begin(), model.begin() + numberOfPoints()}; + result.model_colatRES = {model.begin() + numberOfPoints(), model.end()}; + result.colon_err.resize(numberOfPoints()); + result.colat_err.resize(numberOfPoints()); result.colon_weight = {weights.begin(), weights.begin() + numberOfPoints()}; result.colat_weight = {weights.begin() + numberOfPoints(), weights.end()}; diff --git a/tests/mcc_pcm_test.cpp b/tests/mcc_pcm_test.cpp index eb21077..2dd07b6 100644 --- a/tests/mcc_pcm_test.cpp +++ b/tests/mcc_pcm_test.cpp @@ -3,6 +3,7 @@ #include +#include struct pcm_res_t { double pcmX, pcmY; @@ -14,9 +15,11 @@ int main() { using pcm_t = mcc::impl::MccDefaultPCM; using pcm_const_t = mcc::impl::MccPCMConstructor; + using pcm_fit_t = mcc::impl::MccPCMFitter; pcm_t pcm; pcm_const_t pcm_const; + pcm_fit_t pcm_fitter; pcm_t::pcm_data_t pcm_data{.type = mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY, .siteLatitude = 43.6466666667_degs, @@ -37,9 +40,12 @@ int main() size_t decM = 10; // number of B-spline inner knots along DEC-axis double ha_step = 360.0_degs / (haM - 1); - pcm_data.bspline.knotsX.resize(haM); // [0, 360] - for (size_t i = 0; i < haM; ++i) { - pcm_data.bspline.knotsX[i] = i * ha_step; + pcm_data.bspline.knotsX.resize(haM); + // for (size_t i = 0; i < haM; ++i) { // [0, 360] + // pcm_data.bspline.knotsX[i] = i * ha_step; + // } + for (size_t i = 0; i < haM; ++i) { // [-180, 180] + pcm_data.bspline.knotsX[i] = i * ha_step - std::numbers::pi; } double dec_start = -35.0_degs; @@ -100,7 +106,8 @@ int main() // add 'noise' - double sigma = 24.3248792_arcsecs / 2.355; + // double sigma = 24.3248792_arcsecs / 2.355; + double sigma = 4.3248792_arcsecs; std::random_device rd{}; std::mt19937 gen{rd()}; @@ -115,14 +122,18 @@ int main() resX[i] = ha[i] - encX[i]; resY[i] = dec[i] - encY[i]; - hadec.from(mcc::impl::MccSkyHADEC_OBS{(double)ha[i], (double)dec[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_const.addPoint(hadec, mcc::impl::MccGenXY{(double)encX[i], (double)encY[i]}); + + pcm_fitter.addPoint((double)ha[i], (double)dec[i], (double)encX[i], (double)encY[i], + mcc::impl::MccCelestialCoordEpoch{}); } pcm_t::pcm_data_t fit_pcm_data{.type = pcm_data.type, .siteLatitude = pcm_data.siteLatitude}; - auto r = pcm_const.computeModel(fit_pcm_data); + // auto r = pcm_const.computeModel(fit_pcm_data); + auto r = pcm_fitter.computeModel(fit_pcm_data); if (r.error) { std::println("error: {}", r.error.message()); return 1; @@ -165,8 +176,10 @@ int main() mcc::impl::MccAngle(pcm_data.geomCoefficients.forkFlexure).arcsecs()); - std::println("\n\n{}", r.colon_res); - std::println("\n\n{}", r.colat_res); + // std::println("\n\n{}", r.colon_res); + // std::println("\n\n{}", r.colat_res); + std::println("\n\n{}", r.model_colonRES); + std::println("\n\n{}", r.model_colatRES); return 0; } \ No newline at end of file diff --git a/tests/mcc_pcm_z1000_test.cpp b/tests/mcc_pcm_z1000_test.cpp index 83cb908..a40343d 100644 --- a/tests/mcc_pcm_z1000_test.cpp +++ b/tests/mcc_pcm_z1000_test.cpp @@ -2,6 +2,7 @@ #include #include +#include static constexpr mcc::MccMountType MOUNT_TYPE{mcc::MccMountType::CROSSAXIS_TYPE}; @@ -10,6 +11,8 @@ using namespace mcc::impl; int main(int narg, char* argv[]) { MccPCMConstructor pcm_cstr; + mcc::impl::MccPCMFitter pcm_fitter; + // MccSkyPoint sp; MccSkyHADEC_OBS hadec; MccGenXY xy; @@ -52,6 +55,8 @@ int main(int narg, char* argv[]) std::println("READ DATA:"); + MccCelestialCoordEpoch ep = MccCelestialCoordEpoch::now(); + size_t i = 0; while (!fst.eof()) { // x - degs, y -degs, diff_x - arcsecs, diff_y - arcsecs @@ -63,29 +68,38 @@ int main(int narg, char* argv[]) tag_ha = x + diff_x / 3600.0; tag_dec = y + diff_y / 3600.0; - hadec = MccSkyHADEC_OBS{MccAngleHA_OBS(tag_ha, mcc_degrees), MccAngleDEC_OBS(tag_dec, mcc_degrees)}; + hadec = MccSkyHADEC_OBS{MccAngleHA_OBS(tag_ha, mcc_degrees), MccAngleDEC_OBS(tag_dec, mcc_degrees), ep}; // sp.from(hadec); - xy = MccGenXY(MccAngleX(x, mcc_degrees), MccAngleY(y, mcc_degrees)); + xy = MccGenXY(MccAngleX(x, mcc_degrees), MccAngleY(y, mcc_degrees), ep); std::println("\t[encX = {:10.6f} encY = {:10.6f} HA = {} DEC = {}]", xy.x().degrees(), xy.y().degrees(), hadec.x().sexagesimal(true), hadec.y().sexagesimal()); // pcm_cstr.addPoint(sp, xy); pcm_cstr.addPoint(hadec, xy); + + pcm_fitter.addPoint(hadec, xy); } fst.close(); auto r = pcm_cstr.computeModel(fit_pcm_data); + auto fr = pcm_fitter.computeModel(fit_pcm_data); - if (r.error) { - std::println("error: {}", r.error.message()); + // if (r.error) { + // std::println("error: {}", r.error.message()); + // return 1; + // } + + if (fr.error) { + std::println("error: {}", fr.error.message()); return 1; } std::println("\n\n{:*^40}\n", " FITTED RESULT (TYPE = GEOMETRY) "); - std::println("\tNUM OF ITERS: {}", r.final_iter); + // std::println("\tNUM OF ITERS: {}", r.final_iter); + std::println("\tNUM OF ITERS: {}", fr.final_iter); std::println("FITTED COEFFS:"); std::println("X0 = {}, Y0 = {}", MccAngleFancyString(fit_pcm_data.geomCoefficients.zeroPointX), @@ -100,39 +114,61 @@ int main(int narg, char* argv[]) std::println("\n\n{:*^40}\n", " FITTED DIFFS "); auto tab = pcm_cstr.getTable(); - for (size_t i = 0; i < pcm_cstr.numberOfPoints(); ++i) { + // for (size_t i = 0; i < pcm_cstr.numberOfPoints(); ++i) { + // std::println("{} {} {} {:6.2f}% ({:5.3f}) {} {} {:6.2f}% ({:5.3f}) (HA = {} DEC = {})", i, + // MccAngleFancyString(tab.colon_res[i]), MccAngleFancyString(r.model_colon[i]), + // std::abs((tab.colon_res[i] - r.model_colon[i]) / tab.colon_res[i]) * 100.0, r.colon_weight[i], + // MccAngleFancyString(tab.colat_res[i]), MccAngleFancyString(r.model_colat[i]), + // std::abs((tab.colat_res[i] - r.model_colat[i]) / tab.colat_res[i]) * 100.0, r.colat_weight[i], + // MccAngle(tab.hw_colon[i]).sexagesimal(true), MccAngle(tab.hw_colat[i]).sexagesimal()); + // } + + auto pcm_tab = pcm_fitter.getPCMTable(); + for (size_t i = 0; i < pcm_fitter.numberOfPoints(); ++i) { std::println("{} {} {} {:6.2f}% ({:5.3f}) {} {} {:6.2f}% ({:5.3f}) (HA = {} DEC = {})", i, - MccAngleFancyString(tab.colon_res[i]), MccAngleFancyString(r.model_colon[i]), - std::abs((tab.colon_res[i] - r.model_colon[i]) / tab.colon_res[i]) * 100.0, r.colon_weight[i], - MccAngleFancyString(tab.colat_res[i]), MccAngleFancyString(r.model_colat[i]), - std::abs((tab.colat_res[i] - r.model_colat[i]) / tab.colat_res[i]) * 100.0, r.colat_weight[i], - MccAngle(tab.hw_colon[i]).sexagesimal(true), MccAngle(tab.hw_colat[i]).sexagesimal()); + MccAngleFancyString(pcm_tab[i].res.x()), MccAngleFancyString(fr.model_colonRES[i]), + std::abs(fr.colon_err[i]) / pcm_tab[i].res.x() * 100.0, fr.colon_weight[i], + MccAngleFancyString(pcm_tab[i].res.y()), MccAngleFancyString(fr.model_colatRES[i]), + std::abs(fr.colat_err[i]) / pcm_tab[i].res.y() * 100.0, fr.colat_weight[i], + pcm_tab[i].hw.x().sexagesimal(true), pcm_tab[i].hw.y().sexagesimal()); } - - std::println("\n\n{:*^40}\n", " FITTED RESULT (TYPE = BSPLINE) "); fit_pcm_data.type = MccDefaultPCMType::PCM_TYPE_BSPLINE; r = pcm_cstr.computeModel(fit_pcm_data); + fr = pcm_fitter.computeModel(fit_pcm_data); - if (r.error) { - std::println("error: {}", r.error.message()); - std::println("b-spline error: {}", r.bspline_fit_err); + // if (r.error) { + // std::println("error: {}", r.error.message()); + // std::println("b-spline error: {}", r.bspline_fit_err); + // return 1; + // } + if (fr.error) { + std::println("error: {}", fr.error.message()); + std::println("b-spline error: {}", fr.bspline_fit_err); return 1; } std::println("\n\n{:*^40}\n", " FITTED DIFFS "); - for (size_t i = 0; i < pcm_cstr.numberOfPoints(); ++i) { - std::println("{} {} {} {:6.2f}% {} {} {:6.2f}% (HA = {} DEC = {})", i, - MccAngleFancyString(tab.colon_res[i]), MccAngleFancyString(r.model_colon[i]), - std::abs((tab.colon_res[i] - r.model_colon[i]) / tab.colon_res[i]) * 100.0, - MccAngleFancyString(tab.colat_res[i]), MccAngleFancyString(r.model_colat[i]), - std::abs((tab.colat_res[i] - r.model_colat[i]) / tab.colat_res[i]) * 100.0, - MccAngle(tab.hw_colon[i]).sexagesimal(true), MccAngle(tab.hw_colat[i]).sexagesimal()); - } + // for (size_t i = 0; i < pcm_cstr.numberOfPoints(); ++i) { + // std::println("{} {} {} {:6.2f}% {} {} {:6.2f}% (HA = {} DEC = {})", i, + // MccAngleFancyString(tab.colon_res[i]), MccAngleFancyString(r.model_colon[i]), + // std::abs((tab.colon_res[i] - r.model_colon[i]) / tab.colon_res[i]) * 100.0, + // MccAngleFancyString(tab.colat_res[i]), MccAngleFancyString(r.model_colat[i]), + // std::abs((tab.colat_res[i] - r.model_colat[i]) / tab.colat_res[i]) * 100.0, + // MccAngle(tab.hw_colon[i]).sexagesimal(true), MccAngle(tab.hw_colat[i]).sexagesimal()); + // } + for (size_t i = 0; i < pcm_fitter.numberOfPoints(); ++i) { + std::println("{} {} {} {:6.2f}% ({:5.3f}) {} {} {:6.2f}% ({:5.3f}) (HA = {} DEC = {})", i, + MccAngleFancyString(pcm_tab[i].res.x()), MccAngleFancyString(fr.model_colonRES[i]), + std::abs(fr.colon_err[i]) / pcm_tab[i].res.x() * 100.0, fr.colon_weight[i], + MccAngleFancyString(pcm_tab[i].res.y()), MccAngleFancyString(fr.model_colatRES[i]), + std::abs(fr.colat_err[i]) / pcm_tab[i].res.y() * 100.0, fr.colat_weight[i], + pcm_tab[i].hw.x().sexagesimal(true), pcm_tab[i].hw.y().sexagesimal()); + } return 0;