#pragma once

#include <erfa.h>
#include <erfam.h>
#include <mutex>

#include "mcc_ccte_iers.h"
// #include "mcc_coord.h"
#include "mcc_generics.h"

namespace mcc::ccte::erfa
{


enum class MccCCTE_ERFAErrorCode : int {
    ERROR_OK = 0,
    ERROR_NULLPTR,
    ERROR_INVALID_INPUT_ARG,
    ERROR_julday_INVALID_YEAR,
    ERROR_julday_INVALID_MONTH,
    ERROR_julday_INVALID_DAY,
    ERROR_UNSUPPORTED_COORD_PAIR,
    ERROR_BULLETINA_OUT_OF_RANGE,
    ERROR_LEAPSECONDS_OUT_OF_RANGE,
    ERROR_DUBIOUS_YEAR,
    ERROR_UNACCEPTABLE_DATE,
    ERROR_UPDATE_LEAPSECONDS,
    ERROR_UPDATE_BULLETINA,
    ERROR_UNEXPECTED
};

}  // namespace mcc::ccte::erfa


namespace std
{

template <>
class is_error_code_enum<mcc::ccte::erfa::MccCCTE_ERFAErrorCode> : public true_type
{
};

}  // namespace std



namespace mcc::ccte::erfa
{

/*  error category definition  */

// error category
struct MccCCTE_ERFACategory : public std::error_category {
    MccCCTE_ERFACategory() : std::error_category() {}

    const char* name() const noexcept
    {
        return "CCTE-ERFA";
    }

    std::string message(int ec) const
    {
        MccCCTE_ERFAErrorCode err = static_cast<MccCCTE_ERFAErrorCode>(ec);

        switch (err) {
            case MccCCTE_ERFAErrorCode::ERROR_OK:
                return "OK";
            case MccCCTE_ERFAErrorCode::ERROR_NULLPTR:
                return "input argument is the nullptr";
            case MccCCTE_ERFAErrorCode::ERROR_INVALID_INPUT_ARG:
                return "invalid argument";
            case MccCCTE_ERFAErrorCode::ERROR_julday_INVALID_YEAR:
                return "invalid year number";
            case MccCCTE_ERFAErrorCode::ERROR_julday_INVALID_MONTH:
                return "invalid month number";
            case MccCCTE_ERFAErrorCode::ERROR_julday_INVALID_DAY:
                return "invalid day number";
            case MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR:
                return "unsupported coordinate pair";
            case MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE:
                return "time point is out of range";
            case MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE:
                return "time point is out of range";
            case MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR:
                return "dubious year";
            case MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE:
                return "unacceptable date";
            case MccCCTE_ERFAErrorCode::ERROR_UPDATE_LEAPSECONDS:
                return "leap seconds update error";
            case MccCCTE_ERFAErrorCode::ERROR_UPDATE_BULLETINA:
                return "bulletin A update error";
            case MccCCTE_ERFAErrorCode::ERROR_UNEXPECTED:
                return "unexpected error value";
            default:
                return "UNKNOWN";
        }
    }

    static const MccCCTE_ERFACategory& get()
    {
        static const MccCCTE_ERFACategory constInst;
        return constInst;
    }
};


inline std::error_code make_error_code(MccCCTE_ERFAErrorCode ec)
{
    return std::error_code(static_cast<int>(ec), MccCCTE_ERFACategory::get());
}


class MccCCTE_ERFA
{
    static constexpr double PI_2 = std::numbers::pi / 2.0;

public:
    static constexpr double DEFAULT_WAVELENGTH = 0.55;  // default observed wavelength in mkm

    typedef std::error_code error_t;

    struct refract_model_t {
        static constexpr std::string_view name()
        {
            return "ERFA";
        }

        double refa, refb;
    };


    // meteo parameters (to compute refraction)
    struct meteo_t {
        typedef double temp_t;
        typedef double humid_t;
        typedef double press_t;

