#pragma once

/****************************************************************************************
 *                                                                                      *
 *                          MOUNT CONTROL COMPONENTS LIBRARY                            *
 *                                                                                      *
 *                                                                                      *
 *                      IMPLEMENTATION OF CELESTIAL COORDINATES                         *
 *                                                                                      *
 ****************************************************************************************/



// #include "mcc_angle.h"
#include "mcc_ccte_erfa.h"
#include "mcc_concepts.h"
#include "mcc_epoch.h"

namespace mcc::impl
{


/*                  CLASSES TO REPRESENT COORDINATES PAIR                   */

template <mcc_angle_c CO_LON_T, mcc_angle_c CO_LAT_T>
class MccCoordPair : public mcc_coord_pair_interface_t
{
public:
    typedef CO_LON_T x_t;
    typedef CO_LAT_T y_t;

    static constexpr MccCoordPairKind pairKind =
        !(std::derived_from<CO_LON_T, MccAngle> ||
          std::derived_from<CO_LAT_T, MccAngle>)  // unknown type (possibly just double or float)
            ? MccCoordPairKind::COORDS_KIND_GENERIC
            : (std::same_as<CO_LON_T, MccAngle> || std::same_as<CO_LAT_T, MccAngle>)  // one of the types is MccAngle
                  ? MccCoordPairKind::COORDS_KIND_GENERIC
                  // ICRS RA and DEC
                  : (std::same_as<CO_LON_T, MccAngleRA_ICRS> && std::same_as<CO_LAT_T, MccAngleDEC_ICRS>)
                        ? MccCoordPairKind::COORDS_KIND_RADEC_ICRS
                    // apparent RA and DEC
                    : (std::same_as<CO_LON_T, MccAngleRA_APP> && std::same_as<CO_LAT_T, MccAngleDEC_APP>)
                        ? MccCoordPairKind::COORDS_KIND_RADEC_APP
                    // observed RA and DEC
                    : (std::same_as<CO_LON_T, MccAngleRA_OBS> && std::same_as<CO_LAT_T, MccAngleDEC_OBS>)
                        ? MccCoordPairKind::COORDS_KIND_RADEC_OBS
                    // apparent HA and DEC
                    : (std::same_as<CO_LON_T, MccAngleHA_APP> && std::same_as<CO_LAT_T, MccAngleDEC_APP>)
                        ? MccCoordPairKind::COORDS_KIND_HADEC_APP
                    // observed HA and DEC
                    : (std::same_as<CO_LON_T, MccAngleHA_OBS> && std::same_as<CO_LAT_T, MccAngleDEC_OBS>)
                        ? MccCoordPairKind::COORDS_KIND_HADEC_OBS
                    // apparent AZ and ZD
                    : (std::same_as<CO_LON_T, MccAngleAZ> && std::same_as<CO_LAT_T, MccAngleZD>)
                        ? MccCoordPairKind::COORDS_KIND_AZZD
                    // apparent AZ and ALT
                    : (std::same_as<CO_LON_T, MccAngleAZ> && std::same_as<CO_LAT_T, MccAngleALT>)
                        ? MccCoordPairKind::COORDS_KIND_AZALT
                    // general purpose X and Y
                    : (std::same_as<CO_LON_T, MccAngleX> && std::same_as<CO_LAT_T, MccAngleY>)
                        ? MccCoordPairKind::COORDS_KIND_XY
                    // geographical longitude and latitude
                    : (std::same_as<CO_LON_T, MccAngleLON> && std::same_as<CO_LAT_T, MccAngleLAT>)
                        ? MccCoordPairKind::COORDS_KIND_LONLAT
                        : MccCoordPairKind::COORDS_KIND_UNKNOWN;

    MccCoordPair() : _x(0.0), _y(0.0), _epoch(MccCelestialCoordEpoch::now()) {}

    template <mcc_coord_epoch_c EpT = MccCelestialCoordEpoch>
    MccCoordPair(CO_LON_T const& x, CO_LAT_T const& y, EpT const& epoch = EpT::now()) : _x(x), _y(y), _epoch(epoch)
    {
        normalize();
    }

    MccCoordPair(const MccCoordPair&) = default;
    MccCoordPair(MccCoordPair&&) = default;

    MccCoordPair& operator=(const MccCoordPair&) = default;
    MccCoordPair& operator=(MccCoordPair&&) = default;

