diff --git a/include/teqp/models/multifluid_ancillaries.hpp b/include/teqp/models/multifluid_ancillaries.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5b3eb8f01c5c6db0272a7ffd9246e6421a14afb2
--- /dev/null
+++ b/include/teqp/models/multifluid_ancillaries.hpp
@@ -0,0 +1,49 @@
+#pragma once
+
+#include <valarray>
+#include "nlohmann/json.hpp"
+
+namespace teqp{
+
+	struct VLEAncillary{
+		const double T_r, Tmax, Tmin, reducing_value;
+		const std::string type, description;
+		const std::valarray<double> n, t;
+		const bool using_tau_r, noexp;
+	
+		VLEAncillary(const nlohmann::json &j) :
+			T_r(j.at("T_r")),
+            Tmax(j.at("Tmax")),
+            Tmin(j.at("Tmin")),
+            description(j.at("description")),
+            n(j.at("n").get<std::valarray<double>>()),
+            t(j.at("t").get<std::valarray<double>>()),
+            reducing_value(j.at("reducing_value")),
+            type(j.at("type")),
+            using_tau_r(j.at("using_tau_r")),
+            noexp(type == "rhoLnoexp"){};
+
+		double operator() (double T) const{
+			auto Theta = 1-T/T_r;
+			auto RHS = (pow(Theta, t)*n).sum();
+			if (using_tau_r){
+				RHS *= T_r/T;
+			}
+			if (noexp){
+	            return reducing_value*(1+RHS);
+			}
+	        else{
+	            return exp(RHS)*reducing_value;
+	        }
+		};
+	};
+
+	struct MultiFluidVLEAncillaries {
+		const VLEAncillary rhoL, rhoV, pL, pV;
+		MultiFluidVLEAncillaries(const nlohmann::json& j) :
+			rhoL(VLEAncillary(j.at("rhoL"))),
+			rhoV(VLEAncillary(j.at("rhoV"))),
+			pL((j.contains("pS")) ? VLEAncillary(j.at("pS")) : VLEAncillary(j.at("pL"))),
+			pV((j.contains("pS")) ? VLEAncillary(j.at("pS")) : VLEAncillary(j.at("pV"))){}
+	};
+}
\ No newline at end of file
diff --git a/interface/multifluid.cpp b/interface/multifluid.cpp
index 4f09fa49ad56a74514ed0c4f29be0b081c5bd3be..76e1089d5d9926d99c86cff78e3d50538ea001ef 100644
--- a/interface/multifluid.cpp
+++ b/interface/multifluid.cpp
@@ -1,6 +1,7 @@
 #include "pybind11_wrapper.hpp"
 
 #include "teqp/models/multifluid.hpp"
+#include "teqp/models/multifluid_ancillaries.hpp"
 #include "teqp/derivs.hpp"
 #include "teqp/models/mie/lennardjones.hpp"
 
