#pragma once

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

#include <chrono>
#include <fstream>
#include "mcc_coord.h"

#ifdef VEC_XSIMD
#include <xsimd/xsimd.hpp>
#endif

#include "mcc_traits.h"
#include "mount_astrom_default.h"
#include "utils.h"

namespace erfa
{
#include <erfa.h>
#include <erfam.h>
}  // namespace erfa

namespace mcc::traits
{

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

template <typename T>
concept mcc_real_scalar_or_real_range_c =
    std::floating_point<T> || std::ranges::output_range<T, float> || std::ranges::output_range<T, double>;

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::astrom
{

// a time duration represented in radians (the precision is about 100 nanosecond for 10E8 secs (~1157 days))
typedef std::chrono::duration<double, std::ratio<314159265358979, 4320000000000000000>> mcc_chrono_radians;


// https://gssc.esa.int/navipedia/index.php?title=CEP_to_ITRF
static constexpr double mcc_UT1_to_sideral_ratio = 1.002737909350795;  // UT1/sideral


// modified Julian date (based on ERFA eraCal2jd)
template <traits::mcc_real_scalar_or_real_range_c ResT, traits::mcc_time_duration_c DT = std::chrono::milliseconds>
static int mcc_julday(traits::mcc_systime_c auto const& start_time,
                      ResT& mjd,
                      const DT& step = std::chrono::milliseconds(100))
    requires(!std::is_pointer_v<std::decay_t<ResT>>)
{
    size_t mjd_size = 0;
    if constexpr (std::ranges::range<ResT>) {
        mjd_size = std::ranges::distance(mjd.begin(), mjd.end());
        if (!mjd_size) {
            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;
    }

    int64_t im = (unsigned)ymd.month();
    int64_t id = (unsigned)ymd.day();
    int64_t iy = (int)ymd.year();

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

    // integer part of result MJD
    int64_t mjd_int = (1461LL * (iypmy + 4800LL)) / 4LL + (367LL * (im - 2LL - 12LL * my)) / 12LL -
                      (3LL * ((iypmy + 4900LL) / 100LL)) / 4LL + id - 2432076LL;

    using fd_t = std::chrono::duration<double, std::ratio<86400>>;  // fraction of day
    double mjd_float = static_cast<double>(mjd_int) + std::chrono::duration_cast<fd_t>(start_time - dd).count();

    if constexpr (std::ranges::range<ResT>) {
        double d_step = std::chrono::duration_cast<fd_t>(step).count();

        size_t i = 0;

#ifdef VEC_XSIMD
        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>{});

        alignas(xsimd::batch<double>::arch_type::alignment()) double arr[reg_size];

        auto ptr = mjd.begin();

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

            if constexpr (std::ranges::contiguous_range<ResT>) {
                res_reg.store_unaligned(mjd.data() + i);
                // res_reg.store_aligned(mjd.data() + i);
            } else {
                res_reg.store_aligned(arr);
                std::ranges::copy(arr, ptr);

                if constexpr (std::ranges::random_access_range<ResT>) {
                    ptr += reg_size;
                } else {
                    for (size_t k = 0; k < reg_size; ++k) {
                        ++ptr;
                    }
                }
            }
        }
#endif
        // scalar part
        for (size_t j = i; j < mjd_size; ++j, ++ptr) {
            *ptr = mjd_float + j * d_step;
        }
    } else {  // result is scalar
        mjd = mjd_float;
    }


    return 0;
}

template <typename ResT, traits::mcc_time_duration_c DT = std::chrono::milliseconds>
static int mcc_julday(traits::mcc_systime_c auto const& start_time,
                      ResT& mjd,
                      const DT& step = std::chrono::milliseconds(100))
    requires std::is_pointer_v<std::decay_t<ResT>>
{
    auto sp = std::span(mjd);
    int ret = mcc_julday(start_time, sp, step);
    return ret;
}


/*
 *    Computes a time duration to reach given altitude
 *
 *    the function returns a std::pair:
 *        .first = time duration to reach given altitude BEFORE object upper culmination
 *        .second = time duration to reach given altitude AFTER object upper culmination
 */
std::pair<std::chrono::duration<double>, std::chrono::duration<double>> mcc_time_to_alt(
    const MccAngle& alt,
    const MccAngle& ra,
    const MccAngle& dec,
    const MccAngle& lat,
    const MccAngle& lon,
    traits::mcc_systime_c auto const& now,
    traits::mcc_time_duration_c auto const& dut1,    // UT1-UTC
    traits::mcc_time_duration_c auto const& tt_tai,  // TT-TAI
    // TAI-UTC (leap seconds)
    traits::mcc_time_duration_c auto const& tai_utc)
{
    auto nan_dur = std::chrono::duration<double>(std::numeric_limits<double>::quiet_NaN());
    auto inf_dur = std::chrono::duration<double>(std::numeric_limits<double>::infinity());

    if (alt < 0.0) {
        return {nan_dur, nan_dur};
    }

    if (lat >= 0.0) {                                // north hemisphere
        if (dec < (lat - std::numbers::pi / 2.0)) {  // never rises above horizon
            return {nan_dur, nan_dur};
        }
    } else {                                         // south hemisphere
        if (dec > (lat + std::numbers::pi / 2.0)) {  // never rises above horizon
            return {nan_dur, nan_dur};
        }
    }

    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) {  // it never reach given altitude
        return {inf_dur, inf_dur};
    }


