#include <fstream>
#include <print>

#include <mcc/mcc_pcm_construct.h>
#include <mcc/mcc_pcm_fit.h>

static constexpr mcc::MccMountType MOUNT_TYPE{mcc::MccMountType::CROSSAXIS_TYPE};

using namespace mcc::impl;

int main(int narg, char* argv[])
{
    MccPCMConstructor<MOUNT_TYPE> pcm_cstr;
    mcc::impl::MccPCMFitter<MOUNT_TYPE> pcm_fitter;

    // MccSkyPoint sp;
    MccSkyHADEC_OBS hadec;
    MccGenXY xy;

    MccDefaultPCM<MOUNT_TYPE>::pcm_data_t fit_pcm_data;

    size_t haM = 10;   // number of B-spline inner knots along HA-axis
    size_t decM = 10;  // number of B-spline inner knots along DEC-axis

    double ha_step = 360.0_degs / (haM - 1);
    fit_pcm_data.bspline.knotsX.resize(haM);
    for (size_t i = 0; i < haM; ++i) {
        // fit_pcm_data.bspline.knotsX[i] = i * ha_step; // [0, 360]
        fit_pcm_data.bspline.knotsX[i] = i * ha_step - std::numbers::pi;  // [-180, 180]
    }

    double dec_start = -25.0_degs;
    double dec_step = (90.0_degs - dec_start) / (decM - 1);
    fit_pcm_data.bspline.knotsY.resize(decM);
    for (size_t i = 0; i < decM; ++i) {
        fit_pcm_data.bspline.knotsY[i] = dec_start + i * dec_step;
    }


    std::ifstream fst;

    std::string fname = "z1000_pcm_measu.data";

    if (narg > 1) {
        fname = argv[1];
    }


    fst.open(fname);
    if (!fst.is_open()) {
        std::println("Cannot open file {}", fname);
        return 1;
    }

    double x, y, diff_x, diff_y, tag_ha, tag_dec;

    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
        fst >> x >> y >> diff_x >> diff_y;

        // std::println("\t{}\t{:10.6f} {:10.6f}   {:8.6f} {:8.6f}", i++, x, y, diff_x, diff_y);
        std::print("\t{}\t{:10.6f} {:10.6f}   {:8.6f} {:8.6f}", i++, x, y, diff_x, diff_y);

        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), ep};
        // sp.from(hadec);

        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());
    //     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: {}", fr.final_iter);

    std::println("FITTED COEFFS:");
    std::println("X0 = {}, Y0 = {}", MccAngleFancyString(fit_pcm_data.geomCoefficients.zeroPointX),
                 MccAngleFancyString(fit_pcm_data.geomCoefficients.zeroPointY));

    std::println("collimErr = {}", MccAngleFancyString(fit_pcm_data.geomCoefficients.collimationErr));
    std::println("nonperpErr = {}", MccAngleFancyString(fit_pcm_data.geomCoefficients.nonperpendErr));
    std::println("misalignErr1 = {}", MccAngleFancyString(fit_pcm_data.geomCoefficients.misalignErr1));
    std::println("misalignErr2 = {}", MccAngleFancyString(fit_pcm_data.geomCoefficients.misalignErr2));
    std::println("tubeflexErr = {}", MccAngleFancyString(fit_pcm_data.geomCoefficients.tubeFlexure));
    std::println("DECaxflexErr = {}", MccAngleFancyString(fit_pcm_data.geomCoefficients.DECaxisFlexure));

    std::println("\n\n{:*^40}\n", " FITTED DIFFS ");
    auto tab = pcm_cstr.getTable();
    // 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(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);
    //     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_fitter.numberOfPoints(); ++i) {
        std::println("{}     {} {} {:6.2f}%      {} {} {:6.2f}%    (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, MccAngleFancyString(pcm_tab[i].res.y()),
                     MccAngleFancyString(fr.model_colatRES[i]), std::abs(fr.colat_err[i]) / pcm_tab[i].res.y() * 100.0,
                     pcm_tab[i].hw.x().sexagesimal(true), pcm_tab[i].hw.y().sexagesimal());
    }


    return 0;
}