diff --git a/include/mcc/mcc_pcm_construct.h b/include/mcc/mcc_pcm_construct.h index 7eaef56..41698f4 100644 --- a/include/mcc/mcc_pcm_construct.h +++ b/include/mcc/mcc_pcm_construct.h @@ -409,7 +409,7 @@ protected: Eigen::VectorXd r = mask.select(1.0 - u * u, 0.0); - tukey_weight = r * r; + tukey_weight = r.array() * r.array(); }; auto gastwirth_est = [](Eigen::VectorXd const& resi) { @@ -422,6 +422,10 @@ protected: // add 1.0E-8 to avoid 0.0 return (0.3 * arr(arr.size() / 3) + 0.4 * arr(arr.size() / 2) + 0.3 * arr(2 * arr.size() / 3)) / 0.6745 + 1.0E-8; + + // std::nth_element(arr.begin(), arr.begin() + arr.size() / 2, arr.end()); + + // return arr(arr.size() / 2) / 0.6745 + 1.0E-8; }; compute_result_t result{.pcm_type = pcm_data.type, .error = MccDefaultPCMConstructorErrorCode::ERROR_OK}; @@ -429,15 +433,16 @@ protected: 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); + Eigen::VectorXd y(2 * numberOfPoints()), beta = Eigen::VectorXd::Zero(n_coeffs), new_beta(n_coeffs); Eigen::RowVectorXd row(n_coeffs); - double tgY, cosY, rcosY, cosX, sinX, cosPhi = std::cos(pcm_data.siteLatitude), - sinPhi = std::sin(pcm_data.siteLatitude); + double tgY, sinY, cosY, rcosY, cosX, sinX, cosPhi = std::cos(pcm_data.siteLatitude), + sinPhi = std::sin(pcm_data.siteLatitude); // system matrix and right-hand vector for (size_t i = 0; i < numberOfPoints(); ++i) { cosY = std::cos(_table.target_colat[i]); + sinY = std::sin(_table.target_colat[i]); rcosY = 1.0 / cosY; tgY = std::tan(_table.target_colat[i]); cosX = std::cos(_table.target_colon[i]); @@ -466,15 +471,12 @@ protected: if constexpr (MOUNT_TYPE == MccMountType::FORK_TYPE) { if (utils::isEqual(cosX, 0.0)) { - row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, - cosPhi * cosX * std::sin(_table.target_colat[i]) - sinPhi * cosY, 0.0, 0.0; + row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * sinY - sinPhi * cosY, 0.0, 0.0; } else { - row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, - cosPhi * cosX * std::sin(_table.target_colat[i]) - sinPhi * cosY, 0.0, 1.0 / cosX; + row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * sinY - sinPhi * cosY, 0.0, 1.0 / cosX; } } else { - row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * std::sin(_table.target_colat[i]) - sinPhi * cosY, - 0.0; + row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * sinY - sinPhi * cosY, 0.0; } X.row(i + numberOfPoints()) = row; @@ -504,6 +506,8 @@ protected: } beta = new_beta; + + resi = y - X * beta; } pcm_data.geomCoefficients.zeroPointX = new_beta(0); @@ -521,6 +525,8 @@ protected: Eigen::VectorXd model = X * new_beta; 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()); for (size_t i = 0; i < numberOfPoints(); ++i) { result.colon_res[i] = _table.colon_res[i] - result.model_colon[i]; // = target - model diff --git a/tests/mcc_pcm_test.cpp b/tests/mcc_pcm_test.cpp index 1bb2f15..3a852dd 100644 --- a/tests/mcc_pcm_test.cpp +++ b/tests/mcc_pcm_test.cpp @@ -20,8 +20,10 @@ int main() 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, + // .geomCoefficients = {.zeroPointX = 2.434534_degs, + // .zeroPointY = -3.382549_degs, + .geomCoefficients = {.zeroPointX = 2.434534_arcsecs, + .zeroPointY = -3.382549_arcsecs, .collimationErr = 10.534_arcsecs, .nonperpendErr = 21.86834_arcsecs, .misalignErr1 = 7.345768_arcsecs, @@ -31,7 +33,9 @@ int main() .forkFlexure = 9.342354_arcsecs}}; - size_t N = 50; + pcm.setPCMData(pcm_data); + + size_t N = 100; std::vector encX(N), encY(N), ha(N), dec(N), pcmX(N), pcmY(N); std::vector resX(N), resY(N); @@ -57,11 +61,8 @@ int main() 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 start_encY = -32.2423_degs, stop_encY = 84.896124_degs; double x_step = (stop_encX - start_encX) / (N - 1), y_step = (stop_encY - start_encY) / (N - 1); @@ -80,6 +81,8 @@ int main() // add 'noise' + double sigma = 24.3248792_arcsecs / 2.355; + std::random_device rd{}; std::mt19937 gen{rd()}; std::normal_distribution dist{0.0, sigma}; @@ -106,8 +109,44 @@ int main() return 1; } + std::println("FITTED STATUS:"); + std::println("\tNUM OF ITERS: {}", r.final_iter); + std::println("FITTED COEFFS:"); - std::println("X0 = {}, Y0 = ", fit_pcm_data.geomCoefficients.zeroPointX, fit_pcm_data.geomCoefficients.zeroPointY); + std::println("X0 = {} ({}), Y0 = {} ({})", mcc::impl::MccAngleFancyString(fit_pcm_data.geomCoefficients.zeroPointX), + mcc::impl::MccAngleFancyString(pcm_data.geomCoefficients.zeroPointX), + mcc::impl::MccAngleFancyString(fit_pcm_data.geomCoefficients.zeroPointY), + mcc::impl::MccAngleFancyString(pcm_data.geomCoefficients.zeroPointY)); + + // std::println("X0 = {} ({}), Y0 = {} ({})", + // mcc::impl::MccAngle(fit_pcm_data.geomCoefficients.zeroPointX).degrees(), + // mcc::impl::MccAngle(pcm_data.geomCoefficients.zeroPointX).degrees(), + // mcc::impl::MccAngle(fit_pcm_data.geomCoefficients.zeroPointY).degrees(), + // mcc::impl::MccAngle(pcm_data.geomCoefficients.zeroPointY).degrees()); + + std::println("collimErr = {} ({})", mcc::impl::MccAngle(fit_pcm_data.geomCoefficients.collimationErr).arcsecs(), + mcc::impl::MccAngle(pcm_data.geomCoefficients.collimationErr).arcsecs()); + + std::println("nonperpErr = {} ({})", mcc::impl::MccAngle(fit_pcm_data.geomCoefficients.nonperpendErr).arcsecs(), + mcc::impl::MccAngle(pcm_data.geomCoefficients.nonperpendErr).arcsecs()); + + std::println("misalignErr1 = {} ({}), misalignErr2 = {} ({})", + mcc::impl::MccAngle(fit_pcm_data.geomCoefficients.misalignErr1).arcsecs(), + mcc::impl::MccAngle(pcm_data.geomCoefficients.misalignErr1).arcsecs(), + mcc::impl::MccAngle(fit_pcm_data.geomCoefficients.misalignErr2).arcsecs(), + mcc::impl::MccAngle(pcm_data.geomCoefficients.misalignErr2).arcsecs()); + + std::println("tubeflexErr = {} ({})", mcc::impl::MccAngle(fit_pcm_data.geomCoefficients.tubeFlexure).arcsecs(), + mcc::impl::MccAngle(pcm_data.geomCoefficients.tubeFlexure).arcsecs()); + + std::println("DECaxflexErr = {} ({})", mcc::impl::MccAngle(fit_pcm_data.geomCoefficients.DECaxisFlexure).arcsecs(), + mcc::impl::MccAngle(pcm_data.geomCoefficients.DECaxisFlexure).arcsecs()); + + std::println("forkflexErr = {} ({})", mcc::impl::MccAngle(fit_pcm_data.geomCoefficients.forkFlexure).arcsecs(), + mcc::impl::MccAngle(pcm_data.geomCoefficients.forkFlexure).arcsecs()); + + + std::println("\n\n{}", r.colon_res); return 0; } \ No newline at end of file