multifluid.hpp

        throw std::invalid_argument("Only binary mixtures are currently supported with invariant departure functions");
    }

    auto phiT = red.betaT, lambdaT = red.gammaT, phiV = red.betaV, lambdaV = red.gammaV;
    phiT.setOnes(); phiV.setOnes();
    lambdaV.setZero(); lambdaV.setZero();
    auto Tc = red.Tc, vc = red.vc;

    // Allocate the matrices of default models and F factors
    Eigen::MatrixXd F(N, N); F.setZero();
    std::vector<std::vector<DepartureTerms>> funcs(N);
    for (auto i = 0; i < N; ++i) { funcs[i].resize(N); }

    for (auto i = 0; i < N; ++i) {
        for (auto j = i; j < N; ++j) {
            if (i == j) {
                funcs[i][i].add_term(NullEOSTerm());
            }
            else {
                // Extract the given entry
                auto entry = jj[std::to_string(i)][std::to_string(j)];

                auto BIP = entry["BIP"];
                // Set the reducing function parameters in the copy
                phiT(i, j) = BIP.at("phiT");
                phiT(j, i) = phiT(i, j);
                lambdaT(i, j) = BIP.at("lambdaT"); 
                lambdaT(j, i) = -lambdaT(i, j);

                phiV(i, j) = BIP.at("phiV");
                phiV(j, i) = phiV(i, j);
                lambdaV(i, j) = BIP.at("lambdaV"); 
                lambdaV(j, i) = -lambdaV(i, j);

                // Build the matrix of F and departure functions
                auto dep = entry["departure"];
                F(i, j) = BIP.at("Fij");
                F(j, i) = F(i, j);
                funcs[i][j] = build_departure_function(dep);
                funcs[j][i] = build_departure_function(dep);
            }
        }
    }

    auto newred = MultiFluidInvariantReducingFunction(phiT, lambdaT, phiV, lambdaV, Tc, vc);
    auto newdep = DepartureContribution(std::move(F), std::move(funcs));
    auto mfa = MultiFluidAdapter(model, std::move(newred), std::move(newdep));
    /// Store the departure term in the adapted multifluid class
    mfa.set_meta(jj.dump());
    return mfa;
}


class DummyEOS {
public:
    template<typename TType, typename RhoType> auto alphar(TType tau, const RhoType& delta) const { return tau * delta; }
};
class DummyReducingFunction {
public:
    template<typename MoleFractions> auto get_Tr(const MoleFractions& molefracs) const { return molefracs[0]; }
    template<typename MoleFractions> auto get_rhor(const MoleFractions& molefracs) const { return molefracs[0]; }
};
inline auto build_dummy_multifluid_model(const std::vector<std::string>& components) {
    std::vector<DummyEOS> EOSs(2);
    std::vector<std::vector<DummyEOS>> funcs(2); for (auto i = 0; i < funcs.size(); ++i) { funcs[i].resize(funcs.size()); }
    std::vector<std::vector<double>> F(2); for (auto i = 0; i < F.size(); ++i) { F[i].resize(F.size()); }

    struct Fwrapper {
    private: 
        const std::vector<std::vector<double>> F_;
    public:
        Fwrapper(const std::vector<std::vector<double>> &F) : F_(F){};
        auto operator ()(std::size_t i, std::size_t j) const{ return F_[i][j]; }
    };
    auto ff = Fwrapper(F);
    auto redfunc = DummyReducingFunction();
    return MultiFluid(std::move(redfunc), std::move(CorrespondingStatesContribution(std::move(EOSs))), std::move(DepartureContribution(std::move(ff), std::move(funcs))));
}
inline void test_dummy() {
    auto model = build_dummy_multifluid_model({ "A", "B" });
    std::valarray<double> rhovec = { 1.0, 2.0 };
    auto alphar = model.alphar(300.0, rhovec);
}