        temp_t temperature;  // Temperature in C
        humid_t humidity;    // humidity in % ([0.0, 1.0])
        press_t pressure;    // atmospheric presure in hPa=mB
    };

    // celestial object addition parameters
    struct obj_pars_t {
        double pm_RA = 0.0;   // rads/year
        double pm_DEC = 0.0;  // rads/year
        double parallax;      // in arcsecs
        double radvel;        // radial velocity (signed, km/s)
    };

    struct engine_state_t {
        meteo_t meteo{.temperature = 0.0, .humidity = 0.5, .pressure = 1010.0};

        double wavelength = DEFAULT_WAVELENGTH;  // observed wavelength in mkm

        double lat = 0.0;   // site latitude
        double lon = 0.0;   // site longitude
        double elev = 0.0;  // site elevation (in meters)

        mcc::ccte::iers::MccLeapSeconds _leapSeconds{};
        mcc::ccte::iers::MccIersBulletinA _bulletinA{};
    };

    MccCCTE_ERFA() : _stateMutex(new std::mutex) {}

    MccCCTE_ERFA(engine_state_t state) : _currentState(std::move(state)), _stateMutex(new std::mutex) {}

    MccCCTE_ERFA(const MccCCTE_ERFA&) = delete;
    MccCCTE_ERFA& operator=(const MccCCTE_ERFA&) = delete;

    MccCCTE_ERFA(MccCCTE_ERFA&&) = default;
    MccCCTE_ERFA& operator=(MccCCTE_ERFA&&) = default;

    virtual ~MccCCTE_ERFA() = default;


    std::string_view nameCCTE() const
    {
        return "ERFA-CCTE-ENGINE";
    }

    // engine state related methods

    void setStateERFA(engine_state_t state)
    {
        std::lock_guard lock{*_stateMutex};

        _currentState = std::move(state);
    }

    engine_state_t getStateERFA() const
    {
        std::lock_guard lock{*_stateMutex};

        return _currentState;
    }

    void updateMeteoERFA(meteo_t meteo)
    {
        std::lock_guard lock{*_stateMutex};

        _currentState.meteo = std::move(meteo);
    }

    error_t updateLeapSeconds(std::derived_from<std::basic_istream<char>> auto& stream, char comment_sym = '#')
    {
        std::lock_guard lock{*_stateMutex};

        if (!_currentState._leapSeconds.load(stream, comment_sym)) {
            return MccCCTE_ERFAErrorCode::ERROR_UPDATE_LEAPSECONDS;
        }

        return MccCCTE_ERFAErrorCode::ERROR_OK;
    }


    error_t updateLeapSeconds(traits::mcc_input_char_range auto const& filename, char comment_sym = '#')
    {
        std::lock_guard lock{*_stateMutex};

        if (!_currentState._leapSeconds.load(filename, comment_sym)) {
            return MccCCTE_ERFAErrorCode::ERROR_UPDATE_LEAPSECONDS;
        }

        return MccCCTE_ERFAErrorCode::ERROR_OK;
    }


    error_t updateBulletinA(std::derived_from<std::basic_istream<char>> auto& stream, char comment_sym = '*')
    {
        std::lock_guard lock{*_stateMutex};

        if (!_currentState._bulletinA.load(stream, comment_sym)) {
            return MccCCTE_ERFAErrorCode::ERROR_UPDATE_BULLETINA;
        }

        return MccCCTE_ERFAErrorCode::ERROR_OK;
    }


    error_t updateBulletinA(traits::mcc_input_char_range auto const& filename, char comment_sym = '*')
    {
        std::lock_guard lock{*_stateMutex};

        if (!_currentState._bulletinA.load(filename, comment_sym)) {
            return MccCCTE_ERFAErrorCode::ERROR_UPDATE_BULLETINA;
        }

        return MccCCTE_ERFAErrorCode::ERROR_OK;
    }