@@ -8,6 +9,24 @@
 
 void add_multifluid(py::module& m) {
 
+    // A single ancillary curve
+    py::class_<VLEAncillary>(m, "VLEAncillary")
+        .def(py::init<const nlohmann::json&>())
+        .def("__call__", &VLEAncillary::operator())
+        .def_readonly("T_r", &VLEAncillary::T_r)
+        .def_readonly("Tmax", &VLEAncillary::Tmax)
+        .def_readonly("Tmin", &VLEAncillary::Tmin)
+        ;
+
+    // The collection of VLE ancillary curves
+    py::class_<MultiFluidVLEAncillaries>(m, "MultiFluidVLEAncillaries")
+        .def(py::init<const nlohmann::json&>())
+        .def_readonly("rhoL", &MultiFluidVLEAncillaries::rhoL)
+        .def_readonly("rhoV", &MultiFluidVLEAncillaries::rhoV)
+        .def_readonly("pL", &MultiFluidVLEAncillaries::pL)
+        .def_readonly("pV", &MultiFluidVLEAncillaries::pV)
+        ;
+
     // Multifluid model
     m.def("build_multifluid_model", &build_multifluid_model, py::arg("components"), py::arg("coolprop_root"), py::arg("BIPcollectionpath") = "", py::arg("flags") = nlohmann::json{}, py::arg("departurepath") = "");
     using MultiFluid = decltype(build_multifluid_model(std::vector<std::string>{"",""},"",""));
diff --git a/interface/multifluid_shared.hpp b/interface/multifluid_shared.hpp
index 2aa39f9be8851819ba3f72ead9a9d3bb7e1b4cb8..3830b031ab8a952e58aeca269229105e16d7a8d7 100644
--- a/interface/multifluid_shared.hpp
+++ b/interface/multifluid_shared.hpp
@@ -1,5 +1,7 @@
 #pragma once 
 
+#include "teqp/exceptions.hpp"
+
 template<typename Model, typename Wrapper>
 void add_multifluid_methods(Wrapper &wMF){
     wMF.def("get_Tcvec", [](const Model& c) { return c.redfunc.Tc; })
@@ -7,5 +9,13 @@ void add_multifluid_methods(Wrapper &wMF){
        .def("get_Tr", [](const Model& c, const Eigen::ArrayXd &molefrac) { return c.redfunc.get_Tr(molefrac); })
        .def("get_rhor", [](const Model& c, const Eigen::ArrayXd &molefrac) { return c.redfunc.get_rhor(molefrac); })
        .def("set_meta", [](Model& c, const std::string &s) { return c.set_meta(s); })
-       .def("get_meta", [](const Model& c) { return c.get_meta(); });
+       .def("get_meta", [](const Model& c) { return c.get_meta(); })
+       .def("build_ancillaries", [](const Model& c) { 
+           if (c.redfunc.Tc.size() != 1) {
+               throw teqp::InvalidArgument("Can only build ancillaries for pure fluids");
+           }
+           auto jancillaries = nlohmann::json::parse(c.get_meta()).at("pures")[0].at("ANCILLARIES");
+           return teqp::MultiFluidVLEAncillaries(jancillaries);
+           })
+       ;
 }
\ No newline at end of file
diff --git a/src/test_multifluid_ancillaries.cpp b/src/test_multifluid_ancillaries.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7cecb4f0ddce9fa6f9b82f6427bbe978c6a8d6ae
--- /dev/null
+++ b/src/test_multifluid_ancillaries.cpp
@@ -0,0 +1,16 @@
+
+#include "teqp/models/multifluid_ancillaries.hpp"
+#include "teqp/models/multifluid.hpp"
+#include "teqp/algorithms/VLE.hpp"
+
+int main() {
+	auto model = teqp::build_multifluid_model({"n-Propane"}, "../mycp");
+	auto j = nlohmann::json::parse(model.get_meta());
+	auto jancillaries = j.at("pures")[0].at("ANCILLARIES");
+	teqp::MultiFluidVLEAncillaries anc(jancillaries);
+	double T = 340;
+	auto rhoV = anc.rhoV(T), rhoL = anc.rhoL(T);
+	auto rhovec = teqp::pure_VLE_T(model, T, rhoL, rhoV, 10);
+	//auto [rhoLnew, rhoVnew] = rhovec;
+	int rr = 0;
+}
\ No newline at end of file
diff --git a/src/tests/catch_test_multifluid.cxx b/src/tests/catch_test_multifluid.cxx
index 473b09036f923f590a85b3938d93e437ac49e31f..c78ba0205925ef19a8f46d26633e9f5bfdafcdf3 100644
--- a/src/tests/catch_test_multifluid.cxx
+++ b/src/tests/catch_test_multifluid.cxx
@@ -4,6 +4,7 @@
 using Catch::Approx;
 
 #include "teqp/models/multifluid.hpp"
+#include "teqp/models/multifluid_ancillaries.hpp"
 #include "teqp/algorithms/critical_tracing.hpp"
 #include "teqp/algorithms/VLE.hpp"
 #include "teqp/filesystem.hpp"
@@ -139,6 +140,26 @@ TEST_CASE("Check that all pure fluid models can be instantiated", "[multifluid],
     }    
 }
 
+TEST_CASE("Check that all ancillaries can be instantiated and work properly", "[multifluid],[all]") {
+    std::string root = "../mycp";
+    SECTION("With absolute paths to json file") {
+        int counter = 0;
+        for (auto path : get_files_in_folder(root + "/dev/fluids", ".json")) {
+            if (path.filename().stem() == "Methanol") { continue; }
+            CAPTURE(path.string());
+            auto abspath = std::filesystem::absolute(path).string();
+            auto model = build_multifluid_model({ abspath }, root, root + "/dev/mixtures/mixture_binary_pairs.json");
+            auto jancillaries = nlohmann::json::parse(model.get_meta()).at("pures")[0].at("ANCILLARIES");
+            auto anc = teqp::MultiFluidVLEAncillaries(jancillaries);
+            double T = 0.9*anc.rhoL.T_r;
+            auto rhoV = anc.rhoV(T), rhoL = anc.rhoL(T);
+            auto rhovec = teqp::pure_VLE_T(model, T, rhoL, rhoV, 10);
+            counter += 1;
+        }
+        CHECK(counter > 100);
+    }
+}
+
 TEST_CASE("Check that mixtures can also do absolute paths", "[multifluid],[abspath]") {
     std::string root = "../mycp";
     SECTION("With absolute paths to json file") {