#pragma once

/*********************************
 *    MOUNT CONTROL COMPONENTS    *
 *                                *
 *      astrometry functions      *
 *********************************/

#include <chrono>
#include <unordered_map>
#include <xsimd/xsimd.hpp>
#include "utils.h"


namespace mcc::traits
{

template <typename T>
concept mcc_scalar_or_simd_c = xsimd::is_batch<T>::value || std::is_arithmetic_v<T>;


template <typename T>
concept mcc_real_or_char_range =
    std::floating_point<T> ||
    (std::ranges::contiguous_range<T> && std::same_as<char, std::remove_cvref_t<std::ranges::range_value_t<T>>>);

}  // namespace mcc::traits


namespace mcc::astro
{

// modified Julian date (based on ERFA eraCal2jd)
// template <traits::mcc_scalar_or_simd_c T>
static int mcc_julday(const std::chrono::system_clock::time_point& start_time,
                      const std::chrono::system_clock::duration& step,
                      std::vector<double, xsimd::default_allocator<double>>& mjd)
{
    if (mjd.empty()) {
        return -100;
    }

    using namespace std::literals::chrono_literals;

    auto dd = std::chrono::floor<std::chrono::days>(start_time);
    std::chrono::year_month_day ymd{dd};

    static constexpr std::chrono::year MIN_YEAR = -4799y;

    if (ymd.year() < MIN_YEAR) {
        return -1;
    }
    if (!ymd.month().ok()) {
        return -2;
    }

    // my = (im - 14) / 12;
    // iypmy = (long) (iy + my);

    int64_t my = -(14 - (unsigned)ymd.month()) / 12;
    int64_t iypmy = (int)ymd.year() + my;

    // (1461L * (iypmy + 4800L)) / 4L
    //     + (367L * (long) (im - 2 - 12 * my)) / 12L
    //     - (3L * ((iypmy + 4900L) / 100L)) / 4L
    //     + (long) id - 2432076L

    // integer part of result MJD
    int64_t mjd_int = 1461LL * (iypmy + 480LL) / 4LL +
                      (367LL * ((int64_t)(unsigned)ymd.month() - 2LL - 12LL * my)) / 12LL -
                      (3LL * (iypmy + 4900LL) / 100LL) / 4LL + (int64_t)(unsigned)ymd.day() - 2432076LL;

    constexpr double nanosec = 1.0 / 24.0 / 3600.0 / 1.0E-9;  // 1 nanosecond in days

    double mjd_float = static_cast<double>(mjd_int) +
                       std::chrono::duration_cast<std::chrono::nanoseconds>(start_time - dd).count() * nanosec;
    double d_step = std::chrono::duration_cast<std::chrono::nanoseconds>(step).count() * nanosec;

    constexpr size_t reg_size = xsimd::batch<double>::size;
    size_t vec_size = mjd.size() - mjd.size() % reg_size;

    xsimd::batch<double> res_reg{mjd_float};
    xsimd::batch<double> step_reg = [d_step]<size_t... Is>(std::index_sequence<Is...>) {
        return xsimd::batch<double>{(Is * d_step)...};
    }(std::make_index_sequence<reg_size>{});

    // vectorized part
    size_t i = 0;
    for (; i < vec_size; i += vec_size) {
        res_reg += step_reg;

        res_reg.store_aligned(mjd.data() + i);
    }

    // scalar part
    for (size_t j = i; j < mjd.size(); ++j) {
        mjd[j] = mjd_float + j * d_step;
    }


    return 0;
}

/*
 *  angles are in degrees or sexagimal string form
 */
template <traits::mcc_real_or_char_range AT,
          traits::mcc_real_or_char_range RT,
          traits::mcc_real_or_char_range DT,
          traits::mcc_real_or_char_range LT>
double mcc_time_to_alt_limit(const AT& alt_limit,
                             const RT& RA,
                             const DT& DEC,
                             const LT& LAT,
                             const std::chrono::system_clock::time_point& now)
{
    //     sin(alt) = sin(DEC)*sin(phi) + cos(DEC)*cos(phi)*cos(HA)
    //     HA = LST - RA

    double ra, dec, lat, alt;

    if constexpr (std::floating_point<AT>) {
        alt = alt_limit * utils::deg2radCoeff;
    } else {
        alt = utils::parsAngleString(alt_limit);
        alt *= utils::deg2radCoeff;
    }

    if constexpr (std::floating_point<RT>) {
        ra = RA * utils::deg2radCoeff;
    } else {
        ra = utils::parsAngleString(RA, true);
        ra *= utils::deg2radCoeff;
    }

    if constexpr (std::floating_point<DT>) {
        dec = DEC * utils::deg2radCoeff;
    } else {
        dec = utils::parsAngleString(DEC);
        dec *= utils::deg2radCoeff;
    }

    if constexpr (std::floating_point<LT>) {
        lat = LAT * utils::deg2radCoeff;
    } else {
        lat = utils::parsAngleString(LAT);
        lat *= utils::deg2radCoeff;
    }

    double cos_ha = (std::sin(alt) - std::sin(dec) * std::sin(lat)) / std::cos(dec) / std::cos(lat);
    if (std::abs(cos_ha) > 1.0) {  // circumpolar (it never sets below horizon)
        return std::numeric_limits<double>::infinity();
    }

    double lst = std::acos(cos_ha) + ra;


    return 0.0;
}


/* IERS bulletins
 *
 * BULLETIN A: https://datacenter.iers.org/data/latestVersion/bulletinA.txt
 * leapseconds: https://hpiers.obspm.fr/iers/bul/bulc/Leap_Second.dat
 *
 */

struct leapsecond_db_elem_t {
    double mjd;
    unsigned day, month;
    int year;
    // std::chrono::year_month_day ymd;
    double tai_utc;  // TAI-UTC in seconds
};

typedef std::vector<leapsecond_db_elem_t> leapsecond_db_t;


struct earth_orient_db_elem_t {
    int year;
    unsigned month, day;
    double mjd;
    double x, y;  // Polar coordinates in arcsecs
    double dut1;  // UT1-UTC in seconds
};


// typedef std::vector<earth_orient_db_elem_t> earth_orient_db_t;
struct earth_orient_db_t {
    std::chrono::system_clock::time_point bulletinDate{};
    std::vector<earth_orient_db_elem_t> db{};
};

// init to some known state
static leapsecond_db_t CURRENT_LEAPSECONDS_DB = {
    {41317.0, 1, 1, 1972, 10}, {41499.0, 1, 7, 1972, 11}, {41683.0, 1, 1, 1973, 12},
    {42048.0, 1, 1, 1974, 13}, {42413.0, 1, 1, 1975, 14}, {42778.0, 1, 1, 1976, 15},
    {43144.0, 1, 1, 1977, 16}, {43509.0, 1, 1, 1978, 17}, {43874.0, 1, 1, 1979, 18},
    {44239.0, 1, 1, 1980, 19}, {44786.0, 1, 7, 1981, 20}, {45151.0, 1, 7, 1982, 21},
    {45516.0, 1, 7, 1983, 22}, {46247.0, 1, 7, 1985, 23}, {47161.0, 1, 1, 1988, 24},
    {47892.0, 1, 1, 1990, 25}, {48257.0, 1, 1, 1991, 26}, {48804.0, 1, 7, 1992, 27},
    {49169.0, 1, 7, 1993, 28}, {49534.0, 1, 7, 1994, 29}, {50083.0, 1, 1, 1996, 30},
    {50630.0, 1, 7, 1997, 31}, {51179.0, 1, 1, 1999, 32}, {53736.0, 1, 1, 2006, 33},
    {54832.0, 1, 1, 2009, 34}, {56109.0, 1, 7, 2012, 35}, {57204.0, 1, 7, 2015, 36},
    {57754.0, 1, 1, 2017, 37}

};


static earth_orient_db_t CURRENT_EARTH_ORIENT_DB;

static bool mcc_parse_bulletinA(std::derived_from<std::basic_istream<char>> auto& stream, char comment_sym = '*')
{
    const std::regex bull_date_rx{
        "^ *[0-9]{1,2} +(January|February|March|April|May|June|July|August|September|October|November|December) "
        "+[0-9]{4,} +Vol\\. +[XMLCDVI]+ +No\\. +[0-9]+ *$"};

    const std::regex bull_tab_title_rx{"^ *MJD +x\\(arcsec\\) +y\\(arcsec\\) +UT1-UTC\\(sec\\) *$"};

    // 2025  3  7  60741       0.0663      0.3341      0.04348
    const std::regex bull_tab_vals_rx{
        "^ *[0-9]{4,} +[0-9]{1,2} +[0-9]{1,2} +[0-9]{5,} +[0-9]+\\.[0-9]+ +[0-9]+\\.[0-9]+ +[0-9]+\\.[0-9]+ *$"};

    earth_orient_db_t db;
    enum { TAB_STATE_SEEK, TAB_STATE_START };
    int tab_state = TAB_STATE_SEEK;

    int year;
    unsigned month, day;
    double mjd, x, y, dut1;
    std::istringstream is;

    for (std::string line; std::getline(stream, line);) {
        if (line.empty()) {
            continue;
        }

        auto sv = utils::trimSpaces(line, utils::TrimType::TRIM_LEFT);

        if (sv.size()) {
            if (sv[0] == comment_sym) {  // comment string
                continue;
            }

            if (tab_state == TAB_STATE_START) {
                if (std::regex_match(sv.begin(), sv.end(), bull_tab_vals_rx)) {
                    // is.str({sv.begin(), sv.end()});
                    is.str(line);
                    is >> year >> month >> day >> mjd >> x >> y >> dut1;
                    db.db.emplace_back(year, month, day, mjd, x, y, dut1);
                    is.clear();
                } else {  // end of the table - just stop parsing
                    break;
                }

                continue;
            }

            if (std::regex_match(sv.begin(), sv.end(), bull_date_rx)) {
                is.str({sv.begin(), sv.end()});
                is >> std::chrono::parse("%d %B %Y", db.bulletinDate);
                continue;
            }

            if (std::regex_match(sv.begin(), sv.end(), bull_tab_title_rx)) {
                tab_state = TAB_STATE_START;
                continue;
            }

        } else {  // empty string (only spaces)
            continue;
        }
    }

    if (db.db.empty()) {
        return false;
    }

    CURRENT_EARTH_ORIENT_DB = std::move(db);

    return true;
}


static bool mcc_parse_leapsecs(std::derived_from<std::basic_istream<char>> auto& stream, char comment_sym = '#')
{
    // #  File expires on 28 December 2025
    const std::regex expr_date_rx{
        "^ *# *File +expires +on +[0-8]{1,2} "
        "+(January|February|March|April|May|June|July|August|September|October|November|December) +[0-9]{4} *$"};

    const std::regex data_rx{"^ *[0-9]{5,}(\\.?[0-9]+) +[0-9]{1,2} +[0-9]{1,2} +[0-9]{4} +[0-9]{1,} *$"};

    std::istringstream is;
    double mjd;
    unsigned day, month;
    int year;
    double tai_utc;


    leapsecond_db_t db;
    // leapsecond_db_elem_t db_elem;

    for (std::string line; std::getline(stream, line);) {
        auto sv = utils::trimSpaces(line, utils::TrimType::TRIM_LEFT);

        if (sv.size()) {
            if (sv[0] == comment_sym) {  // comment string
                continue;
            }
        } else {
            continue;
        }

        if (std::regex_match(line, data_rx)) {
            is.str(std::move(line));
            is >> mjd >> day >> month >> year >> tai_utc;
            db.emplace_back(mjd, day, month, year, tai_utc);
            is.clear();
        }
    }

    if (db.empty()) {
        return false;
    }

    CURRENT_LEAPSECONDS_DB = std::move(db);

    return true;
}



}  // namespace mcc::astro