    auto ut1 = now + dut1;
    auto tt = now + tai_utc + tt_tai;

    double ut1_mjd, tt_mjd;
    int ret = mcc_julday(ut1, ut1_mjd);
    if (ret) {
        return {nan_dur, nan_dur};
    }
    ret = mcc_julday(tt, tt_mjd);
    if (ret) {
        return {nan_dur, nan_dur};
    }

    double lst_now = erfa::eraGst06a(ERFA_DJM0, ut1_mjd, ERFA_DJM0, tt_mjd);
    lst_now += lon;

    auto acos_ha = std::acos(cos_ha);

    // before upper culmination
    double lst_before = -acos_ha + ra;
    // after upper culmination
    double lst_after = acos_ha + ra;

    double d1 = (lst_before - lst_now) * mcc_UT1_to_sideral_ratio,
           d2 = (lst_after - lst_now) * mcc_UT1_to_sideral_ratio;

    if (d1 < 0.0) {  // the next day
        d1 += 2.0 * std::numbers::pi;
    }
    if (d2 < 0.0) {  // the next day
        d2 += 2.0 * std::numbers::pi;
    }

    static constexpr double rad2secs = 12.0 / std::numbers::pi * 3600.0;  // radians to time seconds

    return {std::chrono::duration<double>(d1 * rad2secs), std::chrono::duration<double>(d2 * rad2secs)};

    // return std::chrono::duration<double>(result * 12.0 / std::numbers::pi * 3600.0);  // in seconds
}

/*
 *  angles are in degrees or sexagimal string form
 *  RA and DEC are apparent!
 *
 *  returns
 *      NaN if object is non-rising or "alt_limit" < 0, Inf is circumpolar
 */

