#include <print>
#include <random>


#include <mcc/mcc_pcm_construct.h>

struct pcm_res_t {
    double pcmX, pcmY;
};

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

int main()
{
    using pcm_t = mcc::impl::MccDefaultPCM<MOUNT_TYPE>;
    using pcm_const_t = mcc::impl::MccPCMConstructor<MOUNT_TYPE>;

    pcm_t pcm;
    pcm_const_t pcm_const;

    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,
                                                    .collimationErr = 10.534_arcsecs,
                                                    .nonperpendErr = 21.86834_arcsecs,
                                                    .misalignErr1 = 7.345768_arcsecs,
                                                    .misalignErr2 = -12.87934_arcsecs,
                                                    .tubeFlexure = 5.9083_arcsecs,
                                                    .DECaxisFlexure = 3.48379812_arcsecs,
                                                    .forkFlexure = 9.342354_arcsecs}};


    size_t N = 50;

    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);
    pcm_res_t pcm_res;



    auto comp_hadec = [&]() {
        mcc::impl::MccSkyPoint hadec;
        for (size_t i = 0; i < N; ++i) {
            auto err = pcm.computePCM(mcc::impl::MccGenXY{(double)encX[i], (double)encY[i]}, &pcm_res, &hadec);
            if (err) {
                std::println("Error in PCM computing: {}", err.message());
                return 1;
            }

            pcmX[i] = pcm_res.pcmX;
            pcmY[i] = pcm_res.pcmY;

            ha[i] = hadec.co_lon();
            dec[i] = hadec.co_lat();
        }
        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 x_step = (stop_encX - start_encX) / (N - 1), y_step = (stop_encY - start_encY) / (N - 1);

    for (size_t i = 0; i < N; ++i) {
        encX[i] = start_encX + i * x_step;
        encY[i] = start_encY + i * y_step;
    }

    comp_hadec();

    std::println("X Y HA DEC (pcmX pcmY):");
    for (size_t i = 0; i < N; ++i) {
        std::println("{:12.7f} {:12.7f} {:12.7f} {:12.7f} ({:10.7f} {:10.7f})", encX[i].degrees(), encY[i].degrees(),
                     ha[i].degrees(), dec[i].degrees(), pcmX[i].degrees(), pcmY[i].degrees());
    }


    // add 'noise'
    std::random_device rd{};
    std::mt19937 gen{rd()};
    std::normal_distribution dist{0.0, sigma};

    mcc::impl::MccSkyPoint hadec;

    for (size_t i = 0; i < N; ++i) {
        ha[i] += dist(gen);
        dec[i] += dist(gen);

        resX[i] = ha[i] - encX[i];
        resY[i] = dec[i] - encY[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_t::pcm_data_t fit_pcm_data{.type = pcm_data.type, .siteLatitude = pcm_data.siteLatitude};

    auto r = pcm_const.computeModel(fit_pcm_data);
    if (r.error) {
        std::println("error: {}", r.error.message());
        return 1;
    }

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

    return 0;
}