416 lines
16 KiB
C++
416 lines
16 KiB
C++
#include <format>
|
|
#include <fstream>
|
|
#include <print>
|
|
#include <ranges>
|
|
|
|
#include <cxxopts.hpp>
|
|
|
|
#include <mcc/mcc_coordinate.h>
|
|
#include <mcc/mcc_pcm_fit.h>
|
|
|
|
#include "asibfm700_configfile.h"
|
|
|
|
// static constexpr mcc::MccMountType MOUNT_TYPE{mcc::MccMountType::CROSSAXIS_TYPE};
|
|
|
|
int main(int argc, char* argv[])
|
|
{
|
|
/* COMMANDLINE OPTS */
|
|
cxxopts::Options options(argv[0], "Astrosib (c) FM700 mount PCM fitter\n");
|
|
|
|
options.allow_unrecognised_options();
|
|
|
|
options.add_options()("h,help", "Print usage");
|
|
|
|
options.add_options()("v,verbose", "Verbose output");
|
|
|
|
options.add_options()("input_file", "Input encoder-celestial coordinate pairs file",
|
|
cxxopts::value<std::vector<std::string>>()->default_value(""));
|
|
|
|
options.add_options()("niter", "Max number of iterations for robust linear regression method",
|
|
cxxopts::value<size_t>()->default_value("100"));
|
|
|
|
// options.positional_help("[input_encoder-celestial_pair_filename]");
|
|
options.positional_help("mount-server-config-filename input-pcm-data-filename results-filename");
|
|
options.parse_positional({"input_file"});
|
|
|
|
mcc::impl::MccPCMFitter<asibfm700::asibfm700MountType> pcm_fitter;
|
|
mcc::impl::MccPCMFitter<asibfm700::asibfm700MountType>::compute_params_t comp_pars;
|
|
|
|
asibfm700::Asibfm700MountConfig mount_cfg;
|
|
|
|
try {
|
|
auto opt_result = options.parse(argc, argv);
|
|
|
|
auto pos_args = opt_result["input_file"].as<std::vector<std::string>>();
|
|
|
|
if (opt_result["help"].count() || argc == 1 || (pos_args.size() < 3)) {
|
|
std::println("{}", options.help());
|
|
|
|
std::println(
|
|
"\tmount-server-config-filename - filename of the ASIB FM-700 mount server configuration\n"
|
|
"\tinput-pcm-data-filename - filename of the input encoder-celestial coordinate pairs data\n"
|
|
"\tresults-filename - filename of the fitting results\n");
|
|
|
|
// std::println(
|
|
// "[input_encoder-celestial_pair_filename] - Input encoder-celestial coordinate pairs filename. "
|
|
// "It must be in the format: ENCODER_HA ENCODER_DEC ENCODER_EPOCH RA_ICRS DEC_ICRS TEMP(C)
|
|
// PRESSURE(hPa) " "HUMIDITY([0-1])", options.help());
|
|
|
|
return 0;
|
|
}
|
|
|
|
std::string cfg_fname = pos_args[0];
|
|
std::string pcm_data_fname = pos_args[1];
|
|
std::string result_fname = pos_args[2];
|
|
|
|
bool verbose = false;
|
|
if (opt_result["verbose"].count()) {
|
|
verbose = true;
|
|
}
|
|
|
|
auto l_err = mount_cfg.load(cfg_fname);
|
|
if (l_err) {
|
|
std::println("Cannot load mount config: {}", l_err.message());
|
|
return 2;
|
|
}
|
|
|
|
asibfm700::Asibfm700PCM::pcm_data_t pcm_data = mount_cfg.pcmData();
|
|
|
|
std::ifstream fst;
|
|
|
|
fst.open(pcm_data_fname);
|
|
if (!fst.is_open()) {
|
|
std::println("Cannot open input file {}", pcm_data_fname);
|
|
return 2;
|
|
}
|
|
|
|
size_t sz;
|
|
double temp, press, humi;
|
|
std::optional<double> num;
|
|
std::string str;
|
|
std::string fmt_head;
|
|
|
|
std::string_view delim{" "};
|
|
|
|
std::string_view sv;
|
|
std::array<std::string_view, 8> tokens;
|
|
|
|
mcc::impl::MccCelestialCoordEpoch ep;
|
|
mcc::impl::MccAngleX enc_ha;
|
|
mcc::impl::MccAngleY enc_dec;
|
|
mcc::impl::MccAngleRA_ICRS ra;
|
|
mcc::impl::MccAngleDEC_ICRS dec;
|
|
|
|
mcc::impl::MccSkyPoint sp;
|
|
mcc::impl::MccError err;
|
|
|
|
if (verbose) {
|
|
std::format_to(std::back_inserter(fmt_head), "{:^12}{}{:^12}{}{:^11}{}{:^11}{}{:^12}{}{:^6}{}{:^7}{}{:^4}",
|
|
"ENC_HA", delim, "ENC_DEC", delim, "ENC_MJD", delim, "RA_ICRS", delim, "DEC_ICRS", delim,
|
|
"TEMP", delim, "PRESS", delim, "HUMI");
|
|
|
|
str = std::string(fmt_head.size(), '*');
|
|
std::println("{}", str);
|
|
|
|
auto fmt = std::format("*{{:^{}}}*", fmt_head.size() - 2);
|
|
std::println("{}", std::vformat(std::string_view(fmt.begin(), fmt.end()),
|
|
std::make_format_args(" INPUT HARDWARE-CELESTIAL COORDINATE PAIRS ")));
|
|
|
|
fmt = std::format("*{{:<{}}}*", fmt_head.size() - 2);
|
|
auto fmt_sv = std::string_view(fmt.begin(), fmt.end());
|
|
std::println("{}", std::vformat(fmt_sv, std::make_format_args("")));
|
|
auto s = std::format("{} {}", " SITE LAT:", mount_cfg.siteLatitude().sexagesimal());
|
|
std::println("{}", std::vformat(fmt_sv, std::make_format_args(s)));
|
|
s = std::format("{} {}", " SITE LON:", mount_cfg.siteLongitude().sexagesimal());
|
|
std::println("{}", std::vformat(fmt_sv, std::make_format_args(s)));
|
|
s = std::format("{} {} meters", " SITE ELEV:", mount_cfg.siteElevation());
|
|
std::println("{}", std::vformat(fmt_sv, std::make_format_args(s)));
|
|
|
|
s = std::format("{} {} ", " PCM TYPE:", mcc::impl::mccDefaultPCMTypeString(pcm_data.type));
|
|
std::println("{}", std::vformat(fmt_sv, std::make_format_args(s)));
|
|
|
|
std::println("{}", str);
|
|
std::println("{}", fmt_head);
|
|
std::println("{}", std::string(fmt_head.size(), '-'));
|
|
}
|
|
|
|
while (std::getline(fst, str)) {
|
|
// while (!fst.eof()) {
|
|
// std::getline(fst, str);
|
|
|
|
sv = mcc::utils::trimSpaces(str);
|
|
if (!sv.size()) { // an empty string
|
|
continue;
|
|
}
|
|
if (sv[0] == '#') { // comment
|
|
continue;
|
|
}
|
|
|
|
auto toks = std::views::split(sv, std::string_view(" ")) |
|
|
std::views::filter([](auto const& r) { return std::ranges::size(r); });
|
|
|
|
sz = std::ranges::distance(toks.begin(), toks.end());
|
|
if (sz < 8) {
|
|
std::println("Invalid input file format! Number of tokens must be at least 8: {}", str);
|
|
return 3;
|
|
}
|
|
|
|
size_t i = 0;
|
|
for (auto const& t : toks) {
|
|
tokens[i++] = {t.begin(), t.end()};
|
|
}
|
|
|
|
// degrees or sexagesimal hours representations
|
|
enc_ha = {tokens[0], mcc::impl::mcc_hms};
|
|
|
|
// degrees or sexagesimal degrees representations
|
|
enc_dec = {tokens[1]};
|
|
|
|
if (!ep.fromCharRange(tokens[2])) {
|
|
std::println("Invalid input file format! Invalid encoder coordinates epoch representation: {}",
|
|
tokens[2]);
|
|
return 4;
|
|
}
|
|
|
|
// degrees or sexagesimal hours representations
|
|
ra = {tokens[3], mcc::impl::mcc_hms};
|
|
|
|
// degrees or sexagesimal degrees representations
|
|
dec = {tokens[4]};
|
|
|
|
num = mcc::utils::numFromStr<double>(tokens[5]);
|
|
if (!num) {
|
|
std::println("Invalid input file format! Non-numeric representation of temperature: {}", tokens[5]);
|
|
return 5;
|
|
}
|
|
temp = num.value();
|
|
|
|
num = mcc::utils::numFromStr<double>(tokens[6]);
|
|
if (!num) {
|
|
std::println("Invalid input file format! Non-numeric representation of pressure: {}", tokens[6]);
|
|
return 6;
|
|
}
|
|
press = num.value();
|
|
|
|
num = mcc::utils::numFromStr<double>(tokens[7]);
|
|
if (!num) {
|
|
std::println("Invalid input file format! Non-numeric representation of humidity: {}", tokens[7]);
|
|
return 7;
|
|
}
|
|
humi = num.value();
|
|
|
|
// update meteo parameters here
|
|
mcc::impl::MccSkyPoint::cctEngine.updateMeteoERFA(
|
|
{.temperature = temp, .humidity = humi, .pressure = press});
|
|
|
|
sp.from(mcc::impl::MccSkyRADEC_ICRS{ra, dec});
|
|
|
|
// astrometric transformations are here (and it needs meteo updated above!!!)
|
|
err = pcm_fitter.addPoint(sp, mcc::impl::MccGenXY{enc_ha, enc_dec, ep});
|
|
if (err) {
|
|
std::println("An error occured: {}", err.message());
|
|
return 8;
|
|
}
|
|
|
|
if (verbose) {
|
|
std::println("{:>12}{}{:>12}{}{:>11.5f}{}{:>11}{}{:>12}{}{:>5.1f}{}{:>7.2f}{}{:>4.2f}",
|
|
enc_ha.sexagesimal(true), delim, enc_dec.sexagesimal(), delim, ep.MJD(), delim,
|
|
ra.sexagesimal(true), delim, dec.sexagesimal(), delim, temp, delim, press, delim, humi);
|
|
}
|
|
}
|
|
|
|
fst.close();
|
|
|
|
// fitting
|
|
|
|
if (verbose) {
|
|
std::println("\n");
|
|
|
|
fmt_head.clear();
|
|
std::format_to(std::back_inserter(fmt_head), "{:^12}{}{:^12}{}{:^12}{}{:^12}{}{:^14}{}{:^14}", "ENC_HA",
|
|
delim, "ENC_DEC", delim, "OBS_HA", delim, "OBS_DEC", delim, "dHA(OBS-ENC)", delim,
|
|
"dDEC(OBS-ENC)");
|
|
|
|
str = std::string(fmt_head.size(), '*');
|
|
std::println("{}", str);
|
|
|
|
auto fmt = std::format("*{{:^{}}}*", fmt_head.size() - 2);
|
|
std::println("{}",
|
|
std::vformat(std::string_view(fmt.begin(), fmt.end()),
|
|
std::make_format_args(" HARDWARE-CELESTIAL COORDINATE PAIRS TO BE FITTED ")));
|
|
|
|
std::println("{}", str);
|
|
|
|
std::println("{}", fmt_head);
|
|
std::println("{}", std::string(fmt_head.size(), '-'));
|
|
|
|
|
|
fmt = std::format("*{{:<{}}}*", fmt_head.size() - 2);
|
|
|
|
auto tab = pcm_fitter.getPCMTable();
|
|
|
|
for (auto const& el : tab) {
|
|
std::println("{:>12}{}{:>12}{}{:>12}{}{:>12}{}{:>14}{}{:>14}", el.hw.x().sexagesimal(true), delim,
|
|
el.hw.y().sexagesimal(), delim, el.target.x().sexagesimal(true), delim,
|
|
el.target.y().sexagesimal(), delim, mcc::impl::MccAngleFancyString(el.res.x()), delim,
|
|
mcc::impl::MccAngleFancyString(el.res.y()));
|
|
}
|
|
}
|
|
|
|
// return 0;
|
|
|
|
if (pcm_data.type == mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY
|
|
#ifdef USE_BSPLINE_PCM
|
|
|| pcm_data.type == mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE
|
|
#endif
|
|
) {
|
|
if (opt_result["niter"].count()) {
|
|
comp_pars.max_iter = opt_result["niter"].as<size_t>() ? opt_result["niter"].as<size_t>() : 100;
|
|
}
|
|
}
|
|
|
|
auto comp_result = pcm_fitter.computeModel(pcm_data, comp_pars);
|
|
|
|
if (comp_result.error) {
|
|
std::println("An error occured while fit PCM data: {}", comp_result.error.message());
|
|
#ifdef USE_BSPLINE_PCM
|
|
if (comp_result.pcm_type == mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE ||
|
|
comp_result.pcm_type == mcc::impl::MccDefaultPCMType::PCM_TYPE_BSPLINE) {
|
|
std::println("\tB-spline fitting error codes: {}", comp_result.bspline_fit_err);
|
|
}
|
|
#endif
|
|
return 200;
|
|
}
|
|
|
|
auto tab = pcm_fitter.getPCMTable();
|
|
|
|
std::ofstream ofst;
|
|
|
|
ofst.open(result_fname, std::ios_base::trunc);
|
|
if (!ofst.is_open()) {
|
|
std::println("<<< CANNOT OPEN RESULT FILE: {}! >>>", result_fname);
|
|
verbose = true; // print result to the console
|
|
}
|
|
|
|
|
|
fmt_head = std::format("{:^12}{}{:^11}{}{:^12}{}{:^11}{}{:^12}{}{:^12}{}{:^12}{}{:^12}{}{:^7}{}{:^8}",
|
|
"ENC_HA(degs)", delim, "ENC_DEC(degs)", delim, "OBS_HA(degs)", delim, "OBS_DEC(degs)",
|
|
delim, "dHA(OBS-ENC, degs)", delim, "dDEC(OBS-ENC, degs)", delim, "dHA_FIT(degs)", delim,
|
|
"dDEC_FIT(degs)", delim, "EPS_dHA", delim, "EPS_dDEC");
|
|
|
|
if (comp_result.pcm_type == mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY
|
|
#ifdef USE_BSPLINE_PCM
|
|
|| pcm_data.type == mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE
|
|
#endif
|
|
) { // Tukey's weights of the points computed during fitting
|
|
std::format_to(std::back_inserter(fmt_head), "{}{:^7}{}{:^8}", delim, "dHA_WEI", delim, "dDEC_WEI");
|
|
}
|
|
#ifdef USE_BSPLINE_PCM
|
|
else if (comp_result.pcm_type == mcc::impl::MccDefaultPCMType::PCM_TYPE_BSPLINE) { // inverse PCM fitting
|
|
std::format_to(std::back_inserter(fmt_head), "{}{:^17}{}{:^18}{}{:^12}{}{:^12}", delim, "dHA_FIT_INV(degs)",
|
|
delim, "dDEC_FIT_INV(degs)", delim, "EPS_dHA_INV", delim, "EPS_dDEC_INV");
|
|
}
|
|
#endif
|
|
|
|
sv = "POINTING CORRECTION MODEL FITTING RESULTS";
|
|
std::vector<std::string> info{
|
|
std::format("{} {}", " SITE LAT:", mount_cfg.siteLatitude().sexagesimal()),
|
|
std::format("{} {}", " SITE LON:", mount_cfg.siteLongitude().sexagesimal()),
|
|
std::format("{} {} meters", " SITE ELEV:", mount_cfg.siteElevation()),
|
|
std::format("{} {} ", " PCM TYPE:", mcc::impl::mccDefaultPCMTypeString(pcm_data.type))};
|
|
|
|
|
|
if (ofst.is_open()) {
|
|
ofst << "#\n";
|
|
ofst << std::format("# {} ({} UTC)\n", sv, std::chrono::system_clock::now());
|
|
ofst << "#\n";
|
|
for (auto const& s : info) {
|
|
ofst << std::format("# {}\n", s);
|
|
}
|
|
|
|
ofst << "#\n# Format:\n";
|
|
ofst << std::format("# {}\n", fmt_head);
|
|
ofst << "#\n";
|
|
}
|
|
|
|
if (verbose) {
|
|
std::println("\n");
|
|
|
|
str = std::string(fmt_head.size(), '*');
|
|
std::println("{}", str);
|
|
|
|
auto fmt = std::format("*{{:^{}}}*", fmt_head.size() - 2);
|
|
auto fmt_sv = std::string_view(fmt.begin(), fmt.end());
|
|
|
|
std::println("{}", std::vformat(fmt_sv, std::make_format_args(sv)));
|
|
|
|
std::println("{}", std::vformat(fmt_sv, std::make_format_args("")));
|
|
|
|
for (auto const& s : info) {
|
|
std::println("{}", std::vformat(fmt_sv, std::make_format_args(s)));
|
|
}
|
|
|
|
std::println("{}", str);
|
|
|
|
std::println("{}", fmt_head);
|
|
std::println("{}", std::string(fmt_head.size(), '-'));
|
|
}
|
|
|
|
size_t i = 0;
|
|
double dha_eps, ddec_eps;
|
|
for (auto const& el : tab) {
|
|
dha_eps = std::abs(comp_result.colon_err[i] / comp_result.model_colonRES[i]);
|
|
ddec_eps = std::abs(comp_result.colat_err[i] / comp_result.model_colatRES[i]);
|
|
|
|
str = std::format(
|
|
"{:>12.7f}{}{:>11.7f}{}{:>12.7f}{}{:>11.7f}{}{:>12.7f}{}{:>12.7f}{}{:>12.7f}{}{:>12.7f}{}{:>7.3f}",
|
|
el.hw.x().degrees(), delim, el.hw.y().degrees(), delim, el.target.x().degrees(), delim,
|
|
el.target.y().degrees(), delim, el.res.x().degrees(), delim, el.res.y().degrees(), delim,
|
|
comp_result.model_colonRES[i] * mcc::MCC_RADS_TO_DEGRESS, delim,
|
|
comp_result.model_colatRES[i] * mcc::MCC_RADS_TO_DEGRESS, delim, dha_eps, delim, ddec_eps);
|
|
|
|
if (comp_result.pcm_type == mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY
|
|
#ifdef USE_BSPLINE_PCM
|
|
|| pcm_data.type == mcc::impl::MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE
|
|
#endif
|
|
) {
|
|
std::format_to(std::back_inserter(str), "{}{:>7.3f}{}{:>8.3f}", delim, comp_result.colon_weight[i],
|
|
delim, comp_result.colat_weight[i]);
|
|
}
|
|
#ifdef USE_BSPLINE_PCM
|
|
else if (comp_result.pcm_type == mcc::impl::MccDefaultPCMType::PCM_TYPE_BSPLINE) { // inverse PCM fitting
|
|
std::format_to(std::back_inserter(str), "{}{:^17.7}{}{:^18.7}{}{:^12.7}{}{:^12.7}", delim,
|
|
comp_result.inv_model_colonRES[i], delim, comp_result.inv_model_colatRES[i], delim,
|
|
comp_result.inv_colon_err[i], delim, comp_result.inv_colat_err[i]);
|
|
}
|
|
#endif
|
|
|
|
++i;
|
|
|
|
if (ofst.is_open()) {
|
|
ofst << std::format(" {}\n", str);
|
|
}
|
|
|
|
if (verbose) {
|
|
std::println("{}", str);
|
|
}
|
|
}
|
|
|
|
ofst.close();
|
|
} catch (cxxopts::exceptions::parsing& ex) {
|
|
std::println("An error occured while parsing input options: {}", ex.what());
|
|
|
|
return 1;
|
|
} catch (std::exception& ex) {
|
|
std::println("An exception occured: {}", ex.what());
|
|
|
|
return 10;
|
|
|
|
} catch (...) {
|
|
std::println("An unhandled exception occured!");
|
|
|
|
return 100;
|
|
}
|
|
|
|
return 0;
|
|
} |