    template <mcc_coord_pair_c T>
        requires(T::pairKind == pairKind || T::pairKind == MccCoordPairKind::COORDS_KIND_GENERIC ||
                 T::pairKind == MccCoordPairKind::COORDS_KIND_XY)
    MccCoordPair(const T& other)
    {
        setX((double)other.x());
        setY((double)other.y());
        setEpoch(other.epoch());

        normalize();
    }


    template <mcc_coord_pair_c T>
        requires(T::pairKind == pairKind || T::pairKind == MccCoordPairKind::COORDS_KIND_GENERIC ||
                 T::pairKind == MccCoordPairKind::COORDS_KIND_XY)
    MccCoordPair(T&& other)
    {
        setX((double)other.x());
        setY((double)other.y());
        setEpoch(other.epoch());

        normalize();
    }


    template <mcc_coord_pair_c T>
        requires(T::pairKind == pairKind || T::pairKind == MccCoordPairKind::COORDS_KIND_GENERIC ||
                 T::pairKind == MccCoordPairKind::COORDS_KIND_XY)
    MccCoordPair& operator=(const T& other)
    {
        setX((double)other.x());
        setY((double)other.y());
        setEpoch(other.epoch());

        normalize();
    }


    template <mcc_coord_pair_c T>
        requires(T::pairKind == pairKind || T::pairKind == MccCoordPairKind::COORDS_KIND_GENERIC ||
                 T::pairKind == MccCoordPairKind::COORDS_KIND_XY)
    MccCoordPair& operator=(T&& other)
    {
        setX((double)other.x());
        setY((double)other.y());
        setEpoch(other.epoch());

        normalize();
    }


    virtual ~MccCoordPair() = default;

    CO_LON_T x() const
    {
        return _x;
    }

    CO_LAT_T y() const
    {
        return _y;
    }


    MccCelestialCoordEpoch epoch() const
    {
        return _epoch;
    }

    template <mcc_coord_epoch_c EpT>
    EpT epoch() const
    {
        return _epoch;
    }


    double MJD() const
    {
        return _epoch.MJD();
    }

    // for something like:
    //    auto [ra, dec, epoch] = coord_pair;
    operator std::tuple<CO_LON_T, CO_LAT_T, MccCelestialCoordEpoch>() const
    {
        return {_x, _y, _epoch};
    }

    void setX(const CO_LON_T& x)
    {
        _x = x;

        normalize();
    }

    void setY(const CO_LAT_T& y)
    {
        _y = y;

        normalize();
    }

    void setEpoch(mcc_coord_epoch_c auto const& ep)
    {
        _epoch = ep;
    }

protected:
    CO_LON_T _x;
    CO_LAT_T _y;

    MccCelestialCoordEpoch _epoch;

    void normalize()
    {
        if constexpr (pairKind != MccCoordPairKind::COORDS_KIND_GENERIC &&
                      pairKind != MccCoordPairKind::COORDS_KIND_XY) {
            if constexpr (pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP ||
                          pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
                _x = (double)MccAngle(_x).normalize<MccAngle::NORM_KIND_180_180>();
            } else {  // RA, AZ
                _x = (double)MccAngle(_x).normalize<MccAngle::NORM_KIND_0_360>();
            }

            // DEC and ALT is [-90,90] degrees
            if constexpr (pairKind != MccCoordPairKind::COORDS_KIND_AZZD) {
                _y = (double)MccAngle(_y).normalize<MccAngle::NORM_KIND_90_90>();
            } else {  // ZD id [0, 180] degrees
                _y = (double)MccAngle(_y).normalize<MccAngle::NORM_KIND_0_180>();
            }
        }
    }
};


static_assert(mcc_coord_pair_c<MccCoordPair<MccAngleRA_ICRS, MccAngleDEC_ICRS>>, "");


/*                          PREDEFINED COORDINATES PAIR TYPES                            */



struct MccSkyRADEC_ICRS : MccCoordPair<MccAngleRA_ICRS, MccAngleDEC_ICRS> {
    // re-implement constructors to keep the epoch equal to J2000.0

    MccSkyRADEC_ICRS() : MccCoordPair<MccAngleRA_ICRS, MccAngleDEC_ICRS>(0.0, 0.0, MccCelestialCoordEpoch{}) {}

    MccSkyRADEC_ICRS(MccAngleRA_ICRS const& x, MccAngleDEC_ICRS const& y)
        : MccCoordPair<MccAngleRA_ICRS, MccAngleDEC_ICRS>((double)x, (double)y, MccCelestialCoordEpoch{})
    {
    }

