#pragma once
#include <Eigen/Dense>

namespace teqp {

    template <typename MoleFractions>
    auto poly_grad_2(const MoleFractions& x, const double& c1, const double& c2, const double& c3) {
        return c1 + c2 * x + c3 * x * x;
    }

    // reducing function depending on composition, temperature and density
    class Reducing_ECS {

    private:

        Eigen::MatrixXd  tr_coeffs;
        Eigen::MatrixXd dr_coeffs;

    public:

        Eigen::ArrayXd Tc, vc;
        template<typename ArrayLike>
        Reducing_ECS(const ArrayLike& Tc, const ArrayLike& vc, const nlohmann::json& jj) : Tc(Tc), vc(vc) {

            auto json_tr_coeffs = jj.at("tr_coeffs");
            auto json_dr_coeffs = jj.at("dr_coeffs");

            int rows_tr = json_tr_coeffs.size();
            int cols_tr = json_tr_coeffs[0].size();
            tr_coeffs.resize(rows_tr, cols_tr);

            int rows_dr = json_dr_coeffs.size();
            int cols_dr = json_dr_coeffs[0].size();
            dr_coeffs.resize(rows_dr, cols_dr);


            for (int i = 0; i < rows_tr; ++i) {
                for (int j = 0; j < cols_tr; ++j) {
                    tr_coeffs(i, j) = json_tr_coeffs[i][j];
                }
            }

            for (int i = 0; i < rows_dr; ++i) {
                for (int j = 0; j < cols_dr; ++j) {
                    dr_coeffs(i, j) = json_dr_coeffs[i][j];
                }
            }

        }

        template <typename TTYPE, typename RHOTYPE, typename MoleFractions>
        auto get_tr(const TTYPE& T, const RHOTYPE& rho, const MoleFractions& z) const {
            //Eigen::ArrayX<decltype(z[0])> p_ij(6);
            auto p00 = tr_coeffs(0, 0) + z[0] * tr_coeffs(0, 1) + z[0] * z[0] * tr_coeffs(0, 2);
            auto p10 = tr_coeffs(1, 0) + z[0] * tr_coeffs(1, 1) + z[0] * z[0] * tr_coeffs(1, 2);
            auto p01 = tr_coeffs(2, 0) + z[0] * tr_coeffs(2, 1) + z[0] * z[0] * tr_coeffs(2, 2);
            auto p20 = tr_coeffs(3, 0) + z[0] * tr_coeffs(3, 1) + z[0] * z[0] * tr_coeffs(3, 2);
            auto p11 = tr_coeffs(4, 0) + z[0] * tr_coeffs(4, 1) + z[0] * z[0] * tr_coeffs(4, 2);
            auto p02 = tr_coeffs(5, 0) + z[0] * tr_coeffs(5, 1) + z[0] * z[0] * tr_coeffs(5, 2);

            auto dc_scale = 1.0 / (vc[0] * 1E3 + vc[1] * 1E3) * 1E3;
            auto tc_scale = sqrt(Tc[0] * Tc[1]);
            auto tau = T / tc_scale;
            auto delta = rho / dc_scale;

            return pow(z[0], 2.0) * Tc[0] + pow(z[1], 2.0) * Tc[1] - 2.0 * z[0] * z[1] * \
                (p00 + p10 * delta + p01 * tau + p20 * delta * delta + p02 * tau * tau + p11 * delta * tau) * tc_scale;
        }


        template <typename TTYPE, typename RHOTYPE, typename MoleFractions>
        auto get_dr(const TTYPE& T, const RHOTYPE& rho, const MoleFractions& z) const {

            auto p00 = dr_coeffs(0, 0) + z[0] * dr_coeffs(0, 1) + z[0] * z[0] * dr_coeffs(0, 2);
            auto p10 = dr_coeffs(1, 0) + z[0] * dr_coeffs(1, 1) + z[0] * z[0] * dr_coeffs(1, 2);
            auto p01 = dr_coeffs(2, 0) + z[0] * dr_coeffs(2, 1) + z[0] * z[0] * dr_coeffs(2, 2);
            auto p20 = dr_coeffs(3, 0) + z[0] * dr_coeffs(3, 1) + z[0] * z[0] * dr_coeffs(3, 2);
            auto p11 = dr_coeffs(4, 0) + z[0] * dr_coeffs(4, 1) + z[0] * z[0] * dr_coeffs(4, 2);
            auto p02 = dr_coeffs(5, 0) + z[0] * dr_coeffs(5, 1) + z[0] * z[0] * dr_coeffs(5, 2);
            auto dc_scale = 1.0 / (vc[0] * 1E3 + vc[1] * 1E3) * 1E3;
            auto tc_scale = sqrt(Tc[0] * Tc[1]);
            auto tau = T / tc_scale;
            auto delta = rho / dc_scale;

            auto vc_ = pow(z[0], 2.0) * vc[0] * 1E3 + pow(z[1], 2.0) * vc[1] * 1E3 - 2.0 * z[0] * z[1] * \
                (p00 + p10 * delta + p01 * tau + p20 * delta * delta + p02 * tau * tau + p11 * delta * tau) * dc_scale / 1E3;
            return 1E3 * (1.0 / vc_);
        }

    };

    template<typename BaseClass>
    class MultiFluidAdapter_Ecs {

    private:
        std::string meta = "";

    public:
        const BaseClass& base;
        const Reducing_ECS redfunc;

        template<class VecType>
        auto R(const VecType& molefrac) const { return base.R(molefrac); }

        MultiFluidAdapter_Ecs(const BaseClass& base, Reducing_ECS&& redfunc) : base(base), redfunc(redfunc) {};

        /// Store some sort of metadata in string form (perhaps a JSON representation of the model?)
        void set_meta(const std::string& m) { meta = m; }
        /// Get the metadata stored in string form
        auto get_meta() const { return meta; }

        template<typename TType, typename RhoType, typename MoleFracType>
        auto alphar(const TType& T,
            const RhoType& rho,
            const MoleFracType& molefrac) const
        {
            auto Tred = forceeval(redfunc.get_tr(T, rho, molefrac));
            auto rhored = forceeval(redfunc.get_dr(T, rho, molefrac));
            auto delta = forceeval(rho / rhored);
            auto tau = forceeval(Tred / T);
            auto val = base.corr.alphar(tau, delta, molefrac);
            return forceeval(val);
        }
    };

    template<class Model>
    auto build_multifluid_mutant_ecs(const Model& model, const nlohmann::json& jj) {
        auto N = model.redfunc.Tc.size();
        auto red = model.redfunc;
        auto Tc = red.Tc, vc = red.vc;
        auto newred = Reducing_ECS(Tc, vc, jj);
        auto mfa = MultiFluidAdapter_Ecs(model, std::move(newred));
        mfa.set_meta(jj.dump());
        return mfa;
    }

}