    // apparent sideral time (Greenwitch or local)
    error_t apparentSideralTime(mcc_coord_epoch_c auto const& epoch, mcc_angle_c auto* st, bool islocal = false)
    {
        error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;

        if (st == nullptr) {
            return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
        }

        using real_days_t = std::chrono::duration<double, std::ratio<86400>>;

        double ut1 = epoch.MJD();
        double tt = epoch.MJD();

        std::lock_guard lock{*_stateMutex};

        auto dut1 = _currentState._bulletinA.DUT1(epoch.MJD());

        if (dut1.has_value()) {
            ut1 += std::chrono::duration_cast<real_days_t>(dut1.value()).count();
        } else {  // out of range
            return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
        }

        auto tai_utc = _currentState._leapSeconds[epoch.MJD()];
        if (tai_utc.has_value()) {
            tt += std::chrono::duration_cast<real_days_t>(tai_utc.value()).count();
        } else {
            return MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE;
        }


        auto tt_tai = _currentState._bulletinA.TT_TAI();
        tt += std::chrono::duration_cast<real_days_t>(tt_tai).count();

        *st = eraGst06a(ERFA_DJM0, ut1, ERFA_DJM0, tt);

        if (islocal) {
            *st = eraAnp(*st + _currentState.lon);
        }

        return ret;
    }


    // ICRS to observed
    // returned azimuth is counted from the South through the West
    error_t icrsToObs(mcc_angle_c auto const& ra_icrs,
                      mcc_angle_c auto const& dec_icrs,
                      mcc_coord_epoch_c auto const& epoch,
                      mcc_angle_c auto* ra_obs,
                      mcc_angle_c auto* dec_obs,
                      mcc_angle_c auto* ha_obs,
                      mcc_angle_c auto* az,
                      mcc_angle_c auto* zd,
                      obj_pars_t* obj_params = nullptr)
    {
        return icrsTo(true, ra_icrs, dec_icrs, epoch, ra_obs, dec_obs, ha_obs, az, zd, obj_params);
    }

    // error_t icrsToObs(MccSkyRADEC_ICRS const& radec_icrs,
    //                   MccSkyRADEC_OBS* radec_obs,
    //                   MccSkyAZZD* azzd,
    //                   mcc_angle_c auto* ha_obs,
    //                   obj_pars_t* obj_params = nullptr)
    // {
    //     double ra_obs, dec_obs, az, zd, ha;

    //     auto err =
    //         icrsToObs(radec_icrs.x(), radec_icrs.y(), radec_icrs.epoch(), &ra_obs, &dec_obs, &ha, &az, &zd,
    //         obj_params);

    //     if (!err) {
    //         if (radec_obs) {
    //             radec_obs->setX(ra_obs);
    //             radec_obs->setY(dec_obs);
    //         }

    //         if (azzd) {
    //             azzd->setEpoch(radec_obs->epoch());
    //             azzd->setX(az);
    //             azzd->setY(zd);
    //         }

    //         if (ha_obs) {
    //             *ha_obs = ha;
    //         }
    //     }

    //     return err;
    // };


    // ICRS to apparent (in vacuo)
    // returned azimuth is counted from the South through the West
    error_t icrsToApp(mcc_angle_c auto const& ra_icrs,
                      mcc_angle_c auto const& dec_icrs,
                      mcc_coord_epoch_c auto const& epoch,
                      mcc_angle_c auto* ra_app,
                      mcc_angle_c auto* dec_app,
                      mcc_angle_c auto* ha_app,
                      mcc_angle_c auto* az,
                      mcc_angle_c auto* zd,  // should be interpretated as zenithal distance corrected for refraction
                      obj_pars_t* obj_params = nullptr)
    {
        return icrsTo(false, ra_icrs, dec_icrs, epoch, ra_app, dec_app, ha_app, az, zd, obj_params);
    }


    // error_t icrsToApp(MccSkyRADEC_ICRS const& radec_icrs,
    //                   MccSkyRADEC_OBS* radec_app,
    //                   MccSkyAZZD* azzd,
    //                   mcc_angle_c auto* ha_app,
    //                   obj_pars_t* obj_params = nullptr)
    // {
    //     double ra_app, dec_app, az, zd, ha;

    //     auto err =
    //         icrsToApp(radec_icrs.x(), radec_icrs.y(), radec_icrs.epoch(), &ra_app, &dec_app, &ha, &az, &zd,
    //         obj_params);

    //     if (!err) {
    //         if (radec_app) {
    //             radec_app->setX(ra_app);
    //             radec_app->setY(dec_app);
    //         }

    //         if (azzd) {
    //             azzd->setEpoch(radec_app->epoch());
    //             azzd->setX(az);
    //             azzd->setY(zd);
    //         }

    //         if (ha_app) {
    //             *ha_app = ha;
    //         }
    //     }

    //     return err;
    // }


    error_t obsToICRS(MccCoordPairKind obs_type,
                      mcc_coord_epoch_c auto const& epoch,
                      mcc_angle_c auto const& co_lon,
                      mcc_angle_c auto const& co_lat,
                      mcc_angle_c auto* ra_icrs,
                      mcc_angle_c auto* dec_icrs)
    {
        return toICRS(true, obs_type, epoch, co_lon, co_lat, ra_icrs, dec_icrs);
    }


    // error_t obsToICRS(mcc_coord_pair_c auto const& xy_obs, MccSkyRADEC_ICRS* radec_icrs)
    // {
    //     double ra, dec;

    //     auto err = obsToICRS(xy_obs.pair_kind, xy_obs.epoch(), xy_obs.x(), xy_obs.y(), &ra, &dec);
    //     if (err) {
    //         return err;
    //     }

    //     if (radec_icrs) {
    //         radec_icrs->setX(ra);
    //         radec_icrs->setY(dec);
    //     }

    //     return err;
    // }


    error_t appToICRS(MccCoordPairKind app_type,
                      mcc_coord_epoch_c auto const& epoch,
                      mcc_angle_c auto const& co_lon,
                      mcc_angle_c auto const& co_lat,
                      mcc_angle_c auto* ra_icrs,
                      mcc_angle_c auto* dec_icrs)
    {
        return toICRS(false, app_type, epoch, co_lon, co_lat, ra_icrs, dec_icrs);
    }


    // error_t appToICRS(mcc_coord_pair_c auto const& xy_app, MccSkyRADEC_ICRS* radec_icrs)
    // {
    //     double ra, dec;

    //     auto err = appToICRS(xy_app.pair_kind, xy_app.epoch(), xy_app.x(), xy_app.y(), &ra, &dec);
    //     if (!err) {
    //         if (radec_icrs) {
    //             radec_icrs->setX(ra);
    //             radec_icrs->setY(dec);
    //         }
    //     }

    //     return err;
    // }


    error_t equationOrigins(mcc_coord_epoch_c auto const& epoch, mcc_angle_c auto* eo)
    {
        if (eo == nullptr) {
            return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
        }

        error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;

        std::lock_guard lock{*_stateMutex};

        using real_days_t = std::chrono::duration<double, std::ratio<86400>>;

        double mjd = epoch.MJD();

        auto tai_utc = _currentState._leapSeconds[mjd];
        if (tai_utc.has_value()) {
            double tt = mjd;
            tt += std::chrono::duration_cast<real_days_t>(tai_utc.value()).count();

            auto tt_tai = _currentState._bulletinA.TT_TAI();
            tt += +std::chrono::duration_cast<real_days_t>(tt_tai).count();

            *eo = eraEo06a(ERFA_DJM0, tt);
        } else {
            ret = MccCCTE_ERFAErrorCode::ERROR_LEAPSECONDS_OUT_OF_RANGE;
        }

        return ret;
    }

    // refraction

    error_t refractionModel(refract_model_t* model)
    {
        if (model == nullptr) {
            return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
        }

        std::lock_guard lock{*_stateMutex};

        eraRefco(_currentState.meteo.pressure, _currentState.meteo.temperature, _currentState.meteo.humidity,
                 _currentState.wavelength, &model->refa, &model->refb);

        return MccCCTE_ERFAErrorCode::ERROR_OK;
    }

    // Zobs must be observed zenithal distance (Zapp = Zobs + dZ -- corrected (in vacuo) zenithal distance)
    template <typename ZAPP_T = std::nullptr_t>
    error_t refractionCorrection(mcc_angle_c auto Zobs, mcc_angle_c auto* dZ, ZAPP_T Zapp = nullptr)
        requires(std::is_null_pointer_v<ZAPP_T> ||
                 (std::is_pointer_v<ZAPP_T> && mcc_angle_c<std::remove_pointer_t<ZAPP_T>>))
    {
        error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;

        if (dZ == nullptr) {
            return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
        }

        refract_model_t rmodel;
        ret = refractionModel(&rmodel);

        if (!ret) {
            ret = refractionCorrection(rmodel, Zobs, dZ, Zapp);
        }

        return ret;
    }

    // Zobs must be observed zenithal distance (Zapp = Zobs + dZ -- corrected (in vacuo) zenithal distance)
    template <typename ZAPP_T = std::nullptr_t>
    error_t refractionCorrection(const refract_model_t& rmodel,
                                 mcc_angle_c auto Zobs,
                                 mcc_angle_c auto* dZ,
                                 ZAPP_T Zapp = nullptr)
        requires(std::is_null_pointer_v<ZAPP_T> ||
                 (std::is_pointer_v<ZAPP_T> && mcc_angle_c<std::remove_pointer_t<ZAPP_T>>))
    {
        error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;

        if (dZ == nullptr) {
            return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
        }


        if (Zobs >= std::numbers::pi / 2.0) {
            *dZ = 35.4 / 60.0 * std::numbers::pi / 180.0;  // 35.4 arcminutes
        } else {
            auto tanZ = tan(Zobs);
            *dZ = rmodel.refa * tanZ + rmodel.refb * tanZ * tanZ * tanZ;
        }

        if constexpr (!std::is_null_pointer_v<ZAPP_T>) {
            *Zapp = Zobs + *dZ;
        }

        return ret;
    }


    // Zapp must be topocentric (in vacuo) zenithal distance (Zobs = Zapp - dZ -- observed, i.e. affected by refraction,
    // zenithal distance)
    template <typename ZOBS_T = std::nullptr_t>
    error_t refractionReverseCorrection(mcc_angle_c auto Zapp, mcc_angle_c auto* dZ, ZOBS_T Zobs = nullptr)
        requires(std::is_null_pointer_v<ZOBS_T> ||
                 (std::is_pointer_v<ZOBS_T> && mcc_angle_c<std::remove_pointer_t<ZOBS_T>>))
    {
        error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;

        if (dZ == nullptr) {
            return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
        }

        refract_model_t rmodel;
        ret = refractionModel(&rmodel);

        if (!ret) {
            ret = refractionReverseCorrection(rmodel, Zapp, dZ, Zobs);
        }

        return ret;
    }


    // Zapp must be topocentric (in vacuo) zenithal distance (Zobs = Zapp - dZ -- observed, i.e. affected by refraction,
    // zenithal distance)
    template <typename ZOBS_T = std::nullptr_t>
    error_t refractionReverseCorrection(const refract_model_t& rmodel,
                                        mcc_angle_c auto Zapp,
                                        mcc_angle_c auto* dZ,
                                        ZOBS_T Zobs = nullptr)
        requires(std::is_null_pointer_v<ZOBS_T> ||
                 (std::is_pointer_v<ZOBS_T> && mcc_angle_c<std::remove_pointer_t<ZOBS_T>>))
    {
        error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;

        if (dZ == nullptr) {
            return MccCCTE_ERFAErrorCode::ERROR_NULLPTR;
        }


        if (Zapp >= std::numbers::pi / 2.0) {
            *dZ = 35.4 / 60.0 * std::numbers::pi / 180.0;  // 35.4 arcminutes
        } else {
            auto tanZ = tan(Zapp);
            auto tanZ2 = tanZ * tanZ;
            auto b3 = 3.0 * rmodel.refb;

            // with Newton-Raphson correction
            *dZ = (rmodel.refa * tanZ + rmodel.refb * tanZ * tanZ2) /
                  (1.0 + rmodel.refa + tanZ2 * (rmodel.refa + b3) + b3 * tanZ2 * tanZ2);
        }

        if constexpr (!std::is_null_pointer_v<ZOBS_T>) {
            *Zobs = Zapp - *dZ;
        }

        return ret;
    }


    /*  helper mathods  */

    auto leapSecondsExpireDate() const
    {
        return _currentState._leapSeconds.expireDate();
    }

    auto leapSecondsExpireMJD() const
    {
        return _currentState._leapSeconds.expireMJD();
    }


    auto bulletinADateRange() const
    {
        return _currentState._bulletinA.dateRange();
    }

    auto bulletinADateRangeMJD() const
    {
        return _currentState._bulletinA.dateRangeMJD();
    }


protected:
    engine_state_t _currentState{};

    std::unique_ptr<std::mutex> _stateMutex;


    error_t icrsTo(bool observed,  // true - observed, false - apparent
                   mcc_angle_c auto const& ra_icrs,
                   mcc_angle_c auto const& dec_icrs,
                   mcc_coord_epoch_c auto const& epoch,
                   mcc_angle_c auto* ra,
                   mcc_angle_c auto* dec,
                   mcc_angle_c auto* ha,
                   mcc_angle_c auto* az,
                   mcc_angle_c auto* zd,
                   obj_pars_t* obj_params = nullptr)
    {
        int err;
        double r, d, h, a, z, eo;
        double pressure = 0.0;  // 0 for apparent coordinates type (see ERFA's refco.c: if pressure is zero then
        // refraction is also zero)


        error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;

        std::lock_guard lock{*_stateMutex};

        if (observed) {
            pressure = _currentState.meteo.pressure;
        }

        auto dut1 = _currentState._bulletinA.DUT1(epoch.MJD());

        if (!dut1.has_value()) {
            return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
        }

        auto pol_pos = _currentState._bulletinA.polarCoords(epoch.MJD());
        if (!pol_pos.has_value()) {
            return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
        }

        // const auto arcsec2rad = std::numbers::pi / 180 / 3600;
        const auto arcsec2rad = 1.0_arcsecs;
        pol_pos->x *= arcsec2rad;
        pol_pos->y *= arcsec2rad;

        if (obj_params) {
            err = eraAtco13(ra_icrs, dec_icrs, obj_params->pm_RA, obj_params->pm_DEC, obj_params->parallax,
                            obj_params->radvel, ERFA_DJM0, epoch.MJD(), dut1->count(), _currentState.lon,
                            _currentState.lat, _currentState.elev, pol_pos->x, pol_pos->y, pressure,
                            _currentState.meteo.temperature, _currentState.meteo.humidity, _currentState.wavelength, &a,
                            &z, &h, &d, &r, &eo);
        } else {
            err = eraAtco13(ra_icrs, dec_icrs, 0.0, 0.0, 0.0, 0.0, ERFA_DJM0, epoch.MJD(), dut1->count(),
                            _currentState.lon, _currentState.lat, _currentState.elev, pol_pos->x, pol_pos->y, pressure,
                            _currentState.meteo.temperature, _currentState.meteo.humidity, _currentState.wavelength, &a,
                            &z, &h, &d, &r, &eo);
        }

        if (err == 1) {
            ret = MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR;
        } else if (err == -1) {
            ret = MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE;
        }


        if (ra) {
            *ra = r;
        }

        if (dec) {
            *dec = d;
        }

        if (ha) {
            *ha = h;
        }

        if (az) {
            // NOTE: according to definition of astronomical azimuth it is counted from the South through the West, but
            // in the ERFA the azimuth is counted from the North through the East!!!
            //
            *az = MccAngle(a - std::numbers::pi).normalize<MccAngle::NORM_KIND_0_360>();
            // *az = MccAngle(a + std::numbers::pi).normalize<MccAngle::NORM_KIND_0_360>();
        }

        if (zd) {
            *zd = z;
        }

        return ret;
    }

    error_t toICRS(bool observed,  // true - observed, false - apparent
                   MccCoordPairKind pair_type,
                   mcc_coord_epoch_c auto const& epoch,
                   mcc_angle_c auto const& co_lon,
                   mcc_angle_c auto const& co_lat,
                   mcc_angle_c auto* ra_icrs,
                   mcc_angle_c auto* dec_icrs)
    {
        error_t ret = MccCCTE_ERFAErrorCode::ERROR_OK;

        // check coordinate pair consistency
        if (mcc_is_app_coordpair(pair_type) && observed) {
            return MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR;
        }
        if (mcc_is_obs_coordpair(pair_type) && !observed) {
            return MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR;
        }

        std::lock_guard lock{*_stateMutex};

        auto dut1 = _currentState._bulletinA.DUT1(epoch.MJD());

        if (!dut1.has_value()) {
            return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
        }

        auto pol_pos = _currentState._bulletinA.polarCoords(epoch.MJD());
        if (!pol_pos.has_value()) {
            return MccCCTE_ERFAErrorCode::ERROR_BULLETINA_OUT_OF_RANGE;
        }

        // const auto arcsec2rad = std::numbers::pi / 180 / 3600;
        const auto arcsec2rad = 1.0_arcsecs;
        pol_pos->x *= arcsec2rad;
        pol_pos->y *= arcsec2rad;

        std::string type;
        double x, y, ra, dec;
        double pressure = 0.0;

        if (observed) {
            pressure = _currentState.meteo.pressure;
        }

        switch (pair_type) {
            case mcc::MccCoordPairKind::COORDS_KIND_AZZD:
                // NOTE: according to definition of astronomical azimuth it is counted from the South through the West,
                // but in the ERFA the azimuth is counted from the North through the East!!!
                //
                x = co_lon + std::numbers::pi;
                y = co_lat;
                type = "A";
                break;
            case mcc::MccCoordPairKind::COORDS_KIND_AZALT:
                // NOTE: according to definition of astronomical azimuth it is counted from the South through the West,
                // but in the ERFA the azimuth is counted from the North through the East!!!
                //
                x = co_lon + std::numbers::pi;
                y = MccCCTE_ERFA::PI_2 - co_lat;  // altitude to zenithal distance
                type = "A";
                break;
            case mcc::MccCoordPairKind::COORDS_KIND_HADEC_OBS:
                type = "H";
                x = co_lon;
                y = co_lat;
                break;
            case mcc::MccCoordPairKind::COORDS_KIND_RADEC_OBS:
                type = "R";
                x = co_lon;
                y = co_lat;
                break;
            case mcc::MccCoordPairKind::COORDS_KIND_HADEC_APP:
                type = "H";
                x = co_lon;
                y = co_lat;
                break;
            case mcc::MccCoordPairKind::COORDS_KIND_RADEC_APP:
                type = "R";
                x = co_lon;
                y = co_lat;
                break;
            default:
                return MccCCTE_ERFAErrorCode::ERROR_UNSUPPORTED_COORD_PAIR;
        };

        int err =
            eraAtoc13(type.c_str(), x, y, ERFA_DJM0, epoch.MJD(), dut1->count(), _currentState.lon, _currentState.lat,
                      _currentState.elev, pol_pos->x, pol_pos->y, pressure, _currentState.meteo.temperature,
                      _currentState.meteo.humidity, _currentState.wavelength, &ra, &dec);

        if (err == 1) {
            ret = MccCCTE_ERFAErrorCode::ERROR_DUBIOUS_YEAR;
        } else if (err == -1) {
            ret = MccCCTE_ERFAErrorCode::ERROR_UNACCEPTABLE_DATE;
        }

        if (ra) {
            *ra_icrs = ra;
        }

        if (dec) {
            *dec_icrs = dec;
        }

        return ret;
    }
};


}  // namespace mcc::ccte::erfa