From 184ada1b7fa287d9c08de9b6dd4440f20b659510 Mon Sep 17 00:00:00 2001
From: Sven <sven.pohl11991@gmail.com>
Date: Tue, 8 Aug 2023 17:45:15 +0200
Subject: [PATCH] Added properties.hpp for statistics
---
include/teqp/algorithms/properties.hpp | 639 +++++++++++++
include/teqp/algorithms/vle_trend.hpp | 558 ++++++++++++
include/teqp/models/mie/mie.hpp | 1164 ++++++++++++++++++++++--
include/teqp/types.hpp | 6 +
interface/CPP/teqp_impl_factory.cpp | 2 +-
5 files changed, 2305 insertions(+), 64 deletions(-)
create mode 100644 include/teqp/algorithms/properties.hpp
create mode 100644 include/teqp/algorithms/vle_trend.hpp
diff --git a/include/teqp/algorithms/properties.hpp b/include/teqp/algorithms/properties.hpp
new file mode 100644
index 0000000..5813fce
--- /dev/null
+++ b/include/teqp/algorithms/properties.hpp
@@ -0,0 +1,639 @@
+#pragma once
+#pragma once
+#include "teqp/types.hpp"
+#include "teqp/derivs.hpp"
+#include "teqp/algorithms/VLE.hpp"
+#include "teqp/models/multifluid_ancillaries.hpp"
+#include <Eigen/Dense>
+#include <fstream>
+#include <vector>
+#include <iostream>
+#include <string>
+#include <algorithm>
+#include <teqp/models/cubics.hpp>
+#include "teqp/algorithms/vle_trend.hpp"
+#include "teqp/cpp/teqpcpp.hpp"
+#include "teqp/json_builder.hpp"
+#include <teqp/cpp/deriv_adapter.hpp>
+#include "teqp/ideal_eosterms.hpp"
+using namespace std;
+using namespace teqp;
+using namespace Mie;
+
+namespace teqp {
+ namespace cppinterface {
+
+ std::unique_ptr<teqp::cppinterface::AbstractModel> make_SAFTVRMie(const nlohmann::json& j);
+
+ using makefunc = std::function<std::unique_ptr<teqp::cppinterface::AbstractModel>(const nlohmann::json& j)>;
+ using namespace teqp::cppinterface::adapter;
+
+ static std::unordered_map<std::string, makefunc> pointer_factory = {
+ {"MieElong", [](const nlohmann::json& spec) { return make_owned(Mie::MieElong(spec.at("model_path"),spec.at("components"),spec.at("combination"))); }}
+ };
+
+ std::unique_ptr<teqp::cppinterface::AbstractModel> build_model_ptr(const nlohmann::json& json) {
+
+ // Extract the name of the model and the model parameters
+ std::string kind = json.at("kind");
+ auto spec = json.at("model");
+
+ auto itr = pointer_factory.find(kind);
+ if (itr != pointer_factory.end()) {
+ return (itr->second)(spec);
+ }
+ else {
+ throw std::invalid_argument("Don't understand \"kind\" of: " + kind);
+ }
+ }
+
+ std::unique_ptr<AbstractModel> make_multifluid_model(const std::vector<std::string>& components, const std::string& coolprop_root, const std::string& BIPcollectionpath, const nlohmann::json& flags, const std::string& departurepath) {
+ return make_owned(build_multifluid_model(components, coolprop_root, BIPcollectionpath, flags, departurepath));
+ }
+
+ std::unique_ptr<AbstractModel> make_model(const nlohmann::json& j) {
+ return build_model_ptr(j);
+ }
+ }
+}
+
+namespace teqp {
+
+
+ struct data_point {
+ double T, P, D, SND;
+ std::vector<double> x;
+ std::string author, mix_name;
+ int nr;
+ int id;
+ double dense_v, dense_l, snd_l, snd_v;
+ std::vector<double> dense_v_pure;
+ std::vector<double> dense_l_pure;
+ double x1, x2, y1, y2;
+ // excess volumina of vapor and liquid phase
+ double ve_v, ve_l;
+ double deviation;
+ };
+
+
+
+ inline auto get_compare_settings(std::string& path) {
+ std::string filepath = std::filesystem::is_regular_file(path) ? path : path;
+ nlohmann::json j = load_a_JSON_file(filepath);
+ compare_settings set;
+ auto spec = j.at("model");
+ set.model_path = spec.at("model_path");
+ set.srk_path = spec.at("srk_path");
+ set.BIPcollection = spec.at("BIP");
+ set.depcollection = spec.at("departure");
+ set.coolprop_root = spec.at("coolprop_root");
+ std::vector<std::string> mixture = spec.at("components");
+ std::vector<std::string> dev_files = spec.at("dev_files");
+ set.dev_files = dev_files;
+ set.mix = mixture;
+ set.combination = spec.at("combination");
+ set.dev_dir = spec.at("dev_dir");
+ return set;
+ }
+
+
+ template <typename ModelType>
+ inline auto check_phase(const ModelType& model, const double& T, const double& rho, const std::vector<double>& molefrac, teqp::phase& phase_est) {
+
+ // Check if pressure is positive (phase doesnt matter)
+ auto is_correct_phase = false;
+ using tdx = TDXDerivatives<ModelType, double, decltype(molefrac)>;
+ auto p = 0.0; auto dpdd = 0.0; auto d2pdd2 = 0.0;
+ auto ar01 = 0.0; auto ar02 = 0.0; auto ar03 = 0.0;
+ ar01 = tdx::get_Ar01(model, T, rho, molefrac);
+ ar02 = tdx::get_Ar02(model, T, rho, molefrac);
+ ar03 = tdx::get_Ar03(model, T, rho, molefrac);
+ p = (1.0 + ar01) * T * rho * model.R(molefrac);
+ dpdd = (1.0 + 2.0 * ar01 + ar02) * model.R(molefrac) * T;
+ d2pdd2 = (2.0 * ar01 + 4.0 * ar02 + ar03) * model.R(molefrac) * T / rho;
+
+
+ if (p > 0.0 && dpdd > 0.0 && d2pdd2 < 0.0 && phase_est == teqp::VAPOR) {
+ is_correct_phase = true;
+ }
+
+ if (!(is_correct_phase) && p > 0.0 && dpdd > 0.0 && d2pdd2 > 0.0 && phase_est == teqp::LIQUID) {
+ is_correct_phase = true;
+ }
+
+ return is_correct_phase;
+ }
+
+ template <typename ModelType>
+ std::tuple<double, double, double, double> get_rho_bounds(const ModelType& model, const double& T, const double& rho_est, const std::vector<double>& molefrac) {
+
+ using tdx = TDXDerivatives<ModelType, double, decltype(molefrac)>;
+ auto rhomax = 1.100 * rho_est;
+ auto rhomin = 0.900 * rho_est;
+ auto rhomax_allowed = 2.0 * rho_est;
+ auto rhomin_allowed = 0.70 * rho_est;
+ auto ptest = 0.0;
+ ptest = (1.0 + tdx::get_Ar01(model, T, rho_est, molefrac)) * T * rho_est * model.R(molefrac);
+
+ if (ptest < 1E-12) {
+ rhomax = (rhomin_allowed + rhomax_allowed) / 2.0;
+ rhomin = rhomax * 0.80;
+ }
+ else if (isnan(ptest)) {
+ rhomax = rhomax * 0.70;
+ rhomin = rhomax * 0.50;
+ rhomin_allowed = rhomin * 0.10;
+ }
+ return std::make_tuple(rhomin, rhomax, rhomin_allowed, rhomax_allowed);
+ }
+
+ template <typename ModelType>
+ std::tuple<double, double, double, double> get_rho_bounds_extended(const ModelType& model, const double& T, const double& rho_est, const std::vector<double>& molefrac) {
+
+ using tdx = TDXDerivatives<ModelType, double, decltype(molefrac)>;
+ auto rhomax = 2.00 * rho_est;
+ auto rhomin = 0.500 * rho_est;
+ auto rhomax_allowed = 3.0 * rho_est;
+ auto rhomin_allowed = 0.2 * rho_est;
+ auto ptest = 0.0;
+ ptest = (1.0 + tdx::get_Ar01(model, T, rho_est, molefrac)) * T * rho_est * model.R(molefrac);
+
+ if (ptest < 1E-12) {
+ rhomax = (rhomin_allowed + rhomax_allowed) / 2.0;
+ rhomin = rhomax * 0.80;
+ }
+ else if (isnan(ptest)) {
+ rhomax = rhomax * 0.70;
+ rhomin = rhomax * 0.50;
+ rhomin_allowed = rhomin * 0.10;
+ }
+ return std::make_tuple(rhomin, rhomax, rhomin_allowed, rhomax_allowed);
+ }
+
+
+
+ template <typename ModelType, typename ScalarT, typename ScalarP>
+ auto get_rho_from_t_and_p(const ModelType& model, psrk& psrk_s, ScalarT& T, ScalarP& P, const std::vector<double>& molefrac) {
+
+
+ auto rho_est = 0.0;
+ using tdx = TDXDerivatives<ModelType, ScalarT, decltype(molefrac)>;
+
+ // Assume phase is liquid
+ teqp::phase phase = LIQUID;
+
+ // Get estimate from pressure srk model
+ rho_est = psrk_s.get_density(T, P, molefrac, model.R(molefrac), phase);
+
+
+ // Define deviation pressure function
+ auto p_diff = [T, P, molefrac, model](double rho) {
+ auto A01 = 0.0;
+ A01 = tdx::get_Ar01(model, T, rho, molefrac);
+ return P - (1.0 + A01) * T * rho * model.R(molefrac);
+ };
+
+ // Initialize the iteration parameters
+ auto iter = 0;
+ auto max_iter = 50;
+
+ // Get bounds for regular falsi
+ auto b = get_rho_bounds(model, T, rho_est, molefrac);
+ auto rho_l = 0.0;
+ rho_l = regular_falsi(p_diff, std::get<0>(b), std::get<1>(b), 1E-8, std::get<2>(b), std::get<3>(b), max_iter, iter);
+ auto is_correct = check_phase(model, T, rho_l, molefrac, phase);
+
+ //if (!(is_correct)) {
+ phase = VAPOR;
+ rho_est = psrk_s.get_density(T, P, molefrac, model.R(molefrac), phase);
+ b = get_rho_bounds(model, T, rho_est, molefrac);
+ auto rho_v = 0.0;
+ rho_v = regular_falsi(p_diff, std::get<0>(b), std::get<1>(b), 1E-8, std::get<2>(b), std::get<3>(b), max_iter, iter);
+ is_correct = check_phase(model, T, rho_v, molefrac, phase);
+ if (!(is_correct)) {
+ b = get_rho_bounds_extended(model, T, rho_est, molefrac);
+ rho_v = regular_falsi(p_diff, std::get<0>(b), std::get<1>(b), 1E-8, std::get<2>(b), std::get<3>(b), max_iter, iter);
+ is_correct = check_phase(model, T, rho_v, molefrac, phase);
+ }
+ //}
+
+ return std::make_tuple(rho_v, rho_l);
+ }
+
+ auto get_fluid_model(compare_settings& s) {
+ /* if (s.model_type.compare("Pohl") == 0) {*/
+ const auto mix_residual = MieElong(s.model_path, s.mix, s.combination);
+
+ // get all pures a seperate fluids
+ std::vector<MieElong> pures_residual;
+ for (size_t i = 0; i < s.mix.size(); i++) {
+ pures_residual.push_back(MieElong(s.model_path, { s.mix[i] }, s.combination));
+ }
+
+ auto pure_data = nlohmann::json::array();
+ for (auto f : s.mix) {
+ std::string path = s.coolprop_root + "/dev/fluids/" + f + ".json";
+ auto j = convert_CoolProp_idealgas(path, 0);
+ pure_data.push_back(j);
+ }
+
+ IdealHelmholtz mix_ideal(pure_data);
+
+ return std::make_tuple(pures_residual, mix_residual, mix_ideal);
+ //return
+// //nlohmann::json fluid_i = load_a_JSON_file(s.model_path);
+
+ //// get ideal contribution from Multifluid coolprop files
+ //pure_data = ideal_term(s.mix, fluid_i);
+
+ //}
+ //if (s.model_type.compare("Multi-Fluid") == 0) {
+ //const auto mixture_model = build_multifluid_model(s.mix, s.coolprop_root, s.BIPcollection, s.depcollection);
+ //auto jj = nlohmann::json::parse(mixture_model.get_meta());
+ //auto ancillaries = jj.at("pures")[0].at("ANCILLARIES");
+ //teqp::MultiFluidVLEAncillaries anc(ancillaries);
+ ////return build_multifluid_model({ mixture[1] }, coolprop_root, BIPcollection, depcollection);
+ //return mixture_model;
+ //}
+ }
+
+ struct vle_info
+ {
+ std::vector<double> pvap, pliq, xL_0, xV_0;
+ };
+
+
+ class data_set {
+
+ private:
+ std::vector<data_point> points;
+ std::vector<string> head;
+ prop_t prop;
+ std::vector<double> isotherms;
+ std::vector<vle_info> vle;
+ public:
+
+
+ // get dev file PVT for mixtures and add to data points
+ void get_pvt_mix_line(std::string& line) {
+ std::stringstream ss;
+ data_point p;
+ double x1 = 0.0; double x2 = 0.0;
+ std::string name_mix, add_info;
+ ss << line;
+ ss >> p.P >> p.D >> p.T >> x1 >> x2;
+ std::stringstream rest_line;
+ if (line.length() > 6) {
+ add_info = line.substr(51, line.length());
+ rest_line << add_info.substr(7, add_info.length());
+ p.author = add_info.substr(0, 7);
+ rest_line >> p.id >> p.mix_name >> p.nr;
+ p.x.push_back(x1);
+ p.x.push_back(x2);
+ if (p.P > 0.0 || (p.D > 0.0 && p.T > 0.0)) {
+ points.push_back(p);
+ }
+ }
+
+ }
+
+ void get_snd_mix_line(std::string& line) {
+ std::stringstream ss;
+ data_point p;
+ double x1 = 0.0; double x2 = 0.0;
+ std::string name_mix, add_info;
+ ss << line;
+ ss >> p.P >> p.T >> p.SND >> x1 >> x2;
+ std::stringstream rest_line;
+ if (line.length() > 6) {
+ add_info = line.substr(62, line.length());
+ rest_line << add_info.substr(7, add_info.length());
+ p.author = add_info.substr(0, 7);
+ rest_line >> p.id >> p.mix_name >> p.nr;
+ p.x.push_back(x1);
+ p.x.push_back(x2);
+ if (p.P > 0.0 || (p.SND > 0.0 && p.T > 0.0)) {
+ points.push_back(p);
+ }
+ }
+ }
+
+ void get_vle_mix_line(std::string& line) {
+ std::stringstream ss;
+ data_point p;
+ double x1 = 0.0; double x2 = 0.0;
+ std::string name_mix, add_info;
+ ss << line;
+ ss >> p.P >> p.T >> p.x1 >> p.x2 >> p.y1 >> p.y2;
+ std::stringstream rest_line;
+ if (line.length() > 6) {
+ add_info = line.substr(51, line.length());
+ rest_line << add_info.substr(7, add_info.length());
+ p.author = add_info.substr(0, 7);
+ rest_line >> p.id >> p.mix_name >> p.nr;
+ p.x.push_back(x1);
+ p.x.push_back(x2);
+ if (p.P > 0.0 || (p.SND > 0.0 && p.T > 0.0)) {
+ points.push_back(p);
+ }
+ }
+ }
+
+ // read dev file by path and assign data
+ void read_dev_file(std::string& path) {
+ std::ifstream file;
+ std::string line;
+ file.open(path);
+ bool head_end = false;
+
+ prop = PVT;
+
+ if (path.find("PVT") != std::string::npos) {
+ prop = PVT;
+ }
+ else if (path.find("VE") != std::string::npos) {
+ prop = VE;
+ }
+ else if (path.find("SND") != std::string::npos) {
+ prop = SND;
+ }
+ else if (path.find("VLE") != std::string::npos) {
+ prop = VLE;
+ }
+
+
+ while (std::getline(file, line)) {
+ if (line.find("@PCL5") != std::string::npos) {
+ head.push_back(line);
+ head_end = true;
+ }
+ else if (!(head_end)) {
+ head.push_back(line);
+ }
+ else // numbers
+ {
+ if (prop == PVT) {
+ get_pvt_mix_line(line);
+ }
+ else if (prop == VE) {
+ get_pvt_mix_line(line);
+ }
+ else if (prop == SND) {
+ get_snd_mix_line(line);
+ }
+ else if (prop == VLE) {
+ get_vle_mix_line(line);
+ }
+ }
+ }
+ file.close();
+
+ // get isotherms
+ std::vector<double> tmp;
+ std::transform(points.begin(), points.end(), std::back_inserter(tmp),
+ [](data_point const& p) -> double { return p.T; });
+
+ for (size_t i = 0; i < tmp.size(); ++i) {
+ if (count(isotherms.begin(), isotherms.end(), tmp.at(i)) == 0) {
+ isotherms.push_back(tmp.at(i));
+ }
+ }
+
+
+
+ }
+
+ // get the p-x relation for all isotherms in the dev file
+ template <typename ModelType, typename PureType>
+ void get_binary_vle(ModelType& mixture, PureType& pure) {
+
+ for (auto t : isotherms) {
+ auto [rho_vap, rho_liq] = get_vle_pure_start(pure[0], t, 0);
+ auto vle_pure_1 = pure_VLE_T(pure[0], t, 1E3 * rho_liq, 1E3 * rho_vap, 50);
+ auto rhovecV0 = (Eigen::ArrayXd(2) << vle_pure_1(1), 0.0).finished();
+ auto rhovecL0 = (Eigen::ArrayXd(2) << vle_pure_1(0), 0.0).finished();
+ auto res = trace_VLE_isotherm_binary(mixture, t, rhovecL0, rhovecV0);
+
+ auto vec_packed = [](const auto& res, const auto& key) {
+ std::vector<double> x;
+ for (auto res_point : res) {
+ x.push_back(res_point[key]);
+ }
+ return x;
+ };
+ vle_info vle_i;
+ vle_i.pvap = vec_packed(res, "pL / Pa");
+ vle_i.pliq = vec_packed(res, "pV / Pa");
+ vle_i.xL_0 = vec_packed(res, "xL_0 / mole frac.");
+ vle_i.xV_0 = vec_packed(res, "xV_0 / mole frac.");
+ vle.push_back(vle_i);
+ }
+ }
+
+ // Function for PVT mixture stuff
+ // function to calculate the density from temperature and pressure
+ template <typename ModelType>
+ void pvt_fill(const ModelType& model, psrk& helper) {
+ for (size_t i = 0; i < points.size(); i++) {
+ auto p_exp = points[i].P * 1E6;
+ tie(points[i].dense_v, points[i].dense_l) = get_rho_from_t_and_p(model, helper, points[i].T, p_exp, points[i].x);
+ }
+ }
+
+ template <typename ModelType, typename PureT>
+ void ve_fill(const ModelType& model, const PureT& pures, psrk& helper) {
+ for (size_t i = 0; i < points.size(); i++) {
+ auto p_exp = points[i].P * 1E6;
+ tie(points[i].dense_v, points[i].dense_l) = get_rho_from_t_and_p(model, helper, points[i].T, p_exp, points[i].x);
+ }
+ std::vector<double> x_pure = { 1.0 };
+ for (size_t i = 0; i < points.size(); i++) {
+ auto p_exp = points[i].P * 1E6;
+ auto [rho1_v, rho1_l] = get_rho_from_t_and_p(pures[0], helper, points[i].T, p_exp, x_pure);
+ auto [rho2_v, rho2_l] = get_rho_from_t_and_p(pures[1], helper, points[i].T, p_exp, x_pure);
+ points[i].ve_v = 1.0 / points[i].dense_v - points[i].x[0] / rho1_v - points[i].x[1] / rho2_v;
+ points[i].ve_l = 1.0 / points[i].dense_l - points[i].x[0] / rho1_l - points[i].x[1] / rho2_l;
+ }
+ }
+
+ template <typename ModelType, typename IdealT>
+ void snd_fill(const ModelType& model, const IdealT& ideal, psrk& helper) {
+ // First calculate density for all state points
+ for (size_t i = 0; i < points.size(); i++) {
+ auto p_exp = points[i].P * 1E6;
+ tie(points[i].dense_v, points[i].dense_l) = get_rho_from_t_and_p(model, helper, points[i].T, p_exp, points[i].x);
+ }
+ // get speed of sound for each phase
+ auto snd = [model, ideal](const double& T, const double& rho, std::vector<double>& molefrac, const auto& mw) {
+ using tdx = TDXDerivatives<ModelType, double, std::vector<double>>;
+ auto A01 = tdx::get_Ar01(model, T, rho, molefrac);
+ auto A02 = tdx::get_Ar02(model, T, rho, molefrac);
+ auto A11 = tdx::get_Ar11(model, T, rho, molefrac);
+ auto A20 = tdx::get_Ar20(model, T, rho, molefrac);
+
+ auto wih = AlphaCallWrapper<AlphaWrapperOption::idealgas, decltype(ideal)>(ideal);
+ using tdxi = TDXDerivatives<IdealT, double, std::vector<double> >;
+ auto A20_id = 0.0;
+ A20_id = tdxi::get_Agenxy<2, 0, ADBackends::autodiff>(wih, T, rho, molefrac);
+ return sqrt((1.0 + 2.0 * A01 + A02 - pow(1.0 + A01 - A11, 2.0) / (A20_id + A20)) * model.R(molefrac) / mw * T);
+ };
+ for (size_t i = 0; i < points.size(); i++) {
+ auto mw = 0.0;
+ mw = helper.get_mw(points[i].x);
+ points[i].snd_l = snd(points[i].T, points[i].dense_l, points[i].x, mw);
+ points[i].snd_v = snd(points[i].T, points[i].dense_v, points[i].x, mw);
+ }
+ }
+
+ // Selector for which data type should be calculated
+ template <typename ModelType, typename PureT, typename IdealT>
+ void perform_dev_calc(const ModelType& model, const PureT& pures, const IdealT& ideal, psrk& helper) {
+ std::string file_out = "FILE.PLT";
+ switch (prop)
+ {
+ case teqp::PVT:
+ pvt_fill(model, helper);
+ snd_fill(model, ideal, helper);
+ calc_dev_pvt();
+ write_pvt_dev(file_out);
+ break;
+ case teqp::VE:
+ ve_fill(model, pures, helper);
+ calc_dev_ve();
+ write_pvt_dev(file_out);
+ case teqp::SND:
+ snd_fill(model, ideal, helper);
+ calc_dev_snd();
+ write_pvt_dev(file_out);
+ case teqp::VLE:
+ calc_dev_vle();
+ break;
+ default:
+ break;
+ }
+ }
+
+
+
+ auto get_per_dev(double& one, double& two) {
+ return 100 * (one - two) / two;
+ }
+
+ // VLE DEV CALC
+ auto find_iso_idx(std::vector<double>& vec, double& key) {
+ auto itr = std::find(vec.begin(), vec.end(), key);
+ return std::distance(vec.begin(), itr);
+ }
+
+ auto find_closest(std::vector<double>& vec, double& value) {
+ auto itr = std::lower_bound(vec.begin(), vec.end(), value);
+ return std::distance(vec.begin(), itr);
+ }
+
+
+
+ void calc_dev_vle() {
+ for (size_t i = 0; i < points.size(); i++) {
+ auto idx = find_iso_idx(isotherms,points[i].T);
+ // get closes pvap to p
+ auto vle_p = vle[idx];
+ auto p_exp = points[i].P*1E6;
+ auto pvap_idx = find_closest(vle_p.pvap,p_exp);
+ auto x_vap_idx = vle_p.xV_0[pvap_idx];
+ }
+ }
+
+ // calculate deviation of pvt
+ void calc_dev_pvt() {
+ for (size_t i = 0; i < points.size(); i++) {
+ auto dexp = points[i].D * 1E3;
+ auto dev1 = get_per_dev(dexp, points[i].dense_l);
+ auto dev2 = get_per_dev(dexp, points[i].dense_v);
+ if (abs(dev1) < abs(dev2)) {
+ points[i].deviation = dev1;
+ points[i].D = points[i].dense_l;
+ }
+ else
+ {
+ points[i].deviation = dev2;
+ points[i].D = points[i].dense_v;
+ }
+ }
+ }
+
+ void calc_dev_snd() {
+ for (size_t i = 0; i < points.size(); i++) {
+ auto dev1 = get_per_dev(points[i].SND, points[i].snd_l);
+ auto dev2 = get_per_dev(points[i].SND, points[i].snd_v);
+ if (abs(dev1) < abs(dev2)) {
+ points[i].deviation = dev1;
+ }
+ else
+ {
+ points[i].deviation = dev2;
+ }
+ }
+ }
+
+ // calculate deviation of excess volume
+ void calc_dev_ve() {
+ for (size_t i = 0; i < points.size(); i++) {
+ auto dexp = points[i].D;
+ auto dev1 = get_per_dev(dexp, points[i].ve_l);
+ auto dev2 = get_per_dev(dexp, points[i].ve_v);
+ if (abs(dev1) < abs(dev2)) {
+ points[i].deviation = dev1;
+ }
+ else
+ {
+ points[i].deviation = dev2;
+ }
+ }
+ }
+
+ void write_pvt_dev(std::string& file_out) {
+ ofstream myfile;
+ myfile.open(file_out);
+ for (size_t i = 2; i < head.size(); i++) {
+ myfile << head[i] << std::endl;
+ }
+ for (size_t i = 0; i < points.size(); i++) {
+ std::string formatted_str = std::format("{0:>10.5f}{1:>10.5f}{2:>10.4f}{3:>10.4f}{4:>10.4f}{5:>8s}{6:>6d}{7:>12s}{8:>4d}", points[i].P, points[i].D, points[i].T, points[i].x[0], points[i].deviation, points[i].author, points[i].id, points[i].mix_name, points[i].nr);
+ myfile << formatted_str << std::endl;
+ }
+ myfile.close();
+ }
+
+
+ };
+
+ inline auto perform_plt() {
+
+ std::string set_path = "D:/OneDrivePrivat/OneDrive/Work/Diss/PHD/mixtures/compare/settings.json";
+ auto settings = get_compare_settings(set_path);
+ nlohmann::json j = load_a_JSON_file(set_path);
+ //auto m = make_model(j);
+ auto [pures, mixture, mix_ideal] = get_fluid_model(settings);
+ auto helper = psrk(settings);
+
+ std::vector< std::vector<data_point>> data;
+ for (size_t i = 0; i < settings.dev_files.size(); i++) {
+ data_set dat;
+ std::string dev_path = settings.dev_dir + "/" + settings.dev_files[i];
+ dat.read_dev_file(dev_path);
+ dat.get_binary_vle(mixture, pures);
+ dat.perform_dev_calc(mixture, pures, mix_ideal, helper);
+ }
+
+
+
+ return 0;
+ }
+
+
+}
+
+
+
+
+
diff --git a/include/teqp/algorithms/vle_trend.hpp b/include/teqp/algorithms/vle_trend.hpp
new file mode 100644
index 0000000..7c58f7f
--- /dev/null
+++ b/include/teqp/algorithms/vle_trend.hpp
@@ -0,0 +1,558 @@
+#pragma once
+
+#include "teqp/models/model_potentials/2center_ljf.hpp"
+#include "teqp/derivs.hpp"
+#include "teqp/models/multifluid_ancillaries.hpp"
+#include <Eigen/Dense>
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <teqp/models/cubics.hpp>
+#include <teqp/cpp/teqpcpp.hpp>
+#include "teqp/algorithms/VLE.hpp"
+#include "teqp/types.hpp"
+
+using namespace teqp;
+
+namespace teqp {
+
+ bool has_any_digits(const std::string& s) {
+ return std::any_of(s.begin(), s.end(), ::isdigit);
+ }
+
+ template <typename T>
+ typename std::enable_if<std::is_unsigned<T>::value, int>::type
+ inline constexpr signum(T const x) {
+ return T(0) < x;
+ }
+
+ template <typename T>
+ typename std::enable_if<std::is_signed<T>::value, int>::type
+ inline constexpr signum(T const x) {
+ return (T(0) < x) - (x < T(0));
+ }
+
+
+ enum phase { VAPOR, LIQUID };
+ enum prop_t { PVT, VLE , VE , SND};
+
+ struct compare_settings {
+ std::string model_type;
+ std::string model_path;
+ std::string srk_path;
+ std::string BIPcollection;
+ std::string depcollection;
+ std::string coolprop_root;
+ std::vector<std::string> mix;
+ std::vector<std::string> dev_files;
+ std::string dev_dir;
+ std::string combination;
+ };
+
+
+ struct fluid_info {
+ std::vector<double> tc, pc, dc, w;
+ };
+
+ struct phase_result
+ {
+
+ };
+
+ class psrk {
+
+ public:
+ std::vector<double> TC, PC, DC, w , mw;
+ psrk(compare_settings& s) {
+ std::string filepath = std::filesystem::is_regular_file(s.srk_path) ? s.srk_path : s.srk_path;
+ nlohmann::json j = load_a_JSON_file(filepath);
+ for (size_t i = 0; i < s.mix.size(); i++) {
+ TC.push_back(j.at("srk").at(s.mix[i]).at("tc"));
+ PC.push_back(j.at("srk").at(s.mix[i]).at("pc"));
+ DC.push_back(j.at("srk").at(s.mix[i]).at("dc"));
+ w.push_back(j.at("srk").at(s.mix[i]).at("omega"));
+ mw.push_back(j.at("srk").at(s.mix[i]).at("mw"));
+ }
+ }
+
+
+ auto get_mw(std::vector<double>& molefrac) {
+ auto m = 0.0;
+ for (size_t i = 0; i < molefrac.size(); i++)
+ {
+ m = m + molefrac[i] * mw[i];
+ }
+ return m;
+ }
+
+ std::vector<double> root_psrk(std::vector<double>& p) {
+
+ auto pie = 3.14159265358979;
+ std::vector<double> r(3);
+ auto a = p[0];
+ auto b = p[1];
+ auto c = p[2];
+ auto nr = 0;
+ auto nn = 0.0;
+ auto qq = 0.0;
+ auto pp = 0.0;
+ auto aa = 0.0;
+ auto bb = 0.0;
+
+ pp = (a * a - 3. * b) / 9.;
+ qq = (2. * pow(a, 3.0) - 9. * a * b + 27. * c) / 54.0;
+ if (pow(pp, 3.0) > pow(qq, 2.0)) {
+ nn = acos(-qq / pow(pp, 1.5));
+ nr = 3;
+ r[0] = 2. * sqrt(pp) * cos(nn / 3.) - a / 3.;
+ r[1] = 2. * sqrt(pp) * cos((nn + 2. * pie) / 3.) - a / 3.;
+ r[2] = 2. * sqrt(pp) * cos((nn + 4. * pie) / 3.) - a / 3.;
+ }
+ else {
+ nr = -1;
+ aa = -signum(qq) * pow((abs(qq) + sqrt(pow(qq, 2.0) - pow(pp, 3.0))), 1.0 / 3.0);
+ if (aa == 0) {
+ bb = 0.0;
+ }
+ else {
+ bb = pp / aa;
+ }
+
+ r[0] = (aa + bb) - a / 3.0;
+ r[1] = 0.0;
+ r[2] = 0.0;
+ }
+
+ // --creating output vector------------------------------------------------------------------ -
+ std::vector<double> volume(3);
+ if (nr < 0) {
+ volume[0] = -1.0;
+ volume[1] = r[0];
+ volume[2] = 0.0;
+ }
+
+
+ else {
+ volume[0] = 3.0;
+ auto idx = max_element(std::begin(r), std::end(r));
+ volume[1] = r[*idx];
+ auto qq = 1.0;
+ pp = 1.0;
+ nn = 1.0;
+ if (r[0] > 0) { qq = r[0]; }
+
+ if (r[1] > 0) { pp = r[1]; };
+ if (r[2] > 0) { nn = r[2]; };
+ std::vector<double> r_vec = { qq, pp, nn };
+ idx = min_element(std::begin(r_vec), std::end(r_vec));
+ volume[2] = *idx;
+ }
+
+ return volume;
+ }
+
+
+
+ auto get_density(const double& T, const double& P, const std::vector<double>& molefrac, const double& REQ, teqp::phase& phase_est) {
+ auto ncomp = molefrac.size();
+ std::vector<double> aa(ncomp), bb(ncomp), cc(3), msrk(ncomp), alfa_lc(ncomp);
+ std::vector<double> ai(ncomp), bi(ncomp), ci(ncomp);
+ auto a1 = 0.0;
+ auto c = 0.0;
+ auto b = 0.0;
+ for (auto i = 0; i < molefrac.size(); i++) {
+ auto tr = T / TC[i];
+ aa[i] = 0.42748 * REQ * REQ * TC[i] * TC[i] / PC[i];
+ msrk[i] = 0.48 + 1.574 * w[i] - 0.176 * pow(w[i], 2.0);
+ alfa_lc[i] = pow(1.0 + msrk[i] * (1.0 - sqrt(tr)), 2.0);
+ ai[i] = aa[i] * alfa_lc[i];
+ bi[i] = 0.08664 * REQ * TC[i] / PC[i];
+ ci[i] = 0.40768 * REQ * TC[i] / PC[i] * (0.29441 - PC[i] * 1.0 / DC[i] / REQ / TC[i]);
+ a1 = a1 + molefrac[i] * sqrt(ai[i]);
+ b = b + molefrac[i] * bi[i];
+ c = c + molefrac[i] * ci[i];
+ }
+
+ auto a = 0.0;
+ if (ncomp == 1) {
+ a = ai[0];
+ b = bi[0];
+ c = ci[0];
+ cc[0] = -REQ * T / P;
+ cc[1] = -b * b - REQ * T * b / P + a / P;
+ cc[2] = -a * b / P;
+ }
+ else {
+ a = pow(a1, 2.0);
+ cc[0] = -REQ * T / P;
+ cc[1] = -b * b - REQ * T * b / P + a / P; //rmix
+ cc[2] = -a * b / P;
+ }
+
+
+
+ std::vector<double> volume = root_psrk(cc);
+
+ auto psrk_val = 0.0;
+ auto VSRK = 0.0;
+ if (volume[0] < -10.0) {
+ psrk_val = -1.0;
+ }
+ else {
+ if (volume[0] < -0.0) {
+ VSRK = 1. / (volume[1]);
+ }
+ else {
+ if (phase_est == LIQUID) {
+ VSRK = 1. / (volume[2]);
+ }
+ else {
+ VSRK = 1. / (volume[1]);
+ }
+ }
+ }
+
+ auto vf1 = 0.0;
+ if (phase_est == LIQUID) {
+ vf1 = 1.0 / VSRK;
+ vf1 = vf1 - c;
+ VSRK = 1 / vf1;
+ }
+
+ psrk_val = VSRK;
+ return psrk_val;
+ }
+
+ };
+
+ template<typename Callable>
+ inline auto regular_falsi(Callable zerofunction, const double& Xstart_min, const double& Xstart_max, const double& Delta_Allowed, const double& Xmin_allowed, const double& Xmax_allowed, int& Max_iterations, int& Iterations) {
+
+ auto Xmin = Xstart_min;
+ auto Xmax = Xstart_max;
+ auto Xmin_alt = 0.0;
+ auto Xmax_alt = 0.0;
+ auto F_Xmin_alt = 0.0;
+ auto F_Xmax_alt = 0.0;
+ auto F_Xmin = 0.0;
+ auto F_Xmax = 0.0;
+ auto Stop_min = 0;
+ auto Stop_max = 0;
+ auto Int_Errorflag = 0;
+ bool root_in = false;
+ bool Root_Found = false;
+ bool additional_checks = false;
+ auto Count_enlarge = 0;
+ auto Count_reduce = 0;
+ auto result = 0.0;
+ auto Xakt = 0.0;
+ auto DERIV = 0.0;
+ auto F_Xakt = 0.0;
+ auto Deviation_X = 0.0;
+ auto Xroot = 0.0;
+ auto Xtest = 0.0;
+ auto F_Xtest = 0.0;
+ auto xmin = 0.0;
+ auto xmax = 0.0;
+ auto xakt = 0.0;
+
+ if (Xmin < Xmin_allowed) { Int_Errorflag = 1; }
+ if (Xmin < Xmin_allowed) { Xmin = Xmin_allowed; }
+ if (Xmax > Xmax_allowed) { Int_Errorflag = Int_Errorflag + 1; }
+ if (Xmax > Xmax_allowed) { Xmax = Xmax_allowed; }
+ if (Int_Errorflag == 1) { Int_Errorflag = 0; }
+ if (Int_Errorflag == 2) { Int_Errorflag = 1; }
+
+ if (Int_Errorflag == 0) {
+ F_Xmin = zerofunction(Xmin);
+ }
+
+ if (Int_Errorflag == 0) {
+ F_Xmax = zerofunction(Xmax);
+ }
+
+ if ((F_Xmin * F_Xmax) < 0) { root_in = true; }
+
+ // main loop
+ while (!(root_in) && Int_Errorflag == 0 && Count_enlarge <= 10) {
+ Xmin_alt = Xmin;
+ Xmax_alt = Xmax;
+ F_Xmin_alt = F_Xmin;
+ F_Xmax_alt = F_Xmax;
+ if (Stop_min == 0) { Xmin = Xmin - (Xstart_max - Xstart_min) * 0.5; }
+ if (Stop_max == 0) { Xmax = Xmax + (Xstart_max - Xstart_min) * 0.5; }
+ Count_enlarge = Count_enlarge + 1;
+ if ((Xmin - Xmin_allowed) < -1e-11) { Int_Errorflag = 1; }
+ if ((Xmax - Xmax_allowed) > 1e-11) { Int_Errorflag = 1; }
+ if ((Int_Errorflag == 0) && (Stop_min == 0)) { F_Xmin = zerofunction(Xmin); }
+ if ((Int_Errorflag == 0) && (Stop_min == 0)) { F_Xmax = zerofunction(Xmax); }
+ if ((abs(F_Xmin_alt) < abs(F_Xmin)) && ((F_Xmin_alt * F_Xmin) > 0)) { Stop_min = 1; }
+ if (Stop_min == 1) { Xmin = Xmin_alt; }
+ if (Stop_min == 1) { F_Xmin = F_Xmin_alt; }
+ if ((abs(F_Xmax_alt) < abs(F_Xmax)) && ((F_Xmax_alt * F_Xmax) > 0)) { Stop_max = 1; }
+ if (Stop_max == 1) { Xmax = Xmax_alt; }
+ if (Stop_max == 1) { F_Xmax = F_Xmax_alt; }
+ if ((F_Xmin * F_Xmax) < 0) { root_in = true; }
+ }
+
+ // If no root was found by enlarging the interval, there may be two roots within the initial interval.
+ // To check out this posibility, the original interval is subdivided starting in the midle of the original interval.
+ // The size of the new interval starts with 1 / 10 of the original interval and is enlarged then.
+
+ while (!(root_in) && (Int_Errorflag == 0) && (Count_reduce < 9)) {
+ Count_reduce = Count_reduce + 1;
+ Xmin = (Xstart_max + Xstart_min) / 2.0 - (Xstart_max - Xstart_min) / 20.0 * Count_reduce;
+ Xmax = (Xstart_max + Xstart_min) / 2.0 + (Xstart_max - Xstart_min) / 20.0 * Count_reduce;
+ //Check, whether Xmin and Xmax get out of the allowed range
+ if (Xmin < Xmin_allowed) { Int_Errorflag = 1; }
+ if (Xmax > Xmax_allowed) { Int_Errorflag = 1; }
+ if (Int_Errorflag == 0) { F_Xmin = zerofunction(Xmin); }
+ if (Int_Errorflag == 0) { F_Xmax = zerofunction(Xmax); }
+ //Check, whether a root is within the new limits
+ if ((F_Xmin * F_Xmax) < 0.0) { root_in = true; }
+ }
+
+ //If no root could be found in the original interval, by enlarging or by reducing of the interval the
+ //routine returns with Errorflag = 2
+
+ if (!(root_in)) { result = -2; }
+
+ if (result == 0) { // no error occured
+
+
+ while (!(Root_Found) && (Iterations <= Max_iterations)) {
+
+ Iterations = Iterations + 1; //Count iteration steps
+ auto Deriv = (F_Xmax - F_Xmin) / (Xmax - Xmin); //Calculate the "derivative" for the current interval;
+
+ if (abs(Deriv) > 1E-12) {
+ Deviation_X = F_Xakt / Deriv;//Deviation_X is the distance to the linear approximated root from the akt position on the x - axis
+ //SH 10 / 2018: If derivative is large make additional checks(see line 300) if root found or not
+ if (abs(Deriv) > 1E10) {
+ additional_checks = true;
+ }
+ else {
+ additional_checks = false;
+ }
+ }
+ else {
+ Deriv = signum(Deriv);
+ }
+
+ //Calculate Xakt starting from the X value with the smaller function value
+ if (abs(F_Xmin) < abs(F_Xmax)) {
+ // Calculate Xakt starting from the lower limit of the interval
+ if (abs(Deriv) > 0 && Iterations % 5 != 0) {
+ Xakt = Xmin + abs(F_Xmin / Deriv);
+ // Use this algorithm only three out of four times to avoid problems in regions // with strong curvature
+ }
+ else {
+ Xakt = (Xmin + Xmax) / 2.;
+ }
+ }
+ else {
+ if (abs(Deriv) > 0 && Iterations % 5 != 0) { // Calculate Xakt starting from the upper limit of the interval
+ Xakt = Xmax - abs(F_Xmax / Deriv);
+ } // Use this algorithm only three out of four times to avoid problems in regions
+ else { // with strong curvature
+ Xakt = (Xmin + Xmax) / 2.0;
+ }
+ }
+
+
+ F_Xakt = zerofunction(Xakt);
+
+
+ //Check whether the root has been found
+ //The criterion Delta_allowed is applied to the deviation in X
+
+ if (abs(Deriv) > 1E-12) {
+ Deviation_X = F_Xakt / Deriv;
+ if (abs(Deriv) > 1E10) {
+ additional_checks = true;
+ }
+ else {
+ additional_checks = false;
+ }
+ }
+ else {
+ Deriv = signum(Deriv);
+ }
+
+ if (((abs(Deviation_X) < Delta_Allowed) && (!(additional_checks))) || ((abs(Deviation_X) < Delta_Allowed) && (additional_checks) && ((abs((xakt - xmin) / xakt) < Delta_Allowed) || (abs((xakt - xmax) / xakt) < Delta_Allowed)))) {
+ Root_Found = true;
+ Xroot = Xakt + Deviation_X; // Xroot usually is a better guess for the root then Xakt
+ }
+ else //Root has not yet been found, new interval has to be defined
+ {
+ if (abs(Deviation_X) < ((Xmax - Xmin) / 10.0)) {
+
+ if (((Xakt - Xmin) < (Xmax - Xakt)) && ((Xakt - Xmin) > 0.0)) { Deriv = (F_Xakt - F_Xmin) / (Xakt - Xmin); }
+ if (((Xakt - Xmin) > (Xmax - Xakt)) && ((Xmax - Xakt) > 0.0)) { Deriv = (F_Xmax - F_Xakt) / (Xmax - Xakt); }
+
+ if (abs(Deriv) > 1E-12) {
+ Deviation_X = F_Xakt / Deriv;
+ }
+ else {
+ Deriv = signum(Deriv);
+ }
+
+ if ((F_Xmin * F_Xakt) < 0.0) {
+ Xtest = Xakt - 2.0 * abs(Deviation_X);
+ if (Xtest < ((Xmin + Xakt) / 2.0)) { Xtest = (Xmin + Xakt) / 2.0; }
+ F_Xtest = zerofunction(Xtest);
+ if ((F_Xtest * F_Xakt) < 0.0) { Xmin = Xtest; }
+ if ((F_Xtest * F_Xakt) < 0.0) { F_Xmin = F_Xtest; }
+ }
+ else if ((F_Xmax * F_Xakt) < 0.0) {
+ Xtest = Xakt + 2.0 * abs(Deviation_X);
+ if (Xtest > ((Xmax + Xakt) / 2.0)) { Xtest = (Xmax + Xakt) / 2.0; }
+ F_Xtest = zerofunction(Xtest);
+ if ((F_Xtest * F_Xakt) < 0.0) { Xmax = Xtest; }
+ if ((F_Xtest * F_Xakt) < 0.0) { F_Xmax = F_Xtest; }
+ }
+ }
+
+ if ((F_Xmin * F_Xakt) < 0.0) {
+ Xmax = Xakt;
+ F_Xmax = F_Xakt;
+ }
+ else if ((F_Xmax * F_Xakt) < 0.0) {
+ Xmin = Xakt;
+ F_Xmin = F_Xakt;
+ }
+ }
+
+ } // end while
+ } //end if
+
+ if ((Count_enlarge + Count_reduce) > 0) { Xroot = -1000; }
+ if (!(Root_Found)) { Xroot = -1000; }
+
+ return Xroot;
+
+ }
+
+
+ template <typename TypeT, typename TypeX, typename TypeK>
+ inline auto rache_rich(psrk& fldi, TypeK& K_val, TypeX& x_known, TypeT& vapfrac, int errval) {
+ auto ncomp = x_known.size();
+ auto sumx_bubble = 0.0;
+ auto sumx_dew = 0.0;
+ auto sumx_2phase = 0.0;
+ std::vector<double> x_vap(ncomp);
+ std::vector<double> x_liq(ncomp);
+ bool found = false;
+ auto frac_type = 0;
+ auto Delta_allowed = 1E-8;
+ auto frac_min_allowed = -1E-16;
+ auto frac_min = -1E-16;
+ auto frac_max = 0.5;
+ auto frac_max_allowed = 1.0 + 1E-16;
+ auto Max_iterations = 50;
+ auto frac = 0.0;
+
+ auto rache_diff = [frac_type, ncomp, K_val, x_known](double frac) {
+ auto rac_func = 0.0;
+ if (frac_type == 0) {
+ for (size_t i = 0; i < ncomp; i++) {
+ rac_func = rac_func + x_known[i] * (K_val[i] - 1.0) / (1.0 - frac + frac * K_val[i]);
+ }
+ }
+ else {
+ for (size_t i = 0; i < ncomp; i++)
+ {
+ rac_func = rac_func + x_known[i] * (K_val[i] - 1.0) / (frac + (1.0 - frac) * K_val[i]);
+ }
+ }
+ return rac_func;
+ };
+
+ for (size_t i = 0; i < ncomp; i++) {
+ sumx_bubble = sumx_bubble + x_known[i] * K_val[i];
+ sumx_dew = sumx_dew + x_known[i] / K_val[i];
+ sumx_2phase = sumx_2phase + x_known[i] * (K_val[i] - 1.0) / (K_val[i] + 1.0);
+ }
+
+ if (sumx_bubble < 1.0 + 1E-15) {
+ vapfrac = 0;
+ for (size_t i = 0; i < ncomp; i++) {
+ x_vap[i] = x_known[i] * K_val[i] / sumx_bubble;
+ }
+ x_liq = x_known;
+ found = true;
+ }
+
+ if (sumx_dew < 1.0 + 1E-15) {
+ vapfrac = 0;
+ for (size_t i = 0; i < ncomp; i++) {
+ x_liq[i] = x_known[i] / K_val[i] / sumx_dew;
+ }
+ x_vap = x_known;
+ found = true;
+ }
+
+ //The mixture is in the 2 phase region with 0.5 < vapfrac < 1.
+ //Instead of vapfrac the liquid fraction alpha is used, since vapfrac near unity might cause round off errors(see Michelsen& Mollerup)
+ if (!(found)) {
+ //The liquid fraction is chosen as independent variable
+ if (sumx_2phase > 0.0) {
+ frac_type = 1;
+ }
+ else {
+ frac_type = 0;
+ }
+
+ auto Iterations = 0;
+ frac = 0.0;
+ auto res = regular_falsi(rache_diff, frac_min, frac_max, Delta_allowed, frac_min_allowed, frac_max_allowed, Max_iterations, Iterations);
+ if (errval == 3) { errval = 0;};
+ if (errval != 0) { errval = -2111; }
+ if (frac_type == 1) {
+ vapfrac = 1.0 - frac;
+ }
+ else {
+ vapfrac = frac;
+ }
+
+ sumx_dew = 0.0;
+ sumx_bubble = 0.0;
+ for (size_t i = 0; i < x_known.size(); i++) {
+ x_vap[i] = x_known[i] * K_val[i] / (1.0 - vapfrac + vapfrac * K_val[i]);
+ sumx_dew = sumx_dew + x_vap[i];
+ x_liq[i] = x_known[i] / (1.0 - vapfrac + vapfrac * K_val[i]);
+ sumx_bubble = sumx_bubble + x_liq[i];
+ }
+
+ for (size_t i = 0; i < x_vap.size(); i++)
+ {
+ x_vap[i] = x_vap[i] / sumx_dew;
+ x_liq[i] = x_liq[i] / sumx_bubble;
+ }
+
+
+ }
+ return std::make_tuple(x_vap, x_liq, vapfrac);
+ }
+ template <typename TypeT, typename TypeP, typename TypeX>
+ inline auto PTX_startvals_pT(psrk& fldi, TypeT& T, TypeP p, TypeX& x_known, TypeT& vapfrac, int errval) {
+
+ std::vector<double> K_val(x_known.size());
+ for (size_t i = 0; i < x_known.size(); i++) {
+ K_val[i] = fldi.PC[i] / p * exp(5.3730 * (1.0 + fldi.w[i]) * (1.0 - fldi.TC[i] / T));
+ }
+
+ return rache_rich(fldi, K_val, x_known, vapfrac, errval);
+
+ }
+ template <typename TypeT, typename TypeP>
+ inline auto PhaseDet(psrk& fldi, TypeT& T, TypeP& p, const std::vector<double>& x) {
+
+ auto errval = 0;
+ auto vapfrac = 0.0;
+ std::vector<TypeT> x_vap, x_liq2, x_liq;
+ tie(x_vap, x_liq, vapfrac) = PTX_startvals_pT(fldi, T, p, x, vapfrac, errval);
+ return 0.0;
+ }
+
+}
\ No newline at end of file
diff --git a/include/teqp/models/mie/mie.hpp b/include/teqp/models/mie/mie.hpp
index 5d7985b..c80a82b 100644
--- a/include/teqp/models/mie/mie.hpp
+++ b/include/teqp/models/mie/mie.hpp
@@ -2,71 +2,1109 @@
#include "teqp/models/multifluid.hpp"
#include <Eigen/Dense>
+#include <math.h>
+#include "teqp/models/cubics.hpp"
namespace teqp {
-namespace Mie{
-
- /**
- Equation of state for the Mie ($\lambda_{\mathrm{r}}$,6) fluid with a repulsive exponent from 11 to 13
- Pohl, S.; Fingerhut, R., Thol, M.; Vrabec, J.; Span, R.
- */
- class Mie6Pohl2023 {
-
- private:
- using EArray6 = Eigen::Array<double, 6, 1>;
- using EArray4 = Eigen::Array<double, 4, 1>;
-
- const EArray6 c1_pol = (EArray6() << -0.0192944,1.38,-2.2653,1.6291,-1.974,0.40412 ).finished();
- const EArray6 c1_exp = (EArray6() << 0.1845,-0.3227,1.1351,2.232,-2.344,-0.4238 ).finished();
- const EArray4 c1_gbs = (EArray4() << -4.367,0.0371,1.3895,2.835 ).finished();
- const EArray6 c2_pol = (EArray6() << 0.26021,-5.525,8.329,-19.492,25.8,-3.8133 ).finished();
- const EArray6 c2_exp = (EArray6() << -5.05,2.7842,-9.523,-30.383,17.902,2.2264 ).finished();
- const EArray4 c2_gbs = (EArray4() << 48.445,-5.506,-11.643,-24.36 ).finished();
- const EArray6 t_pol = (EArray6() << 1,0.236,0.872,0.313,0.407,0.703 ).finished();
- const EArray6 t_exp = (EArray6() << 1.78,2.99,2.866,1.2,3.06,1.073 ).finished();
- const EArray4 t_gbs = (EArray4() << 1.50,1.03,4.02,1.57 ).finished();
- const EArray6 d_pol = (EArray6() << 4,1,1,2,2,3 ).finished();
- const EArray6 d_exp = (EArray6() << 1,1,3,2,2,5 ).finished();
- const EArray4 d_gbs = (EArray4() << 2,3,2,2 ).finished();
- const EArray6 p = (EArray6() << 1,2,2,1,2,1 ).finished();
- const EArray4 eta = (EArray4() << 0.362,0.313,1.17,0.957 ).finished();
- const EArray4 beta = (EArray4() << 0.0761,0.143,0.63,1.32 ).finished();
- const EArray4 gam = (EArray4() << 1.55,-0.0826,1.505,1.07 ).finished();
- const EArray4 eps = (EArray4() << -1,-1,-0.195,-0.287 ).finished();
-
- const double m_lambda_a;
- const EArray6 n_pol, n_exp;
- const EArray4 n_gbs;
- const double Tc, rhoc; // In simulation units
- public:
-
- Mie6Pohl2023(double lambda_a) : m_lambda_a(lambda_a),
- n_pol(c1_pol + c2_pol / m_lambda_a),
- n_exp(c1_exp + c2_exp / m_lambda_a),
- n_gbs(c1_gbs + c2_gbs / m_lambda_a),
- Tc(0.668 + 6.84 / m_lambda_a + 145 / pow(m_lambda_a, 3)), // T^*
- rhoc(0.2516 + 0.049 * log10(m_lambda_a)) // rho^*
- {}
-
- auto get_lambda_a() const { return m_lambda_a; }
-
- // We are in "simulation units", so R is 1.0, and T and rho that
- // go into alphar are actually T^* and rho^*
- template<typename MoleFracType>
- double R(const MoleFracType &) const { return 1.0; }
-
- template<typename TTYPE, typename RHOTYPE, typename MoleFracType>
- auto alphar(const TTYPE& Tstar, const RHOTYPE& rhostar, const MoleFracType& /*molefrac*/) const {
- auto tau = Tc / Tstar; auto delta = rhostar / rhoc;
- auto alphar_ = (
- (n_pol * pow(tau, t_pol) * pow(delta, d_pol)).sum() +
- (n_exp * pow(tau, t_exp) * pow(delta, d_exp) * exp(-pow(delta, p))).sum() +
- (n_gbs * pow(tau, t_gbs) * pow(delta, d_gbs) * exp(-eta * pow(delta - eps, 2) - beta * pow(tau - gam, 2))).sum()
- );
- return forceeval(alphar_);
+ namespace Mie {
+
+ /**
+ Equation of state for the Mie ($\lambda_{\mathrm{r}}$,6) fluid with a repulsive exponent from 11 to 13
+ Pohl, S.; Fingerhut, R., Thol, M.; Vrabec, J.; Span, R.
+ */
+ class Mie6Pohl2023 {
+
+ private:
+ using EArray6 = Eigen::Array<double, 6, 1>;
+ using EArray4 = Eigen::Array<double, 4, 1>;
+
+ const EArray6 c1_pol = (EArray6() << -0.0192944, 1.38, -2.2653, 1.6291, -1.974, 0.40412).finished();
+ const EArray6 c1_exp = (EArray6() << 0.1845, -0.3227, 1.1351, 2.232, -2.344, -0.4238).finished();
+ const EArray4 c1_gbs = (EArray4() << -4.367, 0.0371, 1.3895, 2.835).finished();
+ const EArray6 c2_pol = (EArray6() << 0.26021, -5.525, 8.329, -19.492, 25.8, -3.8133).finished();
+ const EArray6 c2_exp = (EArray6() << -5.05, 2.7842, -9.523, -30.383, 17.902, 2.2264).finished();
+ const EArray4 c2_gbs = (EArray4() << 48.445, -5.506, -11.643, -24.36).finished();
+ const EArray6 t_pol = (EArray6() << 1, 0.236, 0.872, 0.313, 0.407, 0.703).finished();
+ const EArray6 t_exp = (EArray6() << 1.78, 2.99, 2.866, 1.2, 3.06, 1.073).finished();
+ const EArray4 t_gbs = (EArray4() << 1.50, 1.03, 4.02, 1.57).finished();
+ const EArray6 d_pol = (EArray6() << 4, 1, 1, 2, 2, 3).finished();
+ const EArray6 d_exp = (EArray6() << 1, 1, 3, 2, 2, 5).finished();
+ const EArray4 d_gbs = (EArray4() << 2, 3, 2, 2).finished();
+ const EArray6 p = (EArray6() << 1, 2, 2, 1, 2, 1).finished();
+ const EArray4 eta = (EArray4() << 0.362, 0.313, 1.17, 0.957).finished();
+ const EArray4 beta = (EArray4() << 0.0761, 0.143, 0.63, 1.32).finished();
+ const EArray4 gam = (EArray4() << 1.55, -0.0826, 1.505, 1.07).finished();
+ const EArray4 eps = (EArray4() << -1, -1, -0.195, -0.287).finished();
+
+ const double m_lambda_a;
+ const EArray6 n_pol, n_exp;
+ const EArray4 n_gbs;
+ const double Tc, rhoc; // In simulation units
+ public:
+
+ Mie6Pohl2023(double lambda_a) : m_lambda_a(lambda_a),
+ n_pol(c1_pol + c2_pol / m_lambda_a),
+ n_exp(c1_exp + c2_exp / m_lambda_a),
+ n_gbs(c1_gbs + c2_gbs / m_lambda_a),
+ Tc(0.668 + 6.84 / m_lambda_a + 145 / pow(m_lambda_a, 3)), // T^*
+ rhoc(0.2516 + 0.049 * log10(m_lambda_a)) // rho^*
+ {}
+
+ auto get_lambda_a() const { return m_lambda_a; }
+
+ // We are in "simulation units", so R is 1.0, and T and rho that
+ // go into alphar are actually T^* and rho^*
+ template<typename MoleFracType>
+ double R(const MoleFracType&) const { return 1.0; }
+
+ template<typename TTYPE, typename RHOTYPE, typename MoleFracType>
+ auto alphar(const TTYPE& Tstar, const RHOTYPE& rhostar, const MoleFracType& /*molefrac*/) const {
+ auto tau = Tc / Tstar; auto delta = rhostar / rhoc;
+ auto alphar_ = (
+ (n_pol * pow(tau, t_pol) * pow(delta, d_pol)).sum() +
+ (n_exp * pow(tau, t_exp) * pow(delta, d_exp) * exp(-pow(delta, p))).sum() +
+ (n_gbs * pow(tau, t_gbs) * pow(delta, d_gbs) * exp(-eta * pow(delta - eps, 2) - beta * pow(tau - gam, 2))).sum()
+ );
+ return forceeval(alphar_);
+ }
+
+ };
+
+
+ //
+ inline auto linear_mixing(const double& x, const double& y) {
+ return (x + y) / 2.0;
+ }
+
+ class AncTerm {
+ public:
+ Eigen::ArrayXd n, t;
+ double tred, dred;
+ std::string type;
+ };
+
+ class fluid {
+ public:
+ double lambdas, epsilons, sigmas, m, I;
+ AncTerm anc_dl, anc_dv;
+ bool is_sphere;
+ };
+
+ // Combining rules for one fluid approximation
+ // (1) Simple Van der Waals one fluid combininb rule
+ template<typename RHOTYPE, typename MoleFracType>
+ inline auto combining_rules_one_fluid(RHOTYPE& rhostar, MoleFracType& molefrac, std::vector<fluid> f, const bool& is_spherical) {
+ auto ncomp = f.size();
+ using resulttype = std::common_type_t<decltype(molefrac[0]), decltype(rhostar)>;
+ using rhotype = std::common_type_t<decltype(rhostar)>;
+ std::vector<std::vector<resulttype>> sigma_ij(ncomp, std::vector<resulttype>(ncomp, 0.0));
+ std::vector<std::vector<resulttype>> eps_ij(ncomp, std::vector<resulttype>(ncomp, 0.0));
+ std::vector<std::vector<resulttype>> lambda_ij(ncomp, std::vector<resulttype>(ncomp, 0.0));
+ std::vector<resulttype> m_mixed(ncomp, 0.0);
+ resulttype m_mix = 0.0;
+ resulttype sigma_mean = 0.0;
+ resulttype epsilon_mean = 0.0;
+ resulttype lambda_mean = 0.0;
+
+ if (is_spherical) {
+ m_mix = 1.0;
+ }
+ else
+ {
+ for (auto i = 0; i < ncomp; i++) {
+ m_mix += f[i].m * molefrac[i];
+ }
+ }
+
+ resulttype kij = 0.00;
+ for (size_t i = 0; i < ncomp - 1; i++)
+ {
+ //if(molefrac[0]<1.0) {
+ kij = 0.0;//1.0 - pow(2.0, 7.0) * (sqrt(f[i].I * f[i + 1].I) / (f[i].I + f[i + 1].I)) * (pow(f[i].sigmas, 3.0) * pow(f[i + 1].sigmas, 3.0)) / pow(f[i].sigmas + f[i + 1].sigmas, 6);
+ /*}
+ else
+ {
+ kij = 0.0;
+ }*/
+
+ }
+ //
+
+ // Switch between combing rules for interaction of molecular parameters
+ for (auto i = 0; i < ncomp; i++) {
+ for (auto j = 0; j < ncomp; j++) {
+ lambda_ij[i][j] = linear_mixing(f[i].lambdas, f[j].lambdas); //3.0 + sqrt((f[i].lambdas - 3.0) * (f[j].lambdas - 3.0));
+ sigma_ij[i][j] = linear_mixing(f[i].sigmas, f[j].sigmas);
+ eps_ij[i][j] = (1.0 - kij) * sqrt(f[i].epsilons * f[j].epsilons);//*pow(sqrt(f[i].sigmas*f[j].sigmas)/ sigma_ij[i][j],6.0); //sqrt(pow(f[i].sigmas,3.0)*pow(f[j].sigmas,3.0))/pow(sigma_ij[i][j],3.0)
+ }
+ }
+
+ // Calculate the mean values for the one fluid approximation
+ for (auto i = 0; i < ncomp; i++) { for (auto j = 0; j < ncomp; j++) { sigma_mean += molefrac[i] * molefrac[j] * f[i].m * f[j].m * pow(sigma_ij[i][j], 3.0); } }
+ for (auto i = 0; i < ncomp; i++) { for (auto j = 0; j < ncomp; j++) { epsilon_mean += molefrac[i] * molefrac[j] * f[i].m * f[j].m * pow(sigma_ij[i][j], 3.0) * eps_ij[i][j]; } }
+ for (auto i = 0; i < ncomp; i++) { for (auto j = 0; j < ncomp; j++) { lambda_mean += molefrac[i] * molefrac[j] * lambda_ij[i][j]; } }
+
+
+ epsilon_mean = epsilon_mean / sigma_mean;
+ sigma_mean = pow(sigma_mean / pow(m_mix, 2.0), 1.0 / 3.0);
+ return std::make_tuple(sigma_mean, epsilon_mean, lambda_mean, m_mix);
}
- };
-};
-};
+ inline auto dlv_eq1(const double& T, const AncTerm& anc) {
+ auto tau = abs(1.0 - T / anc.tred);
+ auto sum = 0.0;
+ for (size_t i = 0; i < anc.n.size(); i++) {
+ sum = sum + anc.n(i) * pow(tau, anc.t(i));
+ }
+ return anc.dred * (sum + 1.0);
+ }
+
+ inline auto dlv_eq2(const double& T, const AncTerm& anc) {
+ auto tau = pow(abs(1.0 - T / anc.tred), 1.0 / 3.0);
+ auto sum = 0.0;
+ for (size_t i = 0; i < anc.n.size(); i++) {
+ sum = sum + anc.n(i) * pow(tau, anc.t(i));
+ }
+ return anc.dred * (sum + 1.0);
+ }
+
+ inline auto dlv_eq3(const double& T, const AncTerm& anc) {
+ auto tau = abs(1.0 - T / anc.tred);
+ auto sum = 0.0;
+ for (size_t i = 0; i < anc.n.size(); i++) {
+ sum = sum + anc.n(i) * pow(tau, anc.t(i));
+ }
+ return anc.dred * exp(sum);
+ }
+
+ inline auto dlv_eq4(const double& T, const AncTerm& anc) {
+ auto tau = pow(abs(1.0 - T / anc.tred), 1.0 / 3.0);
+ auto sum = 0.0;
+ for (size_t i = 0; i < anc.n.size(); i++) {
+ sum = sum + anc.n(i) * pow(tau, anc.t(i));
+ }
+ return anc.dred * exp(sum);
+ }
+
+ inline auto dlv_eq5(const double& T, const AncTerm& anc) {
+ auto tau = abs(1.0 - T / anc.tred);
+ auto sum = 0.0;
+ for (size_t i = 0; i < anc.n.size(); i++) {
+ sum = sum + anc.n(i) * pow(tau, anc.t(i));
+ }
+ return anc.dred * exp(anc.tred / T * sum);
+ }
+
+ inline auto dlv_eq6(const double& T, const AncTerm& anc) {
+ auto tau = pow(abs(1.0 - T / anc.tred), 1.0 / 3.0);
+ auto sum = 0.0;
+ for (size_t i = 0; i < anc.n.size(); i++) {
+ sum = sum + anc.n(i) * pow(tau, anc.t(i));
+ }
+ return anc.dred * exp(anc.tred / T * sum);
+ }
+
+ // Ancillaries for pure fluid components to generate starting of vle
+ std::tuple<double, double> get_sat_dense(double& T, fluid f) {
+ double dl = 0.0;
+ double dv = 0.0;
+ if (f.anc_dl.type == "DL1") { dl = dlv_eq1(T, f.anc_dl); }
+ else if (f.anc_dl.type == "DL2") { dl = dlv_eq2(T, f.anc_dl); }
+ else if (f.anc_dl.type == "DL3") { dl = dlv_eq3(T, f.anc_dl); }
+ else if (f.anc_dl.type == "DL4") { dl = dlv_eq4(T, f.anc_dl); }
+ else if (f.anc_dl.type == "DL5") { dl = dlv_eq5(T, f.anc_dl); }
+ else if (f.anc_dl.type == "DL6") { dl = dlv_eq6(T, f.anc_dl); }
+
+ if (f.anc_dv.type == "DV1") { dv = dlv_eq1(T, f.anc_dv); }
+ else if (f.anc_dv.type == "DV2") { dv = dlv_eq2(T, f.anc_dv); }
+ else if (f.anc_dv.type == "DV3") { dv = dlv_eq3(T, f.anc_dv); }
+ else if (f.anc_dv.type == "DV4") { dv = dlv_eq4(T, f.anc_dv); }
+ else if (f.anc_dv.type == "DV5") { dv = dlv_eq5(T, f.anc_dv); }
+ else if (f.anc_dv.type == "DV6") { dv = dlv_eq6(T, f.anc_dv); }
+
+ return { dv,dl };
+ }
+
+ template<typename Model>
+ auto get_vle_pure_start(const Model& model, double& T, const int idx) {
+ return get_sat_dense(T, model.fld[idx]);
+ }
+
+ const double kBoltz = 1.380649E-23;
+ const double NAvo = 6.02214076E23;
+
+
+ class MieElong {
+
+ private:
+
+ const double kBoltz = 1.380649E-23;
+ const double NAvo = 6.02214076E23;
+
+ using EArray6 = Eigen::Array<double, 6, 1>;
+ using EArray4 = Eigen::Array<double, 4, 1>;
+
+
+ std::valarray<double> c1_pol = { -0.0192944, 1.38, -2.2653, 1.6291, -1.974, 0.40412 };
+ std::valarray<double> c1_exp = { 0.1845, -0.3227, 1.1351, 2.232, -2.344, -0.4238 };
+ std::valarray<double> c1_gbs = { -4.367, 0.0371, 1.3895, 2.835 };
+ std::valarray<double> c2_pol = { 0.26021, -5.525, 8.329, -19.492, 25.8, -3.8133 };
+ std::valarray<double> c2_exp = { -5.05, 2.7842, -9.523, -30.383, 17.902, 2.2264 };
+ std::valarray<double> c2_gbs = { 48.445, -5.506, -11.643, -24.36 };
+ std::valarray<double> t_pol = { 1.0, 0.236, 0.872, 0.313, 0.407, 0.703 };
+ std::valarray<double> t_exp = { 1.78, 2.99, 2.866, 1.2, 3.06, 1.073 };
+ std::valarray<double> t_gbs = { 1.50, 1.03, 4.02, 1.57 };
+ std::valarray<double> d_pol = { 4.0, 1.0, 1.0, 2.0, 2.0, 3.0 };
+ std::valarray<double> d_exp = { 1.0, 1.0, 3.0, 2.0, 2.0, 5.0 };
+ std::valarray<double> d_gbs = { 2.0, 3.0, 2.0, 2.0 };
+ std::valarray<double> p = { 1.0, 2.0, 2.0, 1.0, 2.0, 1.0 };
+ std::valarray<double> eta = { 0.362, 0.313, 1.17, 0.957 };
+ std::valarray<double> beta = { 0.0761, 0.143, 0.63, 1.32 };
+ std::valarray<double> gam = { 1.55, -0.0826, 1.505, 1.07 };
+ std::valarray<double> eps = { -1.0, -1.0, -0.195, -0.287 };
+ //std::valarray<double> d1_pol = { -0.0050505, 0.38601, -1.1222, -0.0039209, 0.050028, -0.0013054 };
+ //std::valarray<double> d2_pol = { 0.016623, 0.53098, -1.0766, 0.26569, 0.013785, 0.0013796 };
+ //std::valarray<double> d3_pol = { -0.00042696, 0.18677, 0.094444, 0.01354, 0.0051845, -0.0011878 };
+ //std::valarray<double> d1_exp = { 0.12441, -0.094619, -0.16034, 0.074307, -0.011307, 0.0080606 };
+ //std::valarray<double> d2_exp = { -0.24797, 0.42032, 0.46332, -0.051815, 0.039351, 0.0017919 };
+ //std::valarray<double> d3_exp = { 0.17459, -0.53721, -0.12545, -1.1854, 0.98169, 0.058886 };
+ //std::valarray<double> t_pol_m = { 1, 0.367124560272789, 1.2185476818727, 0.781480878568594, 1.64331019262472, 0.587674327418418 };
+ //std::valarray<double> t_exp_m = { 4.47186031183114, 2.82497155574424, 3.77610727931397, 2.65845878597852, 3.80375972303267, 0.665886883638958 };
+ //std::valarray<double> d_pol_m = { 4.0, 1.0, 1.0, 2.0, 2.0, 3.0 };
+ //std::valarray<double> d_exp_m = { 1.0, 1.0, 3.0, 2.0, 2.0, 5.0 };
+ //std::valarray<double> p_m = { 1.0 , 2.0 , 2.0 , 1.0 , 2.0 , 1.0 };
+ //std::valarray<double> d1_pol = { -0.0050505, 0.38601, -1.1222, -0.0039209, 0.050028, -0.0013054 };
+ //std::valarray<double> d2_pol = { 0.016623, 0.53098, -1.0766, 0.26569, 0.013785, 0.0013796 };
+ //std::valarray<double> d3_pol = { -0.00042696, 0.18677, 0.094444, 0.01354, 0.0051845, -0.0011878 };
+ //std::valarray<double> d1_exp = { 0.12441, -0.094619, -0.16034, 0.074307, -0.011307, 0.0080606 };
+ //std::valarray<double> d2_exp = { -0.24797, 0.42032, 0.46332, -0.051815, 0.039351, 0.0017919 };
+ //std::valarray<double> d3_exp = { 0.17459, -0.53721, -0.12545, -1.1854, 0.98169, 0.058886 };
+ //std::valarray<double> t_pol_m = { 1, 0.367124560272789, 1.2185476818727, 0.781480878568594, 1.64331019262472, 0.587674327418418 };
+ //std::valarray<double> t_exp_m = { 4.47186031183114, 2.82497155574424, 3.77610727931397, 2.65845878597852, 3.80375972303267, 0.665886883638958 };
+ //std::valarray<double> d_pol_m = { 4.0, 1.0, 1.0, 2.0, 2.0, 3.0 };
+ //std::valarray<double> d_exp_m = { 1.0, 1.0, 3.0, 2.0, 2.0, 5.0 };
+ //std::valarray<double> p_m = { 1.0 , 2.0 , 2.0 , 1.0 , 2.0 , 1.0 };
+ std::valarray<double> d1_pol = { -0.044326, -0.9164, 0.54604, 0.93411, -0.36957, -0.20822 };
+ std::valarray<double> d2_pol = { 0.0010372, -0.48659, 0.36519, 2.2719, -2.1598, -0.071921 };
+ std::valarray<double> d3_pol = { 0.045309, -0.080824, 0.0017117, 1.4648, -1.6084, 0.052748 };
+ std::valarray<double> d1_exp = { 0.51394, 0.55711, 0.083787, -0.86286, 0.19885, 0.02952 };
+ std::valarray<double> d2_exp = { 0.24318, 0.68819, 0.36777, -1.0891, 0.17631, -0.0021117 };
+ std::valarray<double> d3_exp = { 0.056978, 0.60157, 0.38942, -0.70837, -0.043016, -0.027741 };
+ std::valarray<double> d1_gbs = { 0.79761, -0.88304, 0.96654, -0.37683 };
+ std::valarray<double> d2_gbs = { 0.60753, -0.9869, 1.5873, -0.50984 };
+ std::valarray<double> d3_gbs = { 0.22228, -0.34791, 0.88312, -0.35055 };
+ std::valarray<double> t_pol_m = { 1, 0.216271942617027, 0.520609388112777, 0.752133081448595, 0.775679663317215, 0.580017747329082 };
+ std::valarray<double> t_exp_m = { 1.96734718455893, 2, 1.38627146150207, 1.4813773563886, 1.2103073201027, 1.01007005697967 };
+ std::valarray<double> t_gbs_m = { 0.418410732352989, 0.691563451339712, 1.31680108926219, 1.99860010157709 };
+ std::valarray<double> d_pol_m = { 4.0, 1.0, 1.0, 2.0, 2.0, 3.0 };
+ std::valarray<double> d_exp_m = { 1.0, 3.0, 2.0, 2.0, 4.0, 7.0 };
+ std::valarray<double> d_gbs_m = { 1,2,2,3 };
+ std::valarray<double> p_m = { 2.0 , 2.0 , 1.0 , 2.0 , 1.0 , 1.0 };
+ std::valarray<double> eta_m = { 0.240229835057927, 0.493986587611359, 1.38196451506317 , 1.18724195184123 };
+ std::valarray<double> beta_m = { 0.814951788107642, 1.72306756221527 ,0.452352097018578 ,0.0656824491597725 };
+ std::valarray<double> gam_m = { -0.403458611209436,0.0471338211178168 ,-0.321710057373199 , 0.494527960329717 };
+ std::valarray<double> eps_m = { -0.482584264913269,-0.435554887414116 ,0.271841899896465 ,0.198724624217954 };
+
+ int nr_fld;
+ bool is_sphere;
+
+ public:
+
+ std::vector<fluid> fld;
+ std::string combining_rule;
+
+ //MieElong(double lambda_a, double ms, double sigma, double eps) { // Repulsive exponent, segment number ms
+ MieElong(const std::string& path, const std::vector<std::string>& fluids, const std::string& combining_rule_in) {
+ std::string filepath = std::filesystem::is_regular_file(path) ? path : path;
+ nlohmann::json j = load_a_JSON_file(filepath);
+
+ fld.resize(fluids.size());
+ int i = 0;
+ for (auto f : fluids)
+ {
+ fld[i].lambdas = j.at("mie").at(f).at("lambda");
+ fld[i].epsilons = j.at("mie").at(f).at("epsilon");
+ fld[i].sigmas = j.at("mie").at(f).at("sigma");
+ fld[i].m = j.at("mie").at(f).at("segment");
+ fld[i].I = j.at("mie").at(f).at("I");
+
+ auto n_anc = static_cast<int>(j.at("mie").at(f).at("dl_coeff").size());
+ fld[i].anc_dl.tred = j.at("mie").at(f).at("tred");
+ fld[i].anc_dl.dred = j.at("mie").at(f).at("dred");
+ fld[i].anc_dl.type = j.at("mie").at(f).at("dl_type");
+ fld[i].anc_dl.n = toeig(j.at("mie").at(f).at("dl_coeff")).head(n_anc);
+ fld[i].anc_dl.t = toeig(j.at("mie").at(f).at("dl_texp")).head(n_anc);
+
+ n_anc = static_cast<int>(j.at("mie").at(f).at("dv_coeff").size());
+ fld[i].anc_dv.tred = j.at("mie").at(f).at("tred");
+ fld[i].anc_dv.dred = j.at("mie").at(f).at("dred");
+ fld[i].anc_dv.type = j.at("mie").at(f).at("dv_type");
+ fld[i].anc_dv.n = toeig(j.at("mie").at(f).at("dv_coeff")).head(n_anc);
+ fld[i].anc_dv.t = toeig(j.at("mie").at(f).at("dv_texp")).head(n_anc);
+
+ i++;
+ }
+
+ nr_fld = fluids.size();
+ is_sphere = true;
+ combining_rule = combining_rule_in;
+ // check if only sphericals are included
+ for (size_t i = 0; i < nr_fld; i++)
+ {
+ if (fld[i].m > 1.0) {
+ is_sphere = false;
+ }
+ }
+ }
+
+ template<typename MoleFracType>
+ double R(const MoleFracType&) const { return NAvo * kBoltz; }
+
+ template<typename TTYPE, typename MoleFracType, typename GammaType>
+ inline auto get_tred_mix(const TTYPE& tc, const MoleFracType& x, const GammaType& gama_t) const {
+ typename TTYPE::value_type tred_mix = 0.0;
+ for (size_t i = 0; i < x.size(); i++) {
+ tred_mix += x[i] * x[i] * tc[i];
+ }
+ auto beta_t = 1.0;
+ for (size_t i = 0; i < x.size() - 1; i++) {
+ for (size_t j = i + 1; j < x.size(); j++) {
+ tred_mix += 2.0 * x[i] * x[j] * beta_t * gama_t[i][j] * (x[i] + x[j]) / (beta_t * beta_t * x[i] + x[j]) * sqrt(tc[i] * tc[j]);
+ }
+ }
+ return tred_mix;
+ }
+
+ template<typename TTYPE, typename MoleFracType, typename GammaType>
+ inline auto get_vred_mix(const TTYPE& dc, const MoleFracType& x, const GammaType& gama_v) const {
+ typename TTYPE::value_type vred_mix1 = 0.0;
+ for (size_t i = 0; i < x.size(); i++) {
+ vred_mix1 += (x[i] * x[i]) / dc[i];
+ }
+ auto beta_v = 1.0;
+ for (size_t i = 0; i < x.size() - 1; i++) {
+ for (size_t j = i + 1; j < x.size(); j++) {
+ vred_mix1 += 2.0 * x[i] * x[j] * beta_v * gama_v[i][j] * (x[i] + x[j]) / (beta_v * beta_v * x[i] + x[j]) * 0.125 * pow(1.0 / pow(dc[i], 1.0 / 3.0) + 1.0 / pow(dc[j], 1.0 / 3.0), 3.0);
+ }
+ }
+
+ return 1.0 / vred_mix1;
+ }
+
+ template<typename ETYPE, typename LTYPE, typename MTYPE>
+ inline auto get_tc(const ETYPE& e, const LTYPE& l, const MTYPE& m) const {
+ return e * (0.668 + 6.84 / l + 145.0 / pow(l, 3.0)) *
+ (1.0 + 0.5135 * (m - 1.0) * l - 0.0344 * (m - 1.0) * pow(l, 2.0) + 0.0037 * pow((m - 1.0), 3.0) * l) /
+ (1.0 + 0.3220 * (m - 1.0) * l - 0.0223 * (m - 1.0) * pow(l, 2.0));
+ }
+
+ template<typename STYPE, typename LTYPE, typename MTYPE>
+ inline auto get_dc(const STYPE& s, const LTYPE& l, const MTYPE& m) const {
+ return 1E3 * m * (0.2516 + 0.049 * log(l) / log(10.0)) *
+ (1.0 - 0.0681 * (m - 1.0) * l + 0.0029 * (m - 1.0) * pow(l, 2.0)) /
+ (1.0 + 0.0462 * (m - 1.0) * l - 0.0019 * (m - 1.0) * pow(l, 2.0)) * 1E27 / (NAvo * pow(s, 3.0)) / m;
+ }
+
+
+ // Approximate the mixture with the one fluid model
+ template<typename TTYPE, typename RHOTYPE, typename MoleFracType>
+ auto one_fluid(TTYPE& Tstar, RHOTYPE& rhostar, MoleFracType& molefrac) const {
+
+ using resulttype = std::common_type_t<decltype(Tstar), decltype(molefrac[0]), decltype(rhostar)>;
+ resulttype l = 0.0;
+ resulttype s = 0.0;
+ resulttype e = 0.0;
+ resulttype m = 0.0;
+
+ std::tie(s, e, l, m) = combining_rules_one_fluid(rhostar, molefrac, fld, is_sphere);
+
+ resulttype tc = get_tc(e, l, m);
+ resulttype dc = get_dc(s, l, m);
+
+ resulttype tau = tc / Tstar;
+ resulttype delta = rhostar / dc;
+
+ // Calculate coefficients
+ std::vector< resulttype> n_pol;
+ std::vector< resulttype> n_exp;
+ std::vector< resulttype> n_gbs;
+ for (size_t i = 0; i < t_pol.size(); i++) { n_pol.push_back(c1_pol[i] + c2_pol[i] / l); }
+ for (size_t i = 0; i < t_exp.size(); i++) { n_exp.push_back(c1_exp[i] + c2_exp[i] / l); }
+ for (size_t i = 0; i < t_gbs.size(); i++) { n_gbs.push_back(c1_gbs[i] + c2_gbs[i] / l); }
+
+
+ std::vector< resulttype> n_pol_m;
+ std::vector< resulttype> n_exp_m;
+ std::vector< resulttype> n_gbs_m;
+ for (size_t i = 0; i < t_pol_m.size(); i++) { n_pol_m.push_back(d1_pol[i] + d2_pol[i] * (1.0 - m) + d3_pol[i] * (1.0 - m) * (1.0 - m)); }
+ for (size_t i = 0; i < t_exp_m.size(); i++) { n_exp_m.push_back(d1_exp[i] + d2_exp[i] * (1.0 - m) + d3_exp[i] * (1.0 - m) * (1.0 - m)); }
+ for (size_t i = 0; i < t_gbs_m.size(); i++) { n_gbs_m.push_back(d1_gbs[i] + d2_gbs[i] * (1.0 - m) + d3_gbs[i] * (1.0 - m) * (1.0 - m)); }
+
+ std::vector< resulttype> pol, exp_, gbs;
+ for (size_t i = 0; i < t_pol.size(); i++) { pol.push_back(n_pol[i] * pow(tau, t_pol[i]) * pow(delta, d_pol[i])); }
+ for (size_t i = 0; i < t_exp.size(); i++) { exp_.push_back(n_exp[i] * pow(tau, t_exp[i]) * pow(delta, d_exp[i]) * exp(-pow(delta, p[i]))); }
+ for (size_t i = 0; i < t_gbs.size(); i++) { gbs.push_back(n_gbs[i] * pow(tau, t_gbs[i]) * pow(delta, d_gbs[i]) * exp(-eta[i] * (delta - eps[i]) * (delta - eps[i]) - beta[i] * (tau - gam[i]) * (tau - gam[i]))); }
+
+ std::vector< resulttype> pol_m, exp_m, gbs_m;
+ for (size_t i = 0; i < t_pol_m.size(); i++) { pol_m.push_back(n_pol_m[i] * pow(tau, t_pol_m[i]) * pow(delta, d_pol_m[i])); }
+ for (size_t i = 0; i < t_exp_m.size(); i++) { exp_m.push_back(n_exp_m[i] * pow(tau, t_exp_m[i]) * pow(delta, d_exp_m[i]) * exp(-pow(delta, p_m[i]))); }
+ for (size_t i = 0; i < t_gbs_m.size(); i++) { gbs_m.push_back(n_gbs_m[i] * pow(tau, t_gbs_m[i]) * pow(delta, d_gbs_m[i]) * exp(-eta_m[i] * (delta - eps_m[i]) * (delta - eps_m[i]) - beta_m[i] * (tau - gam_m[i]) * (tau - gam_m[i]))); }
+
+ auto alpha_r_sphere = std::reduce(pol.begin(), pol.end()) + std::reduce(exp_.begin(), exp_.end()) + std::reduce(gbs.begin(), gbs.end());
+ auto alpha_r_el = std::reduce(pol_m.begin(), pol_m.end()) + std::reduce(exp_m.begin(), exp_m.end()) + std::reduce(gbs_m.begin(), gbs_m.end());
+ resulttype alpha_r_all = 0.0;
+ if (is_sphere) {
+ alpha_r_all = alpha_r_sphere;
+ }
+ else
+ {
+ //alpha_r_all = m * alpha_r_sphere + (1.0 - m) * alpha_r_el;
+ alpha_r_all = alpha_r_sphere + (1.0 - m) * alpha_r_el;
+ }
+
+ return alpha_r_all;
+ }
+
+ template<typename TTYPE, typename RHOTYPE, typename MoleFracType>
+ auto corresponding_states(const TTYPE& Tstar, const RHOTYPE& rhostar, const MoleFracType& molefrac, const std::string combination_rule) const {
+
+ using resulttype = std::common_type_t<decltype(Tstar), decltype(molefrac[0]), decltype(rhostar)>;
+
+ auto ncomp = nr_fld;
+ // get pure fluid reducing
+ std::vector<resulttype> tc(nr_fld, 0.0);
+ std::vector<resulttype> dc(nr_fld, 0.0);
+ resulttype TC_MIX = 0.0;
+ resulttype DC_MIX = 0.0;
+
+ for (size_t i = 0; i < nr_fld; i++) {
+ tc[i] = get_tc(fld[i].epsilons, fld[i].lambdas, fld[i].m);
+ dc[i] = get_dc(fld[i].sigmas, fld[i].lambdas, fld[i].m);
+ }
+
+ std::vector<std::vector<resulttype>> gamma_t(ncomp, std::vector<resulttype>(ncomp, 1.0));
+ std::vector<std::vector<resulttype>> gamma_v(ncomp, std::vector<resulttype>(ncomp, 1.0));
+ if (combining_rule.compare("Lorentz-Berthelot") == 0) {
+ TC_MIX = get_tred_mix(tc, molefrac, gamma_t);
+ DC_MIX = get_vred_mix(dc, molefrac, gamma_v);
+ }
+ else if (combining_rule.compare("Linear") == 0) {
+ for (auto i = 0; i < ncomp; i++) {
+ for (auto j = 0; j < ncomp; j++) {
+ auto dom = pow(1.0 / pow(dc[i], 1.0 / 3.0) + 1.0 / pow(dc[j], 1.0 / 3.0), 3.0);
+ auto nom = 1.0 / dc[i] + 1.0 / dc[j];
+ gamma_v[i][j] = 4.0 * nom / dom;
+ gamma_t[i][j] = 0.5 * (tc[i] + tc[j]) / sqrt(tc[i] * tc[j]);
+ };
+ };
+ TC_MIX = get_tred_mix(tc, molefrac, gamma_t);
+ DC_MIX = get_vred_mix(dc, molefrac, gamma_v);
+ }
+
+ std::vector<std::vector<resulttype>> n_pol(ncomp, std::vector<resulttype>(t_pol.size(), 0.0));
+ std::vector<std::vector<resulttype>> n_exp(ncomp, std::vector<resulttype>(t_exp.size(), 0.0));
+ std::vector<std::vector<resulttype>> n_gbs(ncomp, std::vector<resulttype>(t_gbs.size(), 0.0));
+
+ // Calculate coefficients
+ for (size_t i = 0; i < ncomp; i++) {
+ for (size_t j = 0; j < t_pol.size(); j++) {
+ n_pol[i][j] = c1_pol[j] + c2_pol[j] / fld[i].lambdas;
+ }
+ };
+ for (size_t i = 0; i < ncomp; i++) {
+ for (size_t j = 0; j < t_exp.size(); j++) {
+ n_exp[i][j] = c1_exp[j] + c2_exp[j] / fld[i].lambdas;
+ }
+ };
+ for (size_t i = 0; i < ncomp; i++) {
+ for (size_t j = 0; j < t_gbs.size(); j++) {
+ n_gbs[i][j] = c1_gbs[j] + c2_gbs[j] / fld[i].lambdas;
+ };
+ };
+
+ std::vector<std::vector<resulttype>> n_pol_m(ncomp, std::vector<resulttype>(t_pol_m.size(), 0.0));
+ std::vector<std::vector<resulttype>> n_exp_m(ncomp, std::vector<resulttype>(t_exp_m.size(), 0.0));
+
+ for (size_t i = 0; i < ncomp; i++) {
+ for (size_t j = 0; j < t_pol_m.size(); j++) {
+ n_pol_m[i][j] = d1_pol[j] + d2_pol[j] * (1.0 - fld[i].m) * (1.0 - fld[i].m) + d3_pol[j] * (1.0 - fld[i].m) * (1.0 - fld[i].m) * (1.0 - fld[i].m);
+ };
+ };
+
+ for (size_t i = 0; i < ncomp; i++) {
+ for (size_t j = 0; j < t_exp_m.size(); j++) {
+ n_exp_m[i][j] = d1_exp[j] + d2_exp[j] * (1.0 - fld[i].m) * (1.0 - fld[i].m) + d3_exp[j] * (1.0 - fld[i].m) * (1.0 - fld[i].m) * (1.0 - fld[i].m);
+ };
+ };
+
+ std::vector<resulttype> tau;
+ std::vector<resulttype> delta;
+
+ for (size_t i = 0; i < ncomp; i++) {
+ tau.push_back(TC_MIX / Tstar);
+ delta.push_back(rhostar / DC_MIX);
+ };
+
+ std::vector<resulttype> pol(ncomp, 0.0);
+ std::vector<resulttype> exp_(ncomp, 0.0);
+ std::vector<resulttype> gbs(ncomp, 0.0);
+ std::vector<resulttype> pol_m(ncomp, 0.0);
+ std::vector<resulttype> exp_m(ncomp, 0.0);
+
+ // pure fluid contribution
+ for (size_t i = 0; i < ncomp; i++) {
+ pol[i] = 0.0;
+ for (size_t j = 0; j < t_pol.size(); j++) {
+ pol[i] += n_pol[i][j] * pow(tau[i], t_pol[j]) * pow(delta[i], d_pol[j]);
+ };
+ };
+
+ for (size_t i = 0; i < ncomp; i++) {
+ exp_[i] = 0.0;
+ for (size_t j = 0; j < t_exp.size(); j++) {
+ exp_[i] += n_exp[i][j] * pow(tau[i], t_exp[j]) * pow(delta[i], d_exp[j]) * exp(-pow(delta[i], p[j]));
+ };
+ };
+
+ for (size_t i = 0; i < ncomp; i++) {
+ gbs[i] = 0.0;
+ for (size_t j = 0; j < t_gbs.size(); j++) {
+ gbs[i] += n_gbs[i][j] * pow(tau[i], t_gbs[j]) * pow(delta[i], d_gbs[j]) * exp(-eta[j] * (delta[i] - eps[j]) * (delta[i] - eps[j]) - beta[j] * (tau[i] - gam[j]) * (tau[i] - gam[j]));
+ };
+ };
+
+ if (!(is_sphere)) {
+ for (size_t i = 0; i < ncomp; i++) {
+ pol_m[i] = 0.0;
+ for (size_t j = 0; j < t_pol_m.size(); j++) {
+ pol_m[i] += n_pol_m[i][j] * pow(tau[i], t_pol_m[j]) * pow(delta[i], d_pol_m[j]);
+ };
+ };
+
+ for (size_t i = 0; i < ncomp; i++) {
+ exp_m[i] = 0.0;
+ for (size_t j = 0; j < t_exp_m.size(); j++) {
+ exp_m[i] += n_exp_m[i][j] * pow(tau[i], t_exp_m[j]) * pow(delta[i], d_exp_m[j]) * exp(-pow(delta[i], p_m[j]));
+ };
+ };
+ }
+
+ // Sum up everything
+ resulttype alpha_r = 0.0;
+ resulttype alpha_r_sphere = 0.0;
+ resulttype alpha_r_elongated = 0.0;
+
+ for (size_t i = 0; i < ncomp; i++) {
+ if (fld[i].is_sphere) {
+ // Single sphere part of component i
+ alpha_r_sphere = pol[i] + exp_[i] + gbs[i];
+ alpha_r += molefrac[i] * (fld[i].m * alpha_r_sphere);
+ }
+ else {
+ // Single elongation part of component i
+ alpha_r_sphere = pol[i] + exp_[i] + gbs[i];
+ alpha_r_elongated = pol_m[i] + exp_m[i];
+ alpha_r += molefrac[i] * (fld[i].m * alpha_r_sphere + (1.0 - fld[i].m) * alpha_r_elongated);
+ }
+ };
+
+ return alpha_r;
+ }
+
+
+
+ // Input is temperature in K, density in mol/m^3 and molefractions
+ template<typename TTYPE, typename RHOTYPE, typename MoleFracType>
+ auto alphar(const TTYPE& Tstar, const RHOTYPE& rhostar, const MoleFracType& molefrac) const {
+ using resulttype = std::common_type_t<decltype(Tstar), decltype(molefrac[0]), decltype(rhostar)>;
+
+ resulttype alpha_r_all = 0.0;
+
+ if (combining_rule.compare("One-Fluid") == 0) {
+ alpha_r_all = one_fluid(Tstar, rhostar, molefrac);
+ }
+ else
+ {
+ alpha_r_all = corresponding_states(Tstar, rhostar, molefrac, combining_rule);
+ }
+
+ return forceeval(alpha_r_all);
+ }
+
+ };
+
+
+ }
+}
+
+
+// class MieElongMix2 {
+
+// private:
+
+// using EArray6 = Eigen::Array<double, 6, 1>;
+// using EArray4 = Eigen::Array<double, 4, 1>;
+
+
+// const EArray6 c1_pol = (EArray6() << -0.0192944, 1.38, -2.2653, 1.6291, -1.974, 0.40412).finished();
+// const EArray6 c1_exp = (EArray6() << 0.1845, -0.3227, 1.1351, 2.232, -2.344, -0.4238).finished();
+// const EArray4 c1_gbs = (EArray4() << -4.367, 0.0371, 1.3895, 2.835).finished();
+// const EArray6 c2_pol = (EArray6() << 0.26021, -5.525, 8.329, -19.492, 25.8, -3.8133).finished();
+// const EArray6 c2_exp = (EArray6() << -5.05, 2.7842, -9.523, -30.383, 17.902, 2.2264).finished();
+// const EArray4 c2_gbs = (EArray4() << 48.445, -5.506, -11.643, -24.36).finished();
+// const EArray6 t_pol = (EArray6() << 1, 0.236, 0.872, 0.313, 0.407, 0.703).finished();
+// const EArray6 t_exp = (EArray6() << 1.78, 2.99, 2.866, 1.2, 3.06, 1.073).finished();
+// const EArray4 t_gbs = (EArray4() << 1.50, 1.03, 4.02, 1.57).finished();
+// const EArray6 d_pol = (EArray6() << 4, 1, 1, 2, 2, 3).finished();
+// const EArray6 d_exp = (EArray6() << 1, 1, 3, 2, 2, 5).finished();
+// const EArray4 d_gbs = (EArray4() << 2, 3, 2, 2).finished();
+// const EArray6 p = (EArray6() << 1, 2, 2, 1, 2, 1).finished();
+// const EArray4 eta = (EArray4() << 0.362, 0.313, 1.17, 0.957).finished();
+// const EArray4 beta = (EArray4() << 0.0761, 0.143, 0.63, 1.32).finished();
+// const EArray4 gam = (EArray4() << 1.55, -0.0826, 1.505, 1.07).finished();
+// const EArray4 eps = (EArray4() << -1, -1, -0.195, -0.287).finished();
+
+// // elongation part
+// const EArray6 d1_pol = (EArray6() << -0.0033717, -0.93963, 1.0403, -0.11078, -0.0025363, 0.0029514).finished();
+// const EArray6 d2_pol = (EArray6() << 0.0056524, -1.0549, 1.2465, -0.22323, 0.15836, 0.11005).finished();
+// const EArray6 d3_pol = (EArray6() << -0.0032422, -0.69384, 0.80325, 0.26957, -0.069081, -0.0026843).finished();
+// const EArray6 d1_exp = (EArray6() << -1.0436, 0.25891, -0.040776, -0.89272, -0.12408, -0.053926).finished();
+// const EArray6 d2_exp = (EArray6() << -0.27931, -0.27319, 0.1812, -1.2567, 0.71416, -0.040811).finished();
+// const EArray6 d3_exp = (EArray6() << -0.26861, -0.032778, 0.17816, -1.3529, 0.58892, 0.0034168).finished();
+
+// const EArray6 t_pol_m = (EArray6() << 1, 0.684125778408247, 0.431293949327048, 0.686848959577935, 0.39116753396418, 0.755900521521362).finished();
+// const EArray6 t_exp_m = (EArray6() << 0.994097378429697, 0.81534766545966, 1.96973994708541, 1.35262796041634, 3, 1.41415130674305).finished();
+// const EArray6 d_pol_m = (EArray6() << 4, 1, 1, 2, 2, 3).finished();
+// const EArray6 d_exp_m = (EArray6() << 1, 1, 3, 2, 2, 5).finished();
+// const EArray6 p_m = (EArray6() << 1, 2, 2, 1, 2, 1).finished();
+
+// double m_lambda_a, ms;
+// std::vector<double> lambdas, epsilons, sigmas, ms_s;
+// public:
+
+// //MieElong(double lambda_a, double ms, double sigma, double eps) { // Repulsive exponent, segment number ms
+// MieElongMix2(const std::string& path, const std::vector<std::string>& fluids) {
+// std::string filepath = std::filesystem::is_regular_file(path) ? path : path;
+// nlohmann::json j = load_a_JSON_file(filepath);
+
+// for (auto f : fluids)
+// {
+// lambdas.push_back(j.at("mie").at(f).at("lambda"));
+// epsilons.push_back(j.at("mie").at(f).at("epsilon"));
+// sigmas.push_back(j.at("mie").at(f).at("sigma"));
+// ms_s.push_back(j.at("mie").at(f).at("segment"));
+// };
+
+// };
+
+
+// // We are in "simulation units", so R is 1.0, and T and rho that
+// // go into alphar are actually T^* and rho^*
+// template<typename MoleFracType>
+// double R(const MoleFracType&) const { return NAvo * kBoltz; }
+
+// // Input is temperature in K, density in mol/m^3 and molefractions
+// template<typename TTYPE, typename RHOTYPE, typename MoleFracType>
+// auto alphar(const TTYPE& Tstar, const RHOTYPE& rhostar, const MoleFracType& molefrac) const {
+
+// auto ncomp = ms_s.size();
+// using resulttype = std::common_type_t<decltype(Tstar), decltype(molefrac[0]), decltype(rhostar)>;
+
+
+// std::vector<vector<resulttype>> sigma_ij(ncomp, vector<resulttype>(ncomp, rhostar));
+// std::vector<vector<resulttype>> epsilon_ij(ncomp, vector<resulttype>(ncomp, rhostar));
+// std::vector<vector<resulttype>> lambda_ij(ncomp, vector<resulttype>(ncomp, rhostar));
+// std::vector<vector<resulttype>> gamma_t(ncomp, vector<resulttype>(ncomp, rhostar));
+// std::vector<vector<resulttype>> gamma_v(ncomp, vector<resulttype>(ncomp, rhostar));
+
+
+// for (auto i = 0; i < ncomp; i++) {
+// for (auto j = 0; j < ncomp; j++) {
+// lambda_ij[i][j] = (lambdas[i] + lambdas[j]) / 2.0;
+// sigma_ij[i][j] = (sigmas[i] + sigmas[j]) / 2.0;
+// epsilon_ij[i][j] = sqrt(epsilons[i] * epsilons[j]);
+// auto dom = pow(1.0 / pow(sigmas[i], 1.0 / 3.0) + 1.0 / pow(sigmas[j], 1.0 / 3.0), 3.0);
+// auto nom = 1.0 / sigmas[i] + 1.0 / sigmas[j];
+// gamma_v[i][j] = 1.0;//4.0 * nom / dom;
+// gamma_t[i][j] = 1.0;//0.5 * (epsilons[i] + epsilons[j]) / sqrt(epsilons[i] * epsilons[j]);
+// };
+// };
+
+// std::vector<resulttype> TC;
+// std::vector<resulttype> DC;
+// for (size_t i = 0; i < lambdas.size(); i++) {
+// TC.push_back(get_tc(epsilons[i], lambdas[i], ms_s[i]));
+// DC.push_back(get_dc(sigmas[i], lambdas[i], ms_s[i]));
+// };
+
+
+// std::vector<vector<resulttype>> n_pol(ncomp, vector<resulttype>(t_pol.size(), 0.0));
+// std::vector<vector<resulttype>> n_exp(ncomp, vector<resulttype>(t_exp.size(), 0.0));
+// std::vector<vector<resulttype>> n_gbs(ncomp, vector<resulttype>(t_gbs.size(), 0.0));
+// // Calculate coefficients
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_pol.size(); j++) {
+// n_pol[i][j] = c1_pol[j] + c2_pol[j] / lambdas[i];
+// }
+// };
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_exp.size(); j++) {
+// n_exp[i][j] = c1_exp[j] + c2_exp[j] / lambdas[i];
+// }
+// };
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_gbs.size(); j++) {
+// n_gbs[i][j] = c1_gbs[j] + c2_gbs[j] / lambdas[i];
+// };
+// };
+
+// std::vector<vector<resulttype>> n_pol_m(ncomp, vector<resulttype>(t_pol_m.size(), 0.0));
+// std::vector<vector<resulttype>> n_exp_m(ncomp, vector<resulttype>(t_exp_m.size(), 0.0));
+
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_pol_m.size(); j++) {
+// n_pol_m[i][j] = d1_pol[j] + d2_pol[j] * (1.0 - ms_s[i]) * (1.0 - ms_s[i]) + d3_pol[j] * (1.0 - ms_s[i]) * (1.0 - ms_s[i]) * (1.0 - ms_s[i]);
+// };
+// };
+
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_exp_m.size(); j++) {
+// n_exp_m[i][j] = d1_exp[j] + d2_exp[j] * (1.0 - ms_s[i]) * (1.0 - ms_s[i]) + d3_exp[j] * (1.0 - ms_s[i]) * (1.0 - ms_s[i]) * (1.0 - ms_s[i]);
+// };
+// };
+
+// std::vector<resulttype> tau;
+// std::vector<resulttype> delta;
+// for (size_t i = 0; i < ncomp; i++) {
+// tau.push_back(get_tred_mix(TC, molefrac, gamma_t) / Tstar);
+// delta.push_back(rhostar / get_vred_mix(DC, molefrac, gamma_v));
+// };
+
+// std::vector<resulttype> pol(ncomp, 0.0);
+// std::vector<resulttype> exp_(ncomp, 0.0);
+// std::vector<resulttype> gbs(ncomp, 0.0);
+// std::vector<resulttype> pol_m(ncomp, 0.0);
+// std::vector<resulttype> exp_m(ncomp, 0.0);
+// for (size_t i = 0; i < ncomp; i++) {
+// pol[i] = 0.0;
+// for (size_t j = 0; j < t_pol.size(); j++) {
+// pol[i] += n_pol[i][j] * pow(tau[i], t_pol[j]) * pow(delta[i], d_pol[j]);
+// };
+// };
+// for (size_t i = 0; i < ncomp; i++) {
+// exp_[i] = 0.0;
+// for (size_t j = 0; j < t_exp.size(); j++) {
+// exp_[i] += n_exp[i][j] * pow(tau[i], t_exp[j]) * pow(delta[i], d_exp[j]) * exp(-pow(delta[i], p[j]));
+// };
+// };
+// for (size_t i = 0; i < ncomp; i++) {
+// gbs[i] = 0.0;
+// for (size_t j = 0; j < t_gbs.size(); j++) {
+// gbs[i] += n_gbs[i][j] * pow(tau[i], t_gbs[j]) * pow(delta[i], d_gbs[j]) * exp(-eta[j] * (delta[i] - eps[j]) * (delta[i] - eps[j]) - beta[j] * (tau[i] - gam[j]) * (tau[i] - gam[j]));
+// };
+// };
+
+// for (size_t i = 0; i < ncomp; i++) {
+// pol_m[i] = 0.0;
+// for (size_t j = 0; j < t_pol_m.size(); j++) {
+// pol_m[i] += n_pol_m[i][j] * pow(tau[i], t_pol_m[j]) * pow(delta[i], d_pol_m[j]);
+// };
+// };
+// for (size_t i = 0; i < ncomp; i++) {
+// exp_m[i] = 0.0;
+// for (size_t j = 0; j < t_exp_m.size(); j++) {
+// exp_m[i] += n_exp_m[i][j] * pow(tau[i], t_exp_m[j]) * pow(delta[i], d_exp_m[j]) * exp(-pow(delta[i], p_m[j]));
+// };
+// };
+
+
+// resulttype alpha_r_sphere = 0.0;
+// resulttype alpha_r_el = 0.0;
+// resulttype alpha_r_all = 0.0;
+// for (size_t i = 0; i < ncomp; i++) {
+// alpha_r_sphere = 0.0;
+// alpha_r_el = 0.0;
+// alpha_r_sphere = pol[i] + exp_[i] + gbs[i];
+// alpha_r_el = pol_m[i] + exp_m[i];
+// alpha_r_all += molefrac[i] * (ms_s[i] * alpha_r_sphere + (1.0 - ms_s[i]) * alpha_r_el);
+// };
+// return forceeval(alpha_r_all);
+// }
+
+// };
+
+//};
+
+
+//// new mixing rules try
+//class MieElongMix3 {
+
+//private:
+
+// using EArray6 = Eigen::Array<double, 6, 1>;
+// using EArray4 = Eigen::Array<double, 4, 1>;
+
+
+// const EArray6 c1_pol = (EArray6() << -0.0192944, 1.38, -2.2653, 1.6291, -1.974, 0.40412).finished();
+// const EArray6 c1_exp = (EArray6() << 0.1845, -0.3227, 1.1351, 2.232, -2.344, -0.4238).finished();
+// const EArray4 c1_gbs = (EArray4() << -4.367, 0.0371, 1.3895, 2.835).finished();
+// const EArray6 c2_pol = (EArray6() << 0.26021, -5.525, 8.329, -19.492, 25.8, -3.8133).finished();
+// const EArray6 c2_exp = (EArray6() << -5.05, 2.7842, -9.523, -30.383, 17.902, 2.2264).finished();
+// const EArray4 c2_gbs = (EArray4() << 48.445, -5.506, -11.643, -24.36).finished();
+// const EArray6 t_pol = (EArray6() << 1, 0.236, 0.872, 0.313, 0.407, 0.703).finished();
+// const EArray6 t_exp = (EArray6() << 1.78, 2.99, 2.866, 1.2, 3.06, 1.073).finished();
+// const EArray4 t_gbs = (EArray4() << 1.50, 1.03, 4.02, 1.57).finished();
+// const EArray6 d_pol = (EArray6() << 4, 1, 1, 2, 2, 3).finished();
+// const EArray6 d_exp = (EArray6() << 1, 1, 3, 2, 2, 5).finished();
+// const EArray4 d_gbs = (EArray4() << 2, 3, 2, 2).finished();
+// const EArray6 p = (EArray6() << 1, 2, 2, 1, 2, 1).finished();
+// const EArray4 eta = (EArray4() << 0.362, 0.313, 1.17, 0.957).finished();
+// const EArray4 beta = (EArray4() << 0.0761, 0.143, 0.63, 1.32).finished();
+// const EArray4 gam = (EArray4() << 1.55, -0.0826, 1.505, 1.07).finished();
+// const EArray4 eps = (EArray4() << -1, -1, -0.195, -0.287).finished();
+
+// // elongation part
+// const EArray6 d1_pol = (EArray6() << -0.0033717, -0.93963, 1.0403, -0.11078, -0.0025363, 0.0029514).finished();
+// const EArray6 d2_pol = (EArray6() << 0.0056524, -1.0549, 1.2465, -0.22323, 0.15836, 0.11005).finished();
+// const EArray6 d3_pol = (EArray6() << -0.0032422, -0.69384, 0.80325, 0.26957, -0.069081, -0.0026843).finished();
+// const EArray6 d1_exp = (EArray6() << -1.0436, 0.25891, -0.040776, -0.89272, -0.12408, -0.053926).finished();
+// const EArray6 d2_exp = (EArray6() << -0.27931, -0.27319, 0.1812, -1.2567, 0.71416, -0.040811).finished();
+// const EArray6 d3_exp = (EArray6() << -0.26861, -0.032778, 0.17816, -1.3529, 0.58892, 0.0034168).finished();
+
+// const EArray6 t_pol_m = (EArray6() << 1, 0.684125778408247, 0.431293949327048, 0.686848959577935, 0.39116753396418, 0.755900521521362).finished();
+// const EArray6 t_exp_m = (EArray6() << 0.994097378429697, 0.81534766545966, 1.96973994708541, 1.35262796041634, 3, 1.41415130674305).finished();
+// const EArray6 d_pol_m = (EArray6() << 4, 1, 1, 2, 2, 3).finished();
+// const EArray6 d_exp_m = (EArray6() << 1, 1, 3, 2, 2, 5).finished();
+// const EArray6 p_m = (EArray6() << 1, 2, 2, 1, 2, 1).finished();
+
+// double m_lambda_a, ms;
+// std::vector<double> lambdas, epsilons, sigmas, ms_s;
+//public:
+
+// //MieElong(double lambda_a, double ms, double sigma, double eps) { // Repulsive exponent, segment number ms
+// MieElongMix3(const std::string& path, const std::vector<std::string>& fluids) {
+// std::string filepath = std::filesystem::is_regular_file(path) ? path : path;
+// nlohmann::json j = load_a_JSON_file(filepath);
+
+// for (auto f : fluids)
+// {
+// lambdas.push_back(j.at("mie").at(f).at("lambda"));
+// epsilons.push_back(j.at("mie").at(f).at("epsilon"));
+// sigmas.push_back(j.at("mie").at(f).at("sigma"));
+// ms_s.push_back(j.at("mie").at(f).at("segment"));
+// };
+
+// };
+
+
+// // We are in "simulation units", so R is 1.0, and T and rho that
+// // go into alphar are actually T^* and rho^*
+// template<typename MoleFracType>
+// double R(const MoleFracType&) const { return NAvo * kBoltz; }
+
+// // Input is temperature in K, density in mol/m^3 and molefractions
+// template<typename TTYPE, typename RHOTYPE, typename MoleFracType>
+// auto alphar(const TTYPE& Tstar, const RHOTYPE& rhostar, const MoleFracType& molefrac) const {
+
+// auto ncomp = ms_s.size();
+// using resulttype = std::common_type_t<decltype(Tstar), decltype(molefrac[0]), decltype(rhostar)>;
+
+
+// std::vector<vector<resulttype>> sigma_ij(ncomp, vector<resulttype>(ncomp, rhostar));
+// std::vector<vector<resulttype>> epsilon_ij(ncomp, vector<resulttype>(ncomp, rhostar));
+// std::vector<vector<resulttype>> lambda_ij(ncomp, vector<resulttype>(ncomp, rhostar));
+// std::vector<vector<resulttype>> gamma_t(ncomp, vector<resulttype>(ncomp, rhostar));
+// std::vector<vector<resulttype>> gamma_v(ncomp, vector<resulttype>(ncomp, rhostar));
+
+
+// for (auto i = 0; i < ncomp; i++) {
+// for (auto j = 0; j < ncomp; j++) {
+// lambda_ij[i][j] = (lambdas[i] + lambdas[j]) / 2.0;
+// sigma_ij[i][j] = (sigmas[i] + sigmas[j]) / 2.0;
+// epsilon_ij[i][j] = sqrt(epsilons[i] * epsilons[j]);
+// auto dom = pow(1.0 / pow(sigmas[i], 1.0 / 3.0) + 1.0 / pow(sigmas[j], 1.0 / 3.0), 3.0);
+// auto nom = 1.0 / sigmas[i] + 1.0 / sigmas[j];
+// gamma_v[i][j] = 1.0;//4.0 * nom / dom;
+// gamma_t[i][j] = 1.0;//0.5 * (epsilons[i] + epsilons[j]) / sqrt(epsilons[i] * epsilons[j]);
+// };
+// };
+
+// std::vector<resulttype> TC;
+// std::vector<resulttype> DC;
+// for (size_t i = 0; i < lambdas.size(); i++) {
+// TC.push_back(get_tc(epsilons[i], lambdas[i], ms_s[i]));
+// DC.push_back(get_dc(sigmas[i], lambdas[i], ms_s[i]));
+// };
+
+
+// std::vector<vector<resulttype>> n_pol(ncomp, vector<resulttype>(t_pol.size(), 0.0));
+// std::vector<vector<resulttype>> n_exp(ncomp, vector<resulttype>(t_exp.size(), 0.0));
+// std::vector<vector<resulttype>> n_gbs(ncomp, vector<resulttype>(t_gbs.size(), 0.0));
+// // Calculate coefficients
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_pol.size(); j++) {
+// n_pol[i][j] = c1_pol[j] + c2_pol[j] / lambdas[i];
+// }
+// };
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_exp.size(); j++) {
+// n_exp[i][j] = c1_exp[j] + c2_exp[j] / lambdas[i];
+// }
+// };
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_gbs.size(); j++) {
+// n_gbs[i][j] = c1_gbs[j] + c2_gbs[j] / lambdas[i];
+// };
+// };
+
+// std::vector<vector<resulttype>> n_pol_m(ncomp, vector<resulttype>(t_pol_m.size(), 0.0));
+// std::vector<vector<resulttype>> n_exp_m(ncomp, vector<resulttype>(t_exp_m.size(), 0.0));
+
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_pol_m.size(); j++) {
+// n_pol_m[i][j] = d1_pol[j] + d2_pol[j] * (1.0 - ms_s[i]) * (1.0 - ms_s[i]) + d3_pol[j] * (1.0 - ms_s[i]) * (1.0 - ms_s[i]) * (1.0 - ms_s[i]);
+// };
+// };
+
+// for (size_t i = 0; i < ncomp; i++) {
+// for (size_t j = 0; j < t_exp_m.size(); j++) {
+// n_exp_m[i][j] = d1_exp[j] + d2_exp[j] * (1.0 - ms_s[i]) * (1.0 - ms_s[i]) + d3_exp[j] * (1.0 - ms_s[i]) * (1.0 - ms_s[i]) * (1.0 - ms_s[i]);
+// };
+// };
+
+// std::vector<resulttype> tau;
+// std::vector<resulttype> delta;
+// /*for (size_t i = 0; i < ncomp; i++) {
+// tau.push_back(TC[i] / Tstar);
+// delta.push_back(rhostar / DC[i]);
+// };*/
+// for (size_t i = 0; i < ncomp; i++) {
+// tau.push_back(TC[i] / Tstar);
+// delta.push_back(rhostar / DC[i]);
+// };
+
+// std::vector<resulttype> pol(ncomp, 0.0);
+// std::vector<resulttype> exp_(ncomp, 0.0);
+// std::vector<resulttype> gbs(ncomp, 0.0);
+// std::vector<resulttype> pol_m(ncomp, 0.0);
+// std::vector<resulttype> exp_m(ncomp, 0.0);
+
+// std::vector<resulttype> pol_dep(ncomp - 1, 0.0);
+// std::vector<resulttype> exp_dep(ncomp - 1, 0.0);
+// std::vector<resulttype> gbs_dep(ncomp - 1, 0.0);
+// std::vector<resulttype> pol_m_dep(ncomp - 1, 0.0);
+// std::vector<resulttype> exp_m_dep(ncomp - 1, 0.0);
+
+
+// // pure fluid contribution
+// for (size_t i = 0; i < ncomp; i++) {
+// pol[i] = 0.0;
+// for (size_t j = 0; j < t_pol.size(); j++) {
+// pol[i] += n_pol[i][j] * pow(tau[i], t_pol[j]) * pow(delta[i], d_pol[j]);
+// };
+// };
+
+// for (size_t i = 0; i < ncomp; i++) {
+// exp_[i] = 0.0;
+// for (size_t j = 0; j < t_exp.size(); j++) {
+// exp_[i] += n_exp[i][j] * pow(tau[i], t_exp[j]) * pow(delta[i], d_exp[j]) * exp(-pow(delta[i], p[j]));
+// };
+// };
+
+// for (size_t i = 0; i < ncomp; i++) {
+// gbs[i] = 0.0;
+// for (size_t j = 0; j < t_gbs.size(); j++) {
+// gbs[i] += n_gbs[i][j] * pow(tau[i], t_gbs[j]) * pow(delta[i], d_gbs[j]) * exp(-eta[j] * (delta[i] - eps[j]) * (delta[i] - eps[j]) - beta[j] * (tau[i] - gam[j]) * (tau[i] - gam[j]));
+// };
+// };
+
+// for (size_t i = 0; i < ncomp; i++) {
+// pol_m[i] = 0.0;
+// for (size_t j = 0; j < t_pol_m.size(); j++) {
+// pol_m[i] += n_pol_m[i][j] * pow(tau[i], t_pol_m[j]) * pow(delta[i], d_pol_m[j]);
+// };
+// };
+
+// for (size_t i = 0; i < ncomp; i++) {
+// exp_m[i] = 0.0;
+// for (size_t j = 0; j < t_exp_m.size(); j++) {
+// exp_m[i] += n_exp_m[i][j] * pow(tau[i], t_exp_m[j]) * pow(delta[i], d_exp_m[j]) * exp(-pow(delta[i], p_m[j]));
+// };
+// };
+
+// // depature contribution
+// /*for (size_t i = 0; i < ncomp - 1; i++) {
+// pol_dep[i] = 0.0;
+// for (size_t j = 0; j < t_pol.size(); j++) {
+// pol_dep[i] += (n_pol[i][j] + n_pol[i + 1][j]) / 2.0 * pow(tau[i], t_pol[j]) * pow(delta[i], d_pol[j]);
+// };
+// };
+
+// for (size_t i = 0; i < ncomp - 1; i++) {
+// exp_dep[i] = 0.0;
+// for (size_t j = 0; j < t_exp.size(); j++) {
+// exp_dep[i] += (n_exp[i][j] + n_exp[i + 1][j]) / 2.0 * pow(tau[i], t_exp[j]) * pow(delta[i], d_exp[j]) * exp(-pow(delta[i], p[j]));
+// };
+// };
+
+// for (size_t i = 0; i < ncomp - 1; i++) {
+// gbs_dep[i] = 0.0;
+// for (size_t j = 0; j < t_gbs.size(); j++) {
+// gbs_dep[i] += (n_gbs[i][j] + n_gbs[i + 1][j]) / 2.0 * pow(tau[i], t_gbs[j]) * pow(delta[i], d_gbs[j]) * exp(-eta[j] * (delta[i] - eps[j]) * (delta[i] - eps[j]) - beta[j] * (tau[i] - gam[j]) * (tau[i] - gam[j]));
+// };
+// };
+
+// for (size_t i = 0; i < ncomp - 1; i++) {
+// pol_m_dep[i] = 0.0;
+// for (size_t j = 0; j < t_pol_m.size(); j++) {
+// pol_m_dep[i] += (n_pol_m[i][j] + n_pol_m[i + 1][j]) / 2.0 * pow(tau[i], t_pol_m[j]) * pow(delta[i], d_pol_m[j]);
+// };
+// };
+// for (size_t i = 0; i < ncomp - 1; i++) {
+// exp_m_dep[i] = 0.0;
+// for (size_t j = 0; j < t_exp_m.size(); j++) {
+// exp_m_dep[i] += (n_exp_m[i][j] + n_exp_m[i + 1][j]) / 2.0 * pow(tau[i], t_exp_m[j]) * pow(delta[i], d_exp_m[j]) * exp(-pow(delta[i], p_m[j]));
+// };
+// };*/
+// resulttype m_mix = 0.0;
+// for (size_t i = 0; i < ncomp; i++) {
+// m_mix += ms_s[i] * molefrac[i];
+// }
+
+// resulttype alpha_r_sphere = 0.0;
+// resulttype alpha_r_el = 0.0;
+// resulttype alpha_r_all = 0.0;
+
+// std::vector<resulttype> sphere_part(ncomp, 0.0);
+// std::vector<resulttype> elongation_part(ncomp, 0.0);
+// resulttype sphere_mix = 0.0;
+// resulttype elongation_mix = 0.0;
+// resulttype cross_mix = 0.0;
+
+// for (size_t i = 0; i < ncomp; i++) {
+// // Single sphere part of component i
+// sphere_part[i] = pol[i] + exp_[i] + gbs[i];
+
+// // Single elongation part of component j
+// elongation_part[i] = pol_m[i] + exp_m[i];
+
+// sphere_mix += ms_s[i] * molefrac[i] * sphere_part[i];
+
+// };
+
+// for (size_t i = 0; i < ncomp - 1; i++) {
+// for (size_t j = i + 1; j < ncomp; j++) {
+// /*cross_mix += (molefrac[i] * molefrac[j] * ms_s[i]*(1.0 - ms_s[j])*sphere_part[i]*elongation_part[j] +
+// ms_s[j]*(1.0 - ms_s[i])*sphere_part[j]*elongation_part[i] * molefrac[i] * molefrac[j]);*/
+// elongation_mix = (1.0 - ms_s[i]) * molefrac[i] * elongation_part[i] + (1.0 - ms_s[j]) * molefrac[j] * elongation_part[j];
+// }
+// }
+
+// alpha_r_all = sphere_mix + elongation_mix + cross_mix;
+// return forceeval(alpha_r_all);
+// }
+
+//};
+
+
diff --git a/include/teqp/types.hpp b/include/teqp/types.hpp
index bcebf21..1c296e0 100644
--- a/include/teqp/types.hpp
+++ b/include/teqp/types.hpp
@@ -213,6 +213,12 @@ namespace teqp {
}
#endif
+ template<typename T,typename B>
+ auto pow(const B& BASE, const T& POWER) {
+ return exp(POWER*log(BASE));
+ }
+
+
}; // namespace teqp
#endif
diff --git a/interface/CPP/teqp_impl_factory.cpp b/interface/CPP/teqp_impl_factory.cpp
index 8c0c3f2..ee0909a 100644
--- a/interface/CPP/teqp_impl_factory.cpp
+++ b/interface/CPP/teqp_impl_factory.cpp
@@ -31,7 +31,7 @@ namespace teqp {
{"2CLJF-Dipole", [](const nlohmann::json& spec){ return make_owned(twocenterljf::build_two_center_model_dipole(spec.at("author"), spec.at("L^*"), spec.at("(mu^*)^2")));}},
{"2CLJF-Quadrupole", [](const nlohmann::json& spec){ return make_owned(twocenterljf::build_two_center_model_quadrupole(spec.at("author"), spec.at("L^*"), spec.at("(mu^*)^2")));}},
{"IdealHelmholtz", [](const nlohmann::json& spec){ return make_owned(IdealHelmholtz(spec));}},
-
+ {"MieElong", [](const nlohmann::json& spec) { return make_owned(Mie::MieElong(spec.at("model_path"),spec.at("components"),spec.at("combination"))); }},
// Implemented in its own compilation unit to help with compilation time
{"SAFT-VR-Mie", [](const nlohmann::json& spec){ return make_SAFTVRMie(spec); }}
};
--
GitLab