Mixtures work now too for alphaig

Just need to fill in the remaining terms and test, test, test

include/teqp/ideal_eosterms.hpp

@@ -30,11 +30,11 @@ namespace teqp {
 
     using IdealHelmholtzTerms = std::variant<IdealHelmholtzLead>;
 
-    class IdealHelmholtz {
+    class PureIdealHelmholtz {
     public:
         std::vector<IdealHelmholtzTerms> contributions;
-        IdealHelmholtz(const nlohmann::json& j) {
-            for (auto& term : j) {
+        PureIdealHelmholtz(const nlohmann::json& jpure) {
+            for (auto& term : jpure) {
                 if (term.at("type") == "Lead") {
                     contributions.emplace_back(IdealHelmholtzLead(term.at("a_1"), term.at("a_2")));
                 }
@@ -43,8 +43,8 @@ namespace teqp {
                 }
             }
         }
-        template<typename TType, typename RhoType, typename MoleFrac>
-        auto alphaig(const TType& T, const RhoType &rho, const MoleFrac& /**/) const{
+        template<typename TType, typename RhoType>
+        auto alphaig(const TType& T, const RhoType &rho) const{
             std::common_type_t <TType, RhoType> ig = 0.0;
             for (const auto& term : contributions) {
                 auto contrib = std::visit([&](auto& t) { return t.alphaig(T, rho); }, term);
@@ -54,4 +54,43 @@ namespace teqp {
         }
     };
 
+    /**
+     * @brief Ideal-gas Helmholtz energy container
+     * 
+     * \f[ \alpha^{\rm ig} = \sum_i x_i[\alpha^{\rm ig}_{oi}(T,\rho) + x_i] \f]
+     *
+     * where x_i are mole fractions
+     * 
+     */
+    class IdealHelmholtz {
+        
+        public:
+        
+        std::vector<PureIdealHelmholtz> pures;
+        
+        IdealHelmholtz(const nlohmann::json &jpures){
+            for (auto &jpure : jpures){
+                pures.emplace_back(jpure);
+            }
+        }
+
+        template<typename TType, typename RhoType, typename MoleFrac>
+        auto alphaig(const TType& T, const RhoType &rho, const MoleFrac &molefrac) const {
+            std::common_type_t <TType, RhoType, decltype(molefrac[0])> ig = 0.0;
+            if (molefrac.size() != pures.size()){
+                throw teqp::InvalidArgument("molefrac and pures are not the same length");
+            }
+            std::size_t i = 0;
+            for (auto &pure : pures){
+                if (molefrac[i] != 0){
+                    ig += molefrac[i]*(pure.alphaig(T, rho) + log(molefrac[i]));
+                }
+                else{
+                    // lim_{x\to 0} x*ln(x) => 0
+                }
+                i++;
+            }
+            return ig;
+        }
+    };
 }
\ No newline at end of file

src/tests/catch_test_alphaig.cxx

@@ -10,8 +10,10 @@ using namespace teqp;
 
 TEST_CASE("Simplest case","[alphaig]") {
 	double a_1 = 1, a_2 = 2, T = 300;
+	nlohmann::json j0 = nlohmann::json::array();
+	j0.push_back({ {"type", "Lead"}, { "a_1", 1 }, { "a_2", 2 } });
 	nlohmann::json j = nlohmann::json::array();
-	j.push_back({ {"type", "Lead"}, { "a_1", 1 }, { "a_2", 2 } });
+	j.push_back(j0);
 	IdealHelmholtz ih(j);
 	std::valarray<double> molefrac{1.0};
 	REQUIRE(ih.alphaig(T, 1, molefrac) == log(1) + a_1 + a_2 / T);
@@ -19,8 +21,10 @@ TEST_CASE("Simplest case","[alphaig]") {
 
 TEST_CASE("alphaig derivative", "[alphaig]") {
 	double a_1 = 1, a_2 = 2, T = 300, rho = 1;
+	nlohmann::json j0 = nlohmann::json::array();
+	j0.push_back({ {"type", "Lead"}, { "a_1", 1 }, { "a_2", 2 } }); // For the first component
 	nlohmann::json j = nlohmann::json::array();
-	j.push_back({ {"type", "Lead"}, { "a_1", 1 }, { "a_2", 2 } });
+	j.push_back(j0);
 	IdealHelmholtz ih(j);
 	auto molefrac = (Eigen::ArrayXd(1) << 1.0).finished();
 	auto wih = AlphaCallWrapper<1, decltype(ih)>(ih);