timur/mcc
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

...

Browse Source
This commit is contained in:
Timur A. Fatkhullin
2026-03-24 16:45:37 +03:00
parent 0182844af6
commit 2c561f7259
2 changed files with 63 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
include/mcc/mcc_pcm_construct.h
View File

@@ -409,7 +409,7 @@ protected:
Eigen::VectorXd r = mask.select(1.0 - u * u, 0.0); 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) { auto gastwirth_est = [](Eigen::VectorXd const& resi) {
@@ -422,6 +422,10 @@ protected:
// add 1.0E-8 to avoid 0.0 // 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 + 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; 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}; 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; const size_t n_coeffs = MOUNT_TYPE == MccMountType::FORK_TYPE ? 9 : 8;
Eigen::MatrixXd X(2 * numberOfPoints(), n_coeffs); 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); Eigen::RowVectorXd row(n_coeffs);
double tgY, cosY, rcosY, cosX, sinX, cosPhi = std::cos(pcm_data.siteLatitude), double tgY, sinY, cosY, rcosY, cosX, sinX, cosPhi = std::cos(pcm_data.siteLatitude),
sinPhi = std::sin(pcm_data.siteLatitude); sinPhi = std::sin(pcm_data.siteLatitude);
// system matrix and right-hand vector // system matrix and right-hand vector
for (size_t i = 0; i < numberOfPoints(); ++i) { for (size_t i = 0; i < numberOfPoints(); ++i) {
cosY = std::cos(_table.target_colat[i]); cosY = std::cos(_table.target_colat[i]);
sinY = std::sin(_table.target_colat[i]);
rcosY = 1.0 / cosY; rcosY = 1.0 / cosY;
tgY = std::tan(_table.target_colat[i]); tgY = std::tan(_table.target_colat[i]);
cosX = std::cos(_table.target_colon[i]); cosX = std::cos(_table.target_colon[i]);
@@ -466,15 +471,12 @@ protected:
if constexpr (MOUNT_TYPE == MccMountType::FORK_TYPE) { if constexpr (MOUNT_TYPE == MccMountType::FORK_TYPE) {
if (utils::isEqual(cosX, 0.0)) { if (utils::isEqual(cosX, 0.0)) {
row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * sinY - sinPhi * cosY, 0.0, 0.0;
cosPhi * cosX * std::sin(_table.target_colat[i]) - sinPhi * cosY, 0.0, 0.0;
} else { } else {
row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * sinY - sinPhi * cosY, 0.0, 1.0 / cosX;
cosPhi * cosX * std::sin(_table.target_colat[i]) - sinPhi * cosY, 0.0, 1.0 / cosX;
} }
} else { } else {
row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * std::sin(_table.target_colat[i]) - sinPhi * cosY, row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * sinY - sinPhi * cosY, 0.0;
0.0;
} }
X.row(i + numberOfPoints()) = row; X.row(i + numberOfPoints()) = row;
@@ -504,6 +506,8 @@ protected:
} }
beta = new_beta; beta = new_beta;
resi = y - X * beta;
} }
pcm_data.geomCoefficients.zeroPointX = new_beta(0); pcm_data.geomCoefficients.zeroPointX = new_beta(0);
@@ -521,6 +525,8 @@ protected:
Eigen::VectorXd model = X * new_beta; Eigen::VectorXd model = X * new_beta;
result.model_colon = {model.begin(), model.begin() + numberOfPoints()}; result.model_colon = {model.begin(), model.begin() + numberOfPoints()};
result.model_colat = {model.begin() + numberOfPoints(), model.end()}; 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) { for (size_t i = 0; i < numberOfPoints(); ++i) {
result.colon_res[i] = _table.colon_res[i] - result.model_colon[i]; // = target - model result.colon_res[i] = _table.colon_res[i] - result.model_colon[i]; // = target - model

55
tests/mcc_pcm_test.cpp
View File

@@ -20,8 +20,10 @@ int main()
pcm_t::pcm_data_t pcm_data{.type = mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY, pcm_t::pcm_data_t pcm_data{.type = mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY,
.siteLatitude = 43.6466666667_degs, .siteLatitude = 43.6466666667_degs,
.geomCoefficients = {.zeroPointX = 2.434534_degs, // .geomCoefficients = {.zeroPointX = 2.434534_degs,
.zeroPointY = -3.382549_degs, // .zeroPointY = -3.382549_degs,
.geomCoefficients = {.zeroPointX = 2.434534_arcsecs,
.zeroPointY = -3.382549_arcsecs,
.collimationErr = 10.534_arcsecs, .collimationErr = 10.534_arcsecs,
.nonperpendErr = 21.86834_arcsecs, .nonperpendErr = 21.86834_arcsecs,
.misalignErr1 = 7.345768_arcsecs, .misalignErr1 = 7.345768_arcsecs,
@@ -31,7 +33,9 @@ int main()
.forkFlexure = 9.342354_arcsecs}}; .forkFlexure = 9.342354_arcsecs}};
size_t N = 50; pcm.setPCMData(pcm_data);
size_t N = 100;
std::vector<mcc::impl::MccAngle> encX(N), encY(N), ha(N), dec(N), pcmX(N), pcmY(N); std::vector<mcc::impl::MccAngle> encX(N), encY(N), ha(N), dec(N), pcmX(N), pcmY(N);
std::vector<mcc::impl::MccAngle> resX(N), resY(N); std::vector<mcc::impl::MccAngle> resX(N), resY(N);
@@ -57,11 +61,8 @@ int main()
return 0; return 0;
}; };
pcm.setPCMData(pcm_data);
double start_encX = 1.2423_degs, stop_encX = 354.896124_degs; double start_encX = 1.2423_degs, stop_encX = 354.896124_degs;
double start_encY = -32.2423_degs, stop_encY = 89.896124_degs; double start_encY = -32.2423_degs, stop_encY = 84.896124_degs;
double sigma = 24.3248792_arcsecs;
double x_step = (stop_encX - start_encX) / (N - 1), y_step = (stop_encY - start_encY) / (N - 1); 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' // add 'noise'
double sigma = 24.3248792_arcsecs / 2.355;
std::random_device rd{}; std::random_device rd{};
std::mt19937 gen{rd()}; std::mt19937 gen{rd()};
std::normal_distribution dist{0.0, sigma}; std::normal_distribution dist{0.0, sigma};
@@ -106,8 +109,44 @@ int main()
return 1; return 1;
} }
std::println("FITTED STATUS:");
std::println("\tNUM OF ITERS: {}", r.final_iter);
std::println("FITTED COEFFS:"); 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; return 0;
} }