/*
double mcc_time_to_alt_limit(traits::mcc_real_or_char_range auto const& alt_limit,
                             traits::mcc_real_or_char_range auto const& RA,
                             traits::mcc_real_or_char_range auto const& DEC,
                             traits::mcc_real_or_char_range auto const& LAT,
                             traits::mcc_real_or_char_range auto const& LON,
                             traits::mcc_systime_c auto const& now,
                             traits::mcc_time_duration_c auto const& dut1,    // UT1-UTC
                             traits::mcc_time_duration_c auto const& tt_tai,  // TT-TAI
                             // TAI-UTC (leap seconds)
                             traits::mcc_time_duration_c auto const& tai_utc)
{
    //     sin(alt) = sin(DEC)*sin(phi) + cos(DEC)*cos(phi)*cos(HA)
    //     HA = LST - RA
    //     cos(HA) = cos(LST)*cos(RA) + sin(LST)*sin(RA)

    // using AT = std::decay_t<decltype(alt_limit)>;
    // using RT = std::decay_t<decltype(RA)>;
    // using DT = std::decay_t<decltype(DEC)>;
    // using LT = std::decay_t<decltype(LAT)>;
    // using LGT = std::decay_t<decltype(LON)>;

    double ra, dec, lat, lon, alt;

    auto to_rads = [](const auto& v, bool hms = false) {
        // using v_t = std::remove_cvref<decltype(v)>;
        using v_t = std::remove_cvref_t<decltype(v)>;
        double res;
        if constexpr (!std::floating_point<v_t>) {
            res = utils::parsAngleString(v, hms).value_or(std::numeric_limits<double>::quiet_NaN());
        } else {
            res = v;
        }

        if (!std::isfinite(res)) {
            return res;
        }

        return res * utils::deg2radCoeff;
    };

    alt = to_rads(alt_limit);
    if (!std::isfinite(alt)) {
        return alt;
    }

    if (alt < 0.0) {
        return std::numeric_limits<double>::quiet_NaN();
    }

    ra = to_rads(RA, true);
    if (!std::isfinite(ra)) {
        return ra;
    }

    dec = to_rads(DEC);
    if (!std::isfinite(dec)) {
        return dec;
    }

    lat = to_rads(LAT);
    if (!std::isfinite(lat)) {
        return lat;
    }

    lon = to_rads(LON);
    if (!std::isfinite(lon)) {
        return lon;
    }

    if (lat >= 0.0) {                                // north hemisphere
        if (dec < (lat - std::numbers::pi / 2.0)) {  // never rises above horizon
            return std::numeric_limits<double>::quiet_NaN();
        }
    } else {                                         // south hemisphere
        if (dec > (lat + std::numbers::pi / 2.0)) {  // never rises above horizon
            return std::numeric_limits<double>::quiet_NaN();
        }
    }

    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) {  // it never sets below given altitude
        return std::numeric_limits<double>::infinity();
    }

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

    auto ut1 = now + dut1;
    auto tt = now + tai_utc + tt_tai;

    double ut1_mjd, tt_mjd, result;
    int ret = mcc_julday(ut1, ut1_mjd);
    if (ret) {
        return std::numeric_limits<double>::quiet_NaN();
    }
    ret = mcc_julday(tt, tt_mjd);
    if (ret) {
        return std::numeric_limits<double>::quiet_NaN();
    }

    double lst_now = erfa::eraGst06a(ERFA_DJM0, ut1_mjd, ERFA_DJM0, tt_mjd);
    lst_now += lon;

    result = lst - lst_now;
    if (result < 0.0) {  // the next sideral day
        result += 2.0 * std::numbers::pi;
    }

    if (result > std::numbers::pi) {  // object is already below the limit
        return 0.0;
    }

    result *= mcc_UT1_to_sideral_ratio;  // to UT1 scale

    return result;
}

*/

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


class MccLeapSeconds final
{
public:
    typedef std::chrono::system_clock::time_point time_point_t;
    typedef std::chrono::duration<double> real_secs_t;  // seconds duration in double

    MccLeapSeconds()
    {
        // create default values
        std::istringstream ist(defaults::MCC_DEFAULT_LEAP_SECONDS_FILE);

        load(ist);
    }

    ~MccLeapSeconds() = default;


    time_point_t expireDate() const
    {
        return _expireDate;
    }