    // ignore epoch setting (it is always J2000.0)
    void setEpoch(mcc_coord_epoch_c auto const&)
    {
        // static_assert(false, "CANNOT SET EPOCH FOR ICRS-KIND COORDINATE PAIR!!!");
    }
};

using MccSkyRADEC_APP = MccCoordPair<MccAngleRA_APP, MccAngleDEC_APP>;

using MccSkyRADEC_OBS = MccCoordPair<MccAngleRA_OBS, MccAngleDEC_OBS>;

using MccSkyHADEC_APP = MccCoordPair<MccAngleHA_APP, MccAngleDEC_APP>;

using MccSkyHADEC_OBS = MccCoordPair<MccAngleHA_OBS, MccAngleDEC_OBS>;

// using MccSkyAZZD = MccCoordPair<MccAngleAZ, MccAngleZD>;
struct MccSkyAZZD : MccCoordPair<MccAngleAZ, MccAngleZD> {
    using MccCoordPair<MccAngleAZ, MccAngleZD>::MccCoordPair;

    template <mcc_coord_pair_c AZALT_PAIR_T>
        requires(AZALT_PAIR_T::pairKind == MccCoordPairKind::COORDS_KIND_AZALT)
    MccSkyAZZD(AZALT_PAIR_T const& azalt) : MccSkyAZZD()
    {
        setX((double)azalt.x());
        // setY(std::numbers::pi / 2.0 - (double)azalt.x());
        setY(MCC_HALF_PI - (double)azalt.x());
        setEpoch(azalt.epoch());
    }

    template <mcc_coord_pair_c AZALT_PAIR_T>
        requires(AZALT_PAIR_T::pairKind == MccCoordPairKind::COORDS_KIND_AZALT)
    MccSkyAZZD(AZALT_PAIR_T&& azalt) : MccSkyAZZD()
    {
        setX((double)azalt.x());
        // setY(std::numbers::pi / 2.0 - (double)azalt.x());
        setY(MCC_HALF_PI - (double)azalt.x());
        setEpoch(azalt.epoch());
    }

    template <mcc_coord_pair_c AZALT_PAIR_T>
        requires(AZALT_PAIR_T::pairKind == MccCoordPairKind::COORDS_KIND_AZALT)
    MccSkyAZZD& operator=(AZALT_PAIR_T const& azalt)
    {
        setX((double)azalt.x());
        // setY(std::numbers::pi / 2.0 - (double)azalt.x());
        setY(MCC_HALF_PI - (double)azalt.x());
        setEpoch(azalt.epoch());

        return *this;
    }

    template <mcc_coord_pair_c AZALT_PAIR_T>
        requires(AZALT_PAIR_T::pairKind == MccCoordPairKind::COORDS_KIND_AZALT)
    MccSkyAZZD& operator=(AZALT_PAIR_T&& azalt)
    {
        setX((double)azalt.x());
        // setY(std::numbers::pi / 2.0 - (double)azalt.x());
        setY(MCC_HALF_PI - (double)azalt.x());
        setEpoch(azalt.epoch());

        return *this;
    }
};

// using MccSkyAZALT = MccCoordPair<MccAngleAZ, MccAngleALT>;
struct MccSkyAZALT : MccCoordPair<MccAngleAZ, MccAngleALT> {
    using MccCoordPair<MccAngleAZ, MccAngleALT>::MccCoordPair;

    template <mcc_coord_pair_c AZZD_PAIR_T>
        requires(AZZD_PAIR_T::pairKind == MccCoordPairKind::COORDS_KIND_AZZD)
    MccSkyAZALT(AZZD_PAIR_T const& azzd) : MccSkyAZALT()
    {
        setX((double)azzd.x());
        // setY(std::numbers::pi / 2.0 - (double)azzd.x());
        setY(MCC_HALF_PI - (double)azzd.x());
        setEpoch(azzd.epoch());
    }

    template <mcc_coord_pair_c AZZD_PAIR_T>
        requires(AZZD_PAIR_T::pairKind == MccCoordPairKind::COORDS_KIND_AZZD)
    MccSkyAZALT(AZZD_PAIR_T&& azzd) : MccSkyAZALT()
    {
        setX((double)azzd.x());
        // setY(std::numbers::pi / 2.0 - (double)azzd.x());
        setY(MCC_HALF_PI - (double)azzd.x());
        setEpoch(azzd.epoch());
    }