    // load from stream
    bool load(std::derived_from<std::basic_istream<char>> auto& stream, char comment_sym = '#')
    {
        std::istringstream is;
        double mjd;
        unsigned day, month;
        int year;
        double tai_utc;

        decltype(_expireDate) edate;
        std::vector<leapsecond_db_elem_t> db;

        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
                    if (std::regex_match(line, expr_date_rx)) {
                        auto pos = line.find("on");
                        sv = utils::trimSpaces(std::string_view{line.begin() + pos + 2, line.end()},
                                               utils::TrimType::TRIM_LEFT);
                        is.str({sv.begin(), sv.end()});
                        is >> std::chrono::parse("%d %B %Y", edate);
                        is.clear();
                    }
                    continue;
                }
            } else {
                continue;
            }

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

                continue;
            }
        }

        if (db.empty()) {  // keep previous data
            return false;
        }

        _expireDate = std::move(edate);
        _db = std::move(db);

        return true;
    }


    bool load(traits::mcc_input_char_range auto const& filename, char comment_sym = '#')
    {
        std::ifstream fst(filename);

        bool ok = fst.is_open();
        if (!ok) {
            return false;
        }

        ok = load(fst, comment_sym);

        fst.close();

        return ok;
    }


    // std::optional<double> operator[](const time_point_t& tp) const
    std::optional<real_secs_t> operator[](const time_point_t& tp) const
    {
        if (tp > _expireDate) {  // ???????!!!!!!!!!!!
            return std::nullopt;
            // return _db.back().tai_utc;
        }

        std::chrono::year_month_day ymd{std::chrono::floor<std::chrono::days>(tp)};

        for (auto const& el : _db | std::views::reverse) {
            if (ymd >= el.ymd) {
                // return el.tai_utc;
                return real_secs_t{el.tai_utc};
            }
        }

        return std::nullopt;
    }

    // std::optional<double> operator[](const double& mjd) const
    std::optional<real_secs_t> operator[](const double& mjd) const
    {
        double e_mjd;

        astrom::mcc_julday(_expireDate, e_mjd);
        if (mjd > e_mjd) {  // ???????!!!!!!!!!!!
            return std::nullopt;
            // return _db.back().tai_utc;
        }

        for (auto const& el : _db | std::views::reverse) {
            if (mjd >= el.mjd) {
                return real_secs_t{el.tai_utc};
            }
        }

        return std::nullopt;
    }


    void dump(std::derived_from<std::basic_ostream<char>> auto& stream) const
    {
        stream << std::format("Leap seconds database expire date: {}", _expireDate) << '\n';
        for (auto const& el : _db) {
            stream << std::format("{} {} {}", el.mjd, el.ymd, el.tai_utc) << '\n';
        }
    }


private:
    inline static 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} *$"};

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

    time_point_t _expireDate{};

    struct leapsecond_db_elem_t {
        double mjd;
        std::chrono::year_month_day ymd;
        double tai_utc;  // TAI-UTC in seconds
    };

    std::vector<leapsecond_db_elem_t> _db{};
};



class MccIersBulletinA final
{
public:
    typedef std::chrono::system_clock::time_point time_point_t;
    typedef std::chrono::duration<double> real_secs_t;  // seconds duration in double

    struct pole_pos_t {
        double x, y;
    };

    struct date_range_t {
        std::chrono::year_month_day begin;
        std::chrono::year_month_day end;
    };

    struct date_range_mjd_t {
        double begin;
        double end;
    };

    MccIersBulletinA()
    {
        // create pre-defined (default-state) database
        std::istringstream ist(defaults::MCC_DEFAULT_IERS_BULLETIN_A_FILE);

        load(ist);
    }

    ~MccIersBulletinA() = default;

    std::chrono::system_clock::time_point bulletinDate() const
    {
        return _date;
    }


    date_range_t dateRange() const
    {
        return {_db.front().ymd, _db.back().ymd};
    }

    date_range_mjd_t dateRangeMJD() const
    {
        return {_db.front().mjd, _db.back().mjd};
    }


    // double TT_TAI() const
    real_secs_t TT_TAI() const
    {
        return real_secs_t{_tt_tai};
    }

    // DUT1 = UT1 - UTC
    // std::optional<double> DUT1(const time_point_t& tp) const
    std::optional<real_secs_t> DUT1(const time_point_t& tp) const
    {
        // use of the closest date
        std::chrono::year_month_day ymd{std::chrono::round<std::chrono::days>(tp)};

        if (ymd < _db.front().ymd) {
            return std::nullopt;
        }

        if (ymd > _db.back().ymd) {
            return std::nullopt;
        }

        for (auto const& el : _db) {
            if (ymd <= el.ymd) {
                return real_secs_t{el.dut1};
            }
        }

        return std::nullopt;
    }

    // std::optional<double> DUT1(double mjd) const
    std::optional<real_secs_t> DUT1(double mjd) const
    {
        mjd = std::round(mjd);  // round to closest integer MJD

        if (mjd < _db.front().mjd) {
            return std::nullopt;
        }

        if (mjd > _db.back().mjd) {
            return std::nullopt;
        }

        for (auto const& el : _db) {
            if (mjd <= el.mjd) {
                return real_secs_t{el.dut1};
            }
        }

        return std::nullopt;
    }

    std::optional<pole_pos_t> polarCoords(const time_point_t& tp) const
    {
        std::chrono::year_month_day ymd{std::chrono::round<std::chrono::days>(tp)};

        if (ymd < _db.front().ymd) {
            return std::nullopt;
        }

        if (ymd > _db.back().ymd) {
            return std::nullopt;
        }

        for (auto const& el : _db) {
            if (ymd <= el.ymd) {
                return pole_pos_t{el.x, el.y};
            }
        }

        return std::nullopt;
    }

    std::optional<pole_pos_t> polarCoords(double mjd) const
    {
        mjd = std::round(mjd);  // round to closest integer MJD

        if (mjd < _db.front().mjd) {
            return std::nullopt;
        }

        if (mjd > _db.back().mjd) {
            return std::nullopt;
        }

        for (auto const& el : _db) {
            if (mjd <= el.mjd) {
                return pole_pos_t{el.x, el.y};
            }
        }

        return std::nullopt;
    }

    bool load(std::derived_from<std::basic_istream<char>> auto& stream, char comment_sym = '*')
    {
        std::vector<earth_orient_db_elem_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;
        decltype(_date) bdate;
        double tt_tai;

        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.emplace_back(mjd, std::chrono::year_month_day{std::chrono::year{year} / month / day}, 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", bdate);
                    continue;
                }

                if (std::regex_match(sv.begin(), sv.end(), bull_tt_tai_rx)) {
                    is.str({sv.begin(), sv.end()});
                    std::string dummy;
                    is >> dummy >> dummy >> dummy >> dummy >> tt_tai;
                    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.empty()) {
            return false;
        }

        _date = std::move(bdate);
        _tt_tai = tt_tai;
        _db = std::move(db);

        return true;
    }

    bool load(traits::mcc_input_char_range auto const& filename, char comment_sym = '*')
    {
        std::ifstream fst(filename);

        bool ok = fst.is_open();
        if (!ok) {
            return false;
        }

        ok = load(fst, comment_sym);

        fst.close();

        return ok;
    }


    void dump(std::derived_from<std::basic_ostream<char>> auto& stream) const
    {
        stream << std::format("Bulletin A issue date: {}", _date) << '\n';
        stream << std::format("TT-TAI: {}", _tt_tai) << '\n';
        for (auto const& el : _db) {
            stream << std::format("{} {} {:6.4f} {:6.4f} {:7.5f}", el.mjd, el.ymd, el.x, el.y, el.dut1) << '\n';
        }
    }

private:
    inline static 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]+ *$"};

    inline static const std::regex bull_tt_tai_rx{"^ *TT += +TAI +\\+ +[0-9]+\\.[0-9]+ +seconds *$"};

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

    inline static 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]+ *$"};


    time_point_t _date;
    double _tt_tai;

    struct earth_orient_db_elem_t {
        double mjd;
        std::chrono::year_month_day ymd;
        double x, y;  // Polar coordinates in arcsecs
        double dut1;  // UT1-UTC in seconds
    };

    std::vector<earth_orient_db_elem_t> _db;
};


}  // namespace mcc::astrom