    template <mcc_coord_pair_c AZZD_PAIR_T>
        requires(AZZD_PAIR_T::pairKind == MccCoordPairKind::COORDS_KIND_AZZD)
    MccSkyAZALT& operator=(AZZD_PAIR_T const& azzd)
    {
        setX((double)azzd.x());
        // setY(std::numbers::pi / 2.0 - (double)azzd.x());
        setY(MCC_HALF_PI - (double)azzd.x());
        setEpoch(azzd.epoch());

        return *this;
    }

    template <mcc_coord_pair_c AZZD_PAIR_T>
        requires(AZZD_PAIR_T::pairKind == MccCoordPairKind::COORDS_KIND_AZZD)
    MccSkyAZALT& operator=(AZZD_PAIR_T&& azzd)
    {
        setX((double)azzd.x());
        // setY(std::numbers::pi / 2.0 - (double)azzd.x());
        setY(MCC_HALF_PI - (double)azzd.x());
        setEpoch(azzd.epoch());

        return *this;
    }
};

using MccGenXY = MccCoordPair<MccAngleX, MccAngleY>;

using MccGeoLONLAT = MccCoordPair<MccAngleLON, MccAngleLAT>;

// utility type definition: deduce a coordinate pair kind according to mount type
// mcc_deduced_coord_pair_t = MccSkyHADEC_OBS for equathorial mounts,
// mcc_deduced_coord_pair_t = MccSkyAZZD for altazimuthal ones and
// mcc_deduced_coord_pair_t = std::nullptr_t otherwise
template <MccMountType MOUNT_TYPE>
using mcc_deduced_coord_pair_t =
    std::conditional_t<mccIsEquatorialMount(MOUNT_TYPE),
                       MccSkyHADEC_OBS,
                       std::conditional_t<mccIsAltAzMount(MOUNT_TYPE), MccSkyAZZD, std::nullptr_t>>;



// static MccSkyHADEC_APP hadec = MccGenXY{};
// static MccSkyAZALT azalt{MccSkyAZZD{1.0, 1.1}};


/*                          MCC-LIBRARY DEFAULT GENERIC SKY POINT CLASS IMPLEMENTATION                    */



template <mcc_ccte_c CCTE_T>
class MccGenericSkyPoint : public mcc_skypoint_interface_t
{
public:
    typedef CCTE_T ccte_t;

    static constexpr double MJD0 = 2400000.5;

    inline static CCTE_T cctEngine{};  // celestial coordinates transformation engine

    using error_t = typename CCTE_T::error_t;

    MccGenericSkyPoint() {}

    template <mcc_coord_pair_c PT>
    MccGenericSkyPoint(const PT& coord_pair) : MccGenericSkyPoint()
    {
        auto self = from(coord_pair);
    }

    MccGenericSkyPoint(const MccGenericSkyPoint&) = default;
    MccGenericSkyPoint(MccGenericSkyPoint&&) = default;

    MccGenericSkyPoint& operator=(const MccGenericSkyPoint&) = default;
    MccGenericSkyPoint& operator=(MccGenericSkyPoint&&) = default;


    MccGenericSkyPoint(mcc_skypoint_c auto const& other)
    {
        fromOtherSkyPoint(other);
    }

    MccGenericSkyPoint(mcc_skypoint_c auto&& other)
    {
        fromOtherSkyPoint(other);
    }


    MccGenericSkyPoint& operator=(mcc_skypoint_c auto const& other)
    {
        fromOtherSkyPoint(other);

        return *this;
    }

    MccGenericSkyPoint& operator=(mcc_skypoint_c auto&& other)
    {
        fromOtherSkyPoint(other);

        return *this;
    }


    virtual ~MccGenericSkyPoint() = default;

    MccCoordPairKind pairKind() const
    {
        return _pairKind;
    }

    MccCelestialCoordEpoch epoch() const
    {
        return _epoch;
    }


    double co_lon() const
    {
        return _x;
    }

    double co_lat() const
    {
        return _y;
    }

    template <mcc_coord_pair_c PT>
    MccGenericSkyPoint& from(const PT& coord_pair)
    {
        _x = coord_pair.x();
        _y = coord_pair.y();

        _pairKind = PT::pairKind;

        if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
            _epoch = MccCelestialCoordEpoch();  // J2000.0
        } else {
            _epoch.fromMJD(coord_pair.MJD());
        }

        return *this;
    }

    MccGenericSkyPoint& operator=(mcc_coord_pair_c auto const& coord_pair)
    {
        return from(coord_pair);
    }


    template <mcc_coord_pair_c PT, mcc_coord_pair_c... PTs>
    error_t to(PT& cpair, PTs&... cpairs) const
    {
        auto err = toHelper(cpair);
        if (err) {
            return err;
        }

        if constexpr (sizeof...(PTs)) {
            err = to(cpairs...);
            if (err) {
                return err;
            }
        }

        // according to mcc_error_c concept (see mcc_concepts.h)
        // default-constructed mcc_error_c-like class must be assumed as
        // non-error state
        return error_t{};
    }


    template <mcc_coord_pair_c PT>
    operator PT() const
    {
        if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_LONLAT) {  // returns geographic site coordinates
            std::pair<double, double> pos;
            cctEngine.geoPosition(&pos);

            return MccGeoLONLAT(pos.first, pos.second);
        }


        PT res;
        toAtSameEpoch(res);

        return res;
    }


    template <mcc_coord_pair_c PT, mcc_coord_pair_c... PTs>
    error_t toAtSameEpoch(PT& cpair, PTs&... cpairs) const
    {
        if constexpr (PT::pairKind != MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
            if (_pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {  // from ICRS: set current epoch for result
                cpair.setEpoch(MccCelestialCoordEpoch::now());
            } else {
                cpair.setEpoch(_epoch);
            }
        }


        auto err = toHelper(cpair);
        if (err) {
            return err;
        }

        if constexpr (sizeof...(PTs)) {
            err = toAtSameEpoch(cpairs...);
            if (err) {
                return err;
            }
        }

        // according to mcc_error_c concept (see mcc_concepts.h)
        // default-constructed mcc_error_c-like class must be assumed as
        // non-error state
        return error_t{};
    }


    error_t refractCorrection(mcc_angle_c auto* dZ) const
    {
        if (mcc_is_obs_coordpair(_pairKind)) {
            if (_pairKind == MccCoordPairKind::COORDS_KIND_AZZD) {
                return cctEngine.refractionCorrection(_y, dZ);
            } else if (_pairKind == MccCoordPairKind::COORDS_KIND_AZALT) {
                return cctEngine.refractionCorrection(MCC_HALF_PI - _y, dZ);
            } else {
                MccSkyAZZD azzd;
                auto err = toAtSameEpoch(azzd);
                if (!err) {
                    err = cctEngine.refractionCorrection(azzd.y(), dZ);
                }

                return err;
            }
        } else {
            if (dZ) {
                *dZ = 0.0;
            }

            return {};
        }
    }


    error_t refractInverseCorrection(mcc_angle_c auto* dZ) const
    {
        double phi = cctEngine.getStateERFA().lat;
        double ha = _x, dec = _y;

        if (mcc_is_app_coordpair(_pairKind)) {
            double az, alt;

            if (_pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP) {
                double eo, lst;

                auto ccte_err = cctEngine.equationOrigins(_epoch, &eo);
                if (ccte_err) {
                    return ccte_err;
                }

                ccte_err = cctEngine.apparentSideralTime(_epoch, &lst, true);
                if (ccte_err) {
                    return ccte_err;
                }

                // from RA to HA
                ha = lst + eo - _x;
            }

            hadec2azalt(ha, dec, phi, &az, &alt);

            return cctEngine.refractionInverseCorrection(MCC_HALF_PI - alt, dZ);
        } else {
            if (dZ) {
                *dZ = 0.0;
            }

            return {};
        }
    }


    error_t appSideralTime(mcc_angle_c auto* st, bool is_local = false) const
    {
        return cctEngine.apparentSideralTime(_epoch, st, is_local);
    }


    error_t EO(mcc_angle_c auto* eo) const
    {
        return cctEngine.equationOrigins(_epoch, eo);
    }

protected:
    double _x{0.0}, _y{0.0};
    MccCoordPairKind _pairKind{MccCoordPairKind::COORDS_KIND_RADEC_ICRS};
    MccCelestialCoordEpoch _epoch{};  // J2000.0

    template <mcc_skypoint_c T>
    void fromOtherSkyPoint(T&& other)
    {
        switch (other.pairKind()) {
            case MccCoordPairKind::COORDS_KIND_RADEC_ICRS: {
                MccSkyRADEC_ICRS pt;

                other.to(pt);
                from(pt);
            } break;
            case MccCoordPairKind::COORDS_KIND_RADEC_OBS: {
                MccSkyRADEC_OBS pt;
                pt.setEpoch(other.epoch());

                other.to(pt);
                from(pt);
            } break;
            case MccCoordPairKind::COORDS_KIND_RADEC_APP: {
                MccSkyRADEC_APP pt;
                pt.setEpoch(other.epoch());

                other.to(pt);
                from(pt);
            } break;
            case MccCoordPairKind::COORDS_KIND_HADEC_OBS: {
                MccSkyHADEC_OBS pt;
                pt.setEpoch(other.epoch());

                other.to(pt);
                from(pt);
            } break;
            case MccCoordPairKind::COORDS_KIND_HADEC_APP: {
                MccSkyHADEC_APP pt;
                pt.setEpoch(other.epoch());

                other.to(pt);
                from(pt);
            } break;
            case MccCoordPairKind::COORDS_KIND_AZALT: {
                MccSkyAZALT pt;
                pt.setEpoch(other.epoch());

                other.to(pt);
                from(pt);
            } break;
            case MccCoordPairKind::COORDS_KIND_AZZD: {
                MccSkyAZZD pt;
                pt.setEpoch(other.epoch());

                other.to(pt);
                from(pt);
            } break;
            default:
                // error!!!
                break;
        }
    }

    // HA, DEC to AZ, ALT (AZ from the South through the West)
    void hadec2azalt(double ha, double dec, double phi, double& az, double& alt) const
    {
        // eraHd2ae(ha, dec, phi, &az, &alt);
        // // from ERFA "from N" to "from S"
        // if (az > std::numbers::pi) {
        //     az -= std::numbers::pi;
        // } else {
        //     az += std::numbers::pi;
        // }

        // return;

        const auto cos_phi = std::cos(phi), sin_phi = std::sin(phi);
        const auto cos_dec = std::cos(dec), sin_dec = std::sin(dec);
        const auto cos_ha = std::cos(ha), sin_ha = std::sin(ha);

        auto x = sin_phi * cos_dec * cos_ha - cos_phi * sin_dec;
        auto y = cos_dec * sin_ha;
        auto z = cos_phi * cos_dec * cos_ha + sin_phi * sin_dec;

        auto xx = x * x, yy = y * y;
        decltype(x) r;
        if (xx < yy) {
            r = std::abs(y) * sqrt(1.0 + xx / yy);
        } else {
            r = std::abs(x) * sqrt(1.0 + yy / xx);
        }

        az = utils::isEqual(r, 0.0) ? 0.0 : std::atan2(y, x);
        if (az < 0.0) {
            az += MCC_TWO_PI;  // to range of [0, 2*PI]
        }

        alt = std::atan2(z, r);
    };


    // AZ, ALT to HA, DEC (AZ from the South through the West)
    void azalt2hadec(double az, double alt, double phi, double& ha, double& dec) const
    {
        // az += std::numbers::pi;
        // eraAe2hd(az, alt, phi, &ha, &dec);

        // return;

        const auto cos_phi = std::cos(phi), sin_phi = std::sin(phi);
        const auto cos_az = std::cos(az), sin_az = std::sin(az);
        const auto cos_alt = std::cos(alt), sin_alt = std::sin(alt);

        auto x = sin_phi * cos_alt * cos_az + cos_phi * sin_alt;
        auto y = cos_alt * sin_az;
        auto z = -cos_phi * cos_alt * cos_az + sin_phi * sin_alt;

        auto xx = x * x, yy = y * y;
        decltype(x) r;
        if (xx < yy) {
            r = std::abs(y) * sqrt(1.0 + xx / yy);
        } else {
            r = std::abs(x) * sqrt(1.0 + yy / xx);
        }

        ha = utils::isEqual(r, 0.0) ? 0.0 : std::atan2(y, x);
        dec = std::atan2(z, r);
    };


    template <mcc_coord_pair_c PT>
    error_t toHelper(PT& cpair) const
    {
        typename CCTE_T::error_t ccte_err;

        double phi = cctEngine.getStateERFA().lat;
        double ra_icrs, dec_icrs, ra, dec, ha, az, zd, alt, lst, eo;

        static_assert(PT::pairKind != MccCoordPairKind::COORDS_KIND_GENERIC, "UNSUPPORTED SKY POINT TRANSFORMATION!");
        static_assert(PT::pairKind != MccCoordPairKind::COORDS_KIND_UNKNOWN, "UNSUPPORTED SKY POINT TRANSFORMATION!");


        if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_LONLAT) {  // returns geographic site coordinates
            std::pair<double, double> pos;
            cctEngine.geoPosition(&pos);

            cpair.setX(pos.first);
            cpair.setY(pos.second);

            return error_t{};
        }

        if (_pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS &&
            PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {  // from ICRS to ICRS - just copy and exit
            cpair = PT(typename PT::x_t(_x), typename PT::y_t(_y));

            return error_t{};
        }

        // just copy coordinates and exit
        if (_pairKind == PT::pairKind && utils::isEqual(_epoch.MJD(), cpair.MJD())) {
            // cpair = PT(typename PT::x_t(_x), typename PT::y_t(_y), _epoch);
            cpair.setX(_x);
            cpair.setY(_y);

            return error_t{};
        }


        // if epochs are not the same then
        // 1) convert stored coordinates to ICRS ones
        // 2) convert from the computed ICRS coordinates to required ones
        MccCoordPairKind pkind = _pairKind;
        if (!utils::isEqual(_epoch.MJD(), cpair.MJD())) {  // convert stored pair to ICRS one (ra_icrs, dec_icrs)
            if (_pairKind != MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
                pkind = MccCoordPairKind::COORDS_KIND_RADEC_ICRS;

                if (mcc_is_obs_coordpair(_pairKind)) {
                    ccte_err = cctEngine.obsToICRS(_pairKind, _epoch, _x, _y, &ra_icrs, &dec_icrs);
                } else if (mcc_is_app_coordpair(_pairKind)) {
                    ccte_err = cctEngine.appToICRS(_pairKind, _epoch, _x, _y, &ra_icrs, &dec_icrs);
                } else {  // unsupported transformation!!! silently ignore!
                    return error_t{};
                }

                if (ccte_err) {
                    return ccte_err;
                }
            } else {
                ra_icrs = _x;
                dec_icrs = _y;
            }
        }

        if (utils::isEqual(ra_icrs, MCC_TWO_PI)) {
            ra_icrs = 0.0;
        }


        // here, from APP or OBS to ICRS and exit
        if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS &&
            PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
            cpair = PT(typename PT::x_t(ra_icrs), typename PT::y_t(dec_icrs));

            return error_t{};
        }

        // here, the input coordinates and stored one are at the same epoch

        ccte_err = cctEngine.equationOrigins(cpair.epoch(), &eo);
        if (ccte_err) {
            return ccte_err;
        }

        ccte_err = cctEngine.apparentSideralTime(cpair.epoch(), &lst, true);
        if (ccte_err) {
            return ccte_err;
        }


        if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_APP || pkind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) {
            ra = _x;
            dec = _y;
        } else if (pkind == MccCoordPairKind::COORDS_KIND_HADEC_APP ||
                   pkind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
            ha = _x;
            dec = _y;
        } else if (pkind == MccCoordPairKind::COORDS_KIND_AZZD) {
            az = _x;
            zd = _y;
        } else if (pkind == MccCoordPairKind::COORDS_KIND_AZALT) {
            az = _x;
            alt = _y;
        }
        // else if (pkind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
        //     ra_icrs = _x;
        //     dec_icrs = _y;
        // } else {  // unsupported transformation!!!
        //     return;
        // }

        // coordinate transformation lambda (possibly recursive!!!)
        // "obj = this" to fix GCC compilation crash!!!
        auto comp_func = [&, obj = this](this auto&& self, MccCoordPairKind cp_kind) -> error_t {
            if (cp_kind == MccCoordPairKind::COORDS_KIND_RADEC_ICRS) {
                if constexpr (mccIsAppCoordPairKind<PT::pairKind>) {
                    ccte_err = cctEngine.icrsToApp(ra_icrs, dec_icrs, cpair.epoch(), &ra, &dec, &ha, &az, &zd);
                } else if constexpr (mccIsObsCoordPairKind<PT::pairKind>) {
                    ccte_err = cctEngine.icrsToObs(ra_icrs, dec_icrs, cpair.epoch(), &ra, &dec, &ha, &az, &zd);
                } else {
                    static_assert(true, "UNSUPPORTED SKY POINT TRANSFORMATION!");
                }

                if (ccte_err) {
                    return ccte_err;
                }

                if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP ||
                              PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) {
                    cpair.setX(ra);
                    cpair.setY(dec);
                } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP ||
                                     PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
                    cpair.setX(ha);
                    cpair.setY(dec);
                } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) {
                    cpair.setX(az);
                    cpair.setY(zd);
                } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) {
                    cpair.setX(az);
                    cpair.setY(MCC_HALF_PI - zd);
                } else {
                    static_assert(true, "UNSUPPORTED SKY POINT TRANSFORMATION!");
                }
            } else if (cp_kind == MccCoordPairKind::COORDS_KIND_AZALT) {
                if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) {
                    zd = MCC_HALF_PI - alt;
                    cpair.setX(az);
                    cpair.setY(zd);
                } else {
                    if constexpr (mccIsAppCoordPairKind<PT::pairKind>) {
                        // correct for refraction: alt -= dz_refr
                        double dZ;
                        ccte_err = cctEngine.refractionCorrection(MCC_HALF_PI - alt, &dZ);
                        if (ccte_err) {
                            return ccte_err;
                        }

                        alt -= dZ;
                    }

                    obj->azalt2hadec(az, alt, phi, ha, dec);
                    cpair.setY(dec);

                    if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP) {
                        ra = lst + eo - ha;
                        cpair.setX(ra);
                    } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) {
                        ra = lst + eo - ha;
                        cpair.setX(ra);
                    } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) {
                        cpair.setX(ha);
                    } else if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
                        cpair.setX(ha);
                    } else {  // unsupported transformation!!! silently ignore!!!
                        return error_t{};
                    }
                }
            } else if (cp_kind == MccCoordPairKind::COORDS_KIND_AZZD) {
                alt = MCC_HALF_PI - zd;
                if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZALT) {
                    cpair.setX(az);
                    cpair.setY(alt);
                } else {
                    return self(MccCoordPairKind::COORDS_KIND_AZALT);
                }
            } else if (cp_kind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
                if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) {
                    ra = lst + eo - ha;
                    cpair.setX(ra);
                    cpair.setY(dec);
                } else {
                    obj->hadec2azalt(ha, dec, phi, az, alt);
                    if constexpr (mccIsAppCoordPairKind<PT::pairKind>) {  // RADEC_APP, HADEC_APP
                        return self(MccCoordPairKind::COORDS_KIND_AZALT);
                    } else {  // AZALT, AZZD
                        cpair.setX(az);
                        if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_AZZD) {
                            zd = MCC_HALF_PI - alt;
                            cpair.setY(zd);
                        } else {
                            cpair.setY(alt);
                        }
                    }
                }
            } else if (cp_kind == MccCoordPairKind::COORDS_KIND_HADEC_APP) {
                if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_RADEC_APP) {
                    ra = lst + eo - ha;
                    cpair.setX(ra);
                    cpair.setY(dec);
                } else {
                    obj->hadec2azalt(ha, dec, phi, az, alt);
                    if constexpr (mccIsObsCoordPairKind<PT::pairKind>) {  // RADEC_OBS, HADEC_OBS, AZALT, AZZD
                        // correct for refraction: alt += dz_refr
                        double dZ;
                        ccte_err = cctEngine.refractionInverseCorrection(MCC_HALF_PI - alt, &dZ);
                        alt += dZ;
                        return self(MccCoordPairKind::COORDS_KIND_AZALT);
                    }
                }
            } else if (cp_kind == MccCoordPairKind::COORDS_KIND_RADEC_OBS) {
                ha = lst + eo - ra;
                if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_OBS) {
                    cpair.setX(ha);
                    cpair.setY(dec);
                } else {
                    return self(MccCoordPairKind::COORDS_KIND_HADEC_OBS);
                }
            } else if (cp_kind == MccCoordPairKind::COORDS_KIND_RADEC_APP) {
                ha = lst + eo - ra;
                if constexpr (PT::pairKind == MccCoordPairKind::COORDS_KIND_HADEC_APP) {
                    cpair.setX(ha);
                    cpair.setY(dec);
                } else {
                    return self(MccCoordPairKind::COORDS_KIND_HADEC_APP);
                }
            }

            return error_t{};
        };

        auto err = comp_func(pkind);  // ran transformation

        if (ra < 0.0) {
            ra += MCC_TWO_PI;
            cpair.setX(ra);
        }
        if (az < 0.0) {
            az += MCC_TWO_PI;
            cpair.setX(az);
        }

        // if (utils::isEqual(ra, MCC_TWO_PI)) {
        //     cpair.setX(0.0);
        // }

        // if (utils::isEqual(az, MCC_TWO_PI)) {
        //     cpair.setX(0.0);
        // }

        return err;
    }
};


/*                          MCC-LIBRARY DEFAULT SKY POINT CLASS IMPLEMENTATION BASED ON THE ERFA LIBRARY */

typedef MccGenericSkyPoint<mcc::ccte::erfa::MccCCTE_ERFA> MccSkyPoint;

static_assert(mcc_skypoint_c<MccSkyPoint>, "!!!!");

}  // namespace mcc::impl