diff --git a/include/teqp/models/saft/correlation_integrals.hpp b/include/teqp/models/saft/correlation_integrals.hpp index 8cfa2cb..0812be9 100644 --- a/include/teqp/models/saft/correlation_integrals.hpp +++ b/include/teqp/models/saft/correlation_integrals.hpp @@ -105,6 +105,21 @@ class LuckasKIntegral{ } }; +class KLuckasSidecar{ +public: + + const int n1, n2; + const std::array a; + double a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31, a32, a33; + + KLuckasSidecar(const int n1, const int n2) : n1(n1), n2(n2), a(Luckas_K_coeffs.at({n1, n2})){ + a00 = a[0]; a01 = a[1]; a02 = a[2]; a03 = a[3]; + a10 = a[4]; a11 = a[5]; a12 = a[6]; a13 = a[7]; + a20 = a[8]; a21 = a[9]; a22 = a[10]; a23 = a[11]; + a30 = a[12]; a31 = a[13]; a32 = a[14]; a33 = a[15]; + } +}; + static const std::map> GubbinsTwu_J_coeffs = { {4, {-0.257431, 0.439229, 0.414783, -0.457019, -0.145520, 0.299666}}, {5, {-0.396724, 0.690721, 0.628935, -0.652622, -0.201462, -0.231635}}, @@ -197,7 +212,21 @@ class GubbinsTwuKIntegral{ } }; - +class KGubbinsTwuSidecar{ +public: + + const int n1, n2; + const std::array a; + double A, B, C, D, E, F; + double sign_term; + + KGubbinsTwuSidecar(const int n1, const int n2) : n1(n1), n2(n2), a(GubbinsTwu_K_coeffs.at({n1, n2})){ + A = a[0]; B = a[1]; C = a[2]; D = a[3]; E = a[4]; F = a[5]; + + // The {334, 445} term has opposite sign to the others + sign_term = ((n1==334) && (n2==445) ? -1 : 1); + } +}; //Type R^2 a[0, 0] a[0, 1] a[0,2] a[0,3] a[1, 0] a[1, 1] a[1,2] a[1,3] a[2, 0] a[2, 1] a[2,2] a[2,3] a[3, 0] a[3, 1] a[3,2] a[3,3] a[4, 0] a[4, 1] a[4,2] a[4,3] b[0, 0] b[0, 1] b[0,2] b[1, 0] b[1, 1] b[1,2] b[2, 0] b[2, 1] b[2,2] b[3, 0] b[3, 1] b[3,2] b[4, 0] b[4, 1] b[4,2] static const std::map> Gottschalk_J_coeffs = { @@ -421,6 +450,28 @@ class GottschalkKIntegral{ } }; + +class KGottschalkSidecar{ +private: + std::tuple int2key(int i){ + if (i < 222){ + throw teqp::InvalidArgument(""); + } + if (i > 999){ + throw teqp::InvalidArgument(""); + } + return std::make_tuple(i / 100 % 10, i / 10 % 10, i % 10); + } +public: + const std::tuple k1, k2; + const std::array abc; + + /// Constructor taking two tuples of ints + KGottschalkSidecar(std::tuple k1, std::tuple k2) : k1(k1), k2(k2), abc(Gottschalk_K_coeffs.at({k1, k2})){} + /// Constructor taking two three digit integers, each of which are split into tuples of ints + KGottschalkSidecar(int k1, int k2) : k1(int2key(k1)), k2(int2key(k2)), abc(Gottschalk_K_coeffs.at({this->k1, this->k2})){} +}; + using JSidecar = std::variant; class JIntegral{ @@ -430,7 +481,7 @@ class JIntegral{ JIntegral(const JSidecar& sidecar) : sidecar(sidecar) {}; template - auto call(const TType& Tstar, const RhoType& rhostar) const -> std::common_type_t{ + auto get_J(const TType& Tstar, const RhoType& rhostar) const -> std::common_type_t{ // Runtime type switching if (std::holds_alternative(sidecar)){ @@ -470,12 +521,67 @@ class JIntegral{ throw teqp::InvalidArgument("don't know what to do with this sidecar"); } } +}; + +using KSidecar = std::variant; + +class KIntegral{ + const KSidecar sidecar; +public: + + KIntegral(const KSidecar& sidecar) : sidecar(sidecar) {}; template - auto get_J(const TType& Tstar, const RhoType& rhostar) const -> std::common_type_t{ - return call(Tstar, rhostar); + auto get_K(const TType& Tstar, const RhoType& rhostar) const -> std::common_type_t{ + + // Runtime type switching + if (std::holds_alternative(sidecar)){ + const auto& d = std::get(sidecar); + double Z_1 = 2.0; + double Z_2 = 3.0; + double Z_3 = 4.0; + RhoType b_0 = d.a00 + d.a10*rhostar + d.a20*rhostar*rhostar + d.a30*rhostar*rhostar*rhostar; + RhoType b_1 = d.a01 + d.a11*rhostar + d.a21*rhostar*rhostar + d.a31*rhostar*rhostar*rhostar; + RhoType b_2 = d.a02 + d.a12*rhostar + d.a22*rhostar*rhostar + d.a32*rhostar*rhostar*rhostar; + RhoType b_3 = d.a03 + d.a13*rhostar + d.a23*rhostar*rhostar + d.a33*rhostar*rhostar*rhostar; + std::common_type_t out = b_0 + b_1*Tstar + b_2*pow(exp(pow(1.0-rhostar/sqrt(2.0),Z_3)), Z_1) + b_3*pow(exp(pow(1.0-rhostar/sqrt(2.0),Z_3)), Z_2); + return out; + } + else if (std::holds_alternative(sidecar)){ + const auto& d = std::get(sidecar); + std::common_type_t out = d.sign_term*exp(d.A*rhostar*rhostar*log(Tstar) + d.B*rhostar*rhostar + d.C*rhostar*log(Tstar) + d.D*rhostar + d.E*log(Tstar) + d.F); + return out; + } + else if (std::holds_alternative(sidecar)){ + const auto& d = std::get(sidecar); + const auto& abc = d.abc; + std::common_type_t summer = 0.0; + int N1 = 8, N2 = 8; // N3 = 24 + + for (auto i = 0; i <= 3; ++i){ + for (auto j = 1; j <= 2; ++j){ + auto I = 2*i + (j-1); // 2 entries in j for each i + summer += abc[I]*pow(rhostar, i)*pow(exp((1.0-rhostar/3.0)/Tstar), j); + } + } + for (auto i = 0; i <= 3; ++i){ + for (auto j = 1; j <= 2; ++j){ + auto I = N1 + 2*i + (j-1); // 2 entries in j for each i + summer += abc[I]*pow(rhostar, i)*pow(exp((1.0-rhostar/3.0)*(1.0-rhostar/3.0)/Tstar), j); + } + } + for (auto i = 0; i <= 5; ++i){ + for (auto j = 0; j <= 3; ++j){ + auto I = N1 + N2 + 4*i + j; // 4 entries in j for each i + summer += abc[I]*pow(rhostar, i)*pow(Tstar, j); + } + } + return summer; + } + else{ + throw teqp::InvalidArgument("don't know what to do with this sidecar"); + } } - }; } diff --git a/include/teqp/models/saft/polar_terms.hpp b/include/teqp/models/saft/polar_terms.hpp index 2fbf14d..805fcf5 100644 --- a/include/teqp/models/saft/polar_terms.hpp +++ b/include/teqp/models/saft/polar_terms.hpp @@ -38,36 +38,30 @@ auto get_Kijk_334445(const KType& Kint, const RhoType& rhostar, const Txy &Tstar }; /** - \tparam KIntegral A type that can be indexed with a two integers a and b to give the K(a,b) integral - - The flexibility was added to include J and K integrals from either Luckas et al. or Gubbins and Twu (or any others following the interface) +The flexibility was added to include J and K integrals from either Luckas et al. or Gubbins and Twu (or any others following the interface) */ -template + class MultipolarContributionGubbinsTwu { public: static constexpr multipolar_argument_spec arg_spec = multipolar_argument_spec::TK_rhoNm3_rhostar_molefractions; private: + const JIntegral J6, J8, J10, J11, J13, J15; + const KIntegral K222_333, K233_344, K334_445, K444_555; + const Eigen::ArrayXd sigma_m, epsilon_over_k, mubar2, Qbar2; const bool has_a_polar; const Eigen::ArrayXd sigma_m3, sigma_m5; - const JIntegral J6{JSidecar{6}}; - const JIntegral J8{JSidecar{8}}; - const JIntegral J10{JSidecar{10}}; - const JIntegral J11{JSidecar{11}}; - const JIntegral J13{JSidecar{13}}; - const JIntegral J15{JSidecar{15}}; - const KIntegral K222_333{222, 333}; - const KIntegral K233_344{233, 344}; - const KIntegral K334_445{334, 445}; - const KIntegral K444_555{444, 555}; - const double epsilon_0 = 8.8541878128e-12; // https://en.wikipedia.org/wiki/Vacuum_permittivity, in F/m, or C^2⋅N^−1⋅m^−2 const double PI_ = static_cast(EIGEN_PI); Eigen::MatrixXd SIGMAIJ, EPSKIJ; multipolar_rhostar_approach approach = multipolar_rhostar_approach::use_packing_fraction; public: - MultipolarContributionGubbinsTwu(const Eigen::ArrayX &sigma_m, const Eigen::ArrayX &epsilon_over_k, const Eigen::ArrayX &mubar2, const Eigen::ArrayX &Qbar2, multipolar_rhostar_approach approach) : sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), mubar2(mubar2), Qbar2(Qbar2), has_a_polar(mubar2.cwiseAbs().sum() > 0 || Qbar2.cwiseAbs().sum() > 0), sigma_m3(sigma_m.pow(3)), sigma_m5(sigma_m.pow(5)), approach(approach) { + + MultipolarContributionGubbinsTwu(const auto& sidecarJ, const auto& sidecarK, const Eigen::ArrayX &sigma_m, const Eigen::ArrayX &epsilon_over_k, const Eigen::ArrayX &mubar2, const Eigen::ArrayX &Qbar2, multipolar_rhostar_approach approach) : + J6(decltype(sidecarJ){6}), J8(decltype(sidecarJ){8}), J10(decltype(sidecarJ){10}), J11(decltype(sidecarJ){11}), J13(decltype(sidecarJ){13}), J15(decltype(sidecarJ){15}), + K222_333(decltype(sidecarK){222,333}), K233_344(decltype(sidecarK){233,344}), K334_445(decltype(sidecarK){334,445}), K444_555(decltype(sidecarK){444,555}), + sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), mubar2(mubar2), Qbar2(Qbar2), has_a_polar(mubar2.cwiseAbs().sum() > 0 || Qbar2.cwiseAbs().sum() > 0), sigma_m3(sigma_m.pow(3)), sigma_m5(sigma_m.pow(5)), approach(approach) { // Check lengths match if (sigma_m.size() != mubar2.size()){ throw teqp::InvalidArgument("bad size of mubar2"); @@ -111,15 +105,15 @@ class MultipolarContributionGubbinsTwu { double sigmaij = SIGMAIJ(i,j); { double dbl = sigma_m3[i]*sigma_m3[j]/powi(sigmaij,3)*mubar2[i]*mubar2[j]; - alpha2_112 += leading*dbl*J6.call(Tstarij, rhostar); + alpha2_112 += leading*dbl*J6.get_J(Tstarij, rhostar); } { double dbl = sigma_m3[i]*sigma_m5[j]/powi(sigmaij,5)*mubar2[i]*Qbar2[j]; - alpha2_123 += leading*dbl*J8.call(Tstarij, rhostar); + alpha2_123 += leading*dbl*J8.get_J(Tstarij, rhostar); } { double dbl = sigma_m5[i]*sigma_m5[j]/powi(sigmaij,7)*Qbar2[i]*Qbar2[j]; - alpha2_224 += leading*dbl*J10.call(Tstarij, rhostar); + alpha2_224 += leading*dbl*J10.get_J(Tstarij, rhostar); } } } @@ -274,31 +268,21 @@ struct PolarizableArrays{ The flexibility was added to include J and K integrals from either Luckas et al. or Gubbins and Twu (or any others following the interface) */ -template class MultipolarContributionGrayGubbins { public: static constexpr multipolar_argument_spec arg_spec = multipolar_argument_spec::TK_rhoNm3_rhostar_molefractions; private: + const JIntegral J6, J8, J10, J11, J13, J15; + const KIntegral K222_333, K233_344, K334_445, K444_555; + const Eigen::ArrayXd sigma_m, epsilon_over_k; Eigen::MatrixXd SIGMAIJ, EPSKIJ; const Eigen::ArrayXd mu, Q, mu2, Q2, Q3; const bool has_a_polar; const Eigen::ArrayXd sigma_m3, sigma_m5; - const JIntegral J6{JSidecar{6}}; - const JIntegral J8{JSidecar{8}}; - const JIntegral J10{JSidecar{10}}; - const JIntegral J11{JSidecar{11}}; - const JIntegral J13{JSidecar{13}}; - const JIntegral J15{JSidecar{15}}; - const KIntegral K222_333{222, 333}; - const KIntegral K233_344{233, 344}; - const KIntegral K334_445{334, 445}; - const KIntegral K444_555{444, 555}; - const double PI_ = static_cast(EIGEN_PI); const double PI3 = PI_*PI_*PI_; - const double epsilon_0 = 8.8541878128e-12; // https://en.wikipedia.org/wiki/Vacuum_permittivity, in F/m, or C^2⋅N^−1⋅m^−2 const double k_e = teqp::constants::k_e; // coulomb constant, with units of N m^2 / C^2 const double k_B = teqp::constants::k_B; // Boltzmann constant, with units of J/K @@ -337,9 +321,12 @@ class MultipolarContributionGrayGubbins { } public: - MultipolarContributionGrayGubbins(const Eigen::ArrayX &sigma_m, const Eigen::ArrayX &epsilon_over_k, const Eigen::MatrixXd& SIGMAIJ, const Eigen::MatrixXd& EPSKIJ, const Eigen::ArrayX &mu, const Eigen::ArrayX &Q, const std::optional& flags) + MultipolarContributionGrayGubbins(const auto& sidecarJ, const auto& sidecarK, const Eigen::ArrayX &sigma_m, const Eigen::ArrayX &epsilon_over_k, const Eigen::MatrixXd& SIGMAIJ, const Eigen::MatrixXd& EPSKIJ, const Eigen::ArrayX &mu, const Eigen::ArrayX &Q, const std::optional& flags) - : sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), SIGMAIJ(SIGMAIJ), EPSKIJ(EPSKIJ), mu(mu), Q(Q), mu2(mu.pow(2)), Q2(Q.pow(2)), Q3(Q.pow(3)), has_a_polar(Q.cwiseAbs().sum() > 0 || mu.cwiseAbs().sum() > 0), sigma_m3(sigma_m.pow(3)), sigma_m5(sigma_m.pow(5)), approach(get_approach(flags)), C3(get_C3(flags)), C3b(get_C3b(flags)), polarizable(get_polarizable(flags)) { + : + J6(decltype(sidecarJ){6}), J8(decltype(sidecarJ){8}), J10(decltype(sidecarJ){10}), J11(decltype(sidecarJ){11}), J13(decltype(sidecarJ){13}), J15(decltype(sidecarJ){15}), + K222_333(decltype(sidecarK){222,333}), K233_344(decltype(sidecarK){233,344}), K334_445(decltype(sidecarK){334,445}), K444_555(decltype(sidecarK){444,555}), + sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), SIGMAIJ(SIGMAIJ), EPSKIJ(EPSKIJ), mu(mu), Q(Q), mu2(mu.pow(2)), Q2(Q.pow(2)), Q3(Q.pow(3)), has_a_polar(Q.cwiseAbs().sum() > 0 || mu.cwiseAbs().sum() > 0), sigma_m3(sigma_m.pow(3)), sigma_m5(sigma_m.pow(5)), approach(get_approach(flags)), C3(get_C3(flags)), C3b(get_C3b(flags)), polarizable(get_polarizable(flags)) { // Check lengths match if (sigma_m.size() != mu.size()){ throw teqp::InvalidArgument("bad size of mu"); @@ -685,13 +672,8 @@ class MultipolarContributionGrayGubbins { /// The variant containing the multipolar types that can be provided using multipolar_contributions_variant = std::variant< teqp::saft::polar_terms::GrossVrabec::MultipolarContributionGrossVrabec, - MultipolarContributionGrayGubbins, - MultipolarContributionGrayGubbins, - MultipolarContributionGrayGubbins, - - MultipolarContributionGubbinsTwu, - MultipolarContributionGubbinsTwu, - MultipolarContributionGubbinsTwu + MultipolarContributionGrayGubbins, + MultipolarContributionGubbinsTwu >; } diff --git a/include/teqp/models/saftvrmie.hpp b/include/teqp/models/saftvrmie.hpp index 2aff17b..fd3ecdb 100644 --- a/include/teqp/models/saftvrmie.hpp +++ b/include/teqp/models/saftvrmie.hpp @@ -1282,30 +1282,26 @@ inline auto SAFTVRMiefactory(const nlohmann::json & spec){ return SAFTVRMieMixture(std::move(chain), coeffs, std::move(polar)); } if (polar_model == "GubbinsTwu+Luckas"){ - using MCGTL = MultipolarContributionGubbinsTwu; auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval(); auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval(); - auto polar = MCGTL(sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::use_packing_fraction); + auto polar = MultipolarContributionGubbinsTwu(JLuckasSidecar{6}, KLuckasSidecar{222,333}, sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::use_packing_fraction); return SAFTVRMieMixture(std::move(chain), coeffs, std::move(polar)); } if (polar_model == "GubbinsTwu+GubbinsTwu"){ - using MCGG = MultipolarContributionGubbinsTwu; auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval(); auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval(); - auto polar = MCGG(sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::use_packing_fraction); + auto polar = MultipolarContributionGubbinsTwu(JGubbinsTwuSidecar{6}, KGubbinsTwuSidecar{222,333}, sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::use_packing_fraction); return SAFTVRMieMixture(std::move(chain), coeffs, std::move(polar)); } if (polar_model == "GubbinsTwu+Gottschalk"){ - using MCGG = MultipolarContributionGubbinsTwu; auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval(); auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval(); - auto polar = MCGG(sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::use_packing_fraction); + auto polar = MultipolarContributionGubbinsTwu(JGottschalkSidecar{6}, KGottschalkSidecar{222,333}, sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::use_packing_fraction); return SAFTVRMieMixture(std::move(chain), coeffs, std::move(polar)); } if (polar_model == "GrayGubbins+GubbinsTwu"){ - using MCGG = MultipolarContributionGrayGubbins; - auto polar = MCGG(sigma_ms, epsks, SIGMAIJ, EPSKIJ, mu_Cm, Q_Cm2, polar_flags); + auto polar = MultipolarContributionGrayGubbins(JGubbinsTwuSidecar{6}, KGubbinsTwuSidecar{222,333}, sigma_ms, epsks, SIGMAIJ, EPSKIJ, mu_Cm, Q_Cm2, polar_flags); return SAFTVRMieMixture(std::move(chain), coeffs, std::move(polar)); } // if (polar_model == "GrayGubbins+Gottschalk"){ @@ -1314,24 +1310,23 @@ inline auto SAFTVRMiefactory(const nlohmann::json & spec){ // return SAFTVRMieMixture(std::move(chain), std::move(polar)); // } if (polar_model == "GrayGubbins+Luckas"){ - using MCGG = MultipolarContributionGrayGubbins; - auto polar = MCGG(sigma_ms, epsks, SIGMAIJ, EPSKIJ, mu_Cm, Q_Cm2, polar_flags); + auto polar = MultipolarContributionGrayGubbins(JLuckasSidecar{6}, KLuckasSidecar{222,333}, sigma_ms, epsks, SIGMAIJ, EPSKIJ, mu_Cm, Q_Cm2, polar_flags); return SAFTVRMieMixture(std::move(chain), coeffs, std::move(polar)); } if (polar_model == "GubbinsTwu+Luckas+GubbinsTwuRhostar"){ - using MCGTL = MultipolarContributionGubbinsTwu; + using MCGTL = MultipolarContributionGubbinsTwu; auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval(); auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval(); - auto polar = MCGTL(sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::calculate_Gubbins_rhostar); + auto polar = MCGTL(JLuckasSidecar{6}, KLuckasSidecar{222,333}, sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::calculate_Gubbins_rhostar); return SAFTVRMieMixture(std::move(chain), coeffs, std::move(polar)); } if (polar_model == "GubbinsTwu+GubbinsTwu+GubbinsTwuRhostar"){ - using MCGG = MultipolarContributionGubbinsTwu; + using MCGG = MultipolarContributionGubbinsTwu; auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval(); auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval(); - auto polar = MCGG(sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::calculate_Gubbins_rhostar); + auto polar = MCGG(JGubbinsTwuSidecar{6}, KGubbinsTwuSidecar{222,333}, sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::calculate_Gubbins_rhostar); return SAFTVRMieMixture(std::move(chain), coeffs, std::move(polar)); } throw teqp::InvalidArgument("didn't understand this polar_model:"+polar_model); diff --git a/src/tests/catch_test_Grayiteration.cxx b/src/tests/catch_test_Grayiteration.cxx index 5547446..45d064d 100644 --- a/src/tests/catch_test_Grayiteration.cxx +++ b/src/tests/catch_test_Grayiteration.cxx @@ -212,7 +212,7 @@ TEST_CASE("Test calculation of polarized dipole moment with polarizable Lennard- } )"_json; flags["polarizable"]["alpha_symm / m^3"][0] = alpha_symm; - MultipolarContributionGrayGubbins GG{sigma_m_, epsilon_over_kB_, SIGMAIJ, EPSKBIJ, mu_, Q_, flags}; + MultipolarContributionGrayGubbins GG{JGubbinsTwuSidecar{6}, KGubbinsTwuSidecar{222,333}, sigma_m_, epsilon_over_kB_, SIGMAIJ, EPSKBIJ, mu_, Q_, flags}; double muprime = GG.iterate_muprime_SS(T, rhoN, rhostar, molefracs_, mu_, 20)[0]; double mustareffcalc = muprime/(sqrt(epsilon_over_kB*k_B*pow(sigma_m,3))/sqrt(k_e)); CAPTURE(ref); diff --git a/src/tests/catch_test_SAFTpolar.cxx b/src/tests/catch_test_SAFTpolar.cxx index ca7ad09..ded885e 100644 --- a/src/tests/catch_test_SAFTpolar.cxx +++ b/src/tests/catch_test_SAFTpolar.cxx @@ -114,8 +114,8 @@ TEST_CASE("Evaluate higher derivatives of K", "[GTK]") } -using MCGTL = MultipolarContributionGubbinsTwu; -using MCGG = MultipolarContributionGubbinsTwu; +using MCGTL = MultipolarContributionGubbinsTwu; +using MCGG = MultipolarContributionGubbinsTwu; TEST_CASE("Evaluation of Gubbins and Twu combos ", "[GTLPolar]") { @@ -124,13 +124,13 @@ TEST_CASE("Evaluation of Gubbins and Twu combos ", "[GTLPolar]") auto mubar2 = (Eigen::ArrayXd(2) << 0.0, 0.5).finished(); auto Qbar2 = (Eigen::ArrayXd(2) << 0.5, 0).finished(); SECTION("+ Luckas"){ - MCGTL GTL{sigma_m, epsilon_over_k, mubar2, Qbar2, multipolar_rhostar_approach::calculate_Gubbins_rhostar}; + MCGTL GTL{JLuckasSidecar{6}, KLuckasSidecar{222,333}, sigma_m, epsilon_over_k, mubar2, Qbar2, multipolar_rhostar_approach::calculate_Gubbins_rhostar}; auto z = (Eigen::ArrayXd(2) << 0.1, 0.9).finished(); auto rhoN = std::complex(300, 1e-100); GTL.eval(300.0, rhoN, rhoN, z); } SECTION("+ Gubbins&Twu"){ - MCGG GTL{sigma_m, epsilon_over_k, mubar2, Qbar2, multipolar_rhostar_approach::calculate_Gubbins_rhostar}; + MCGG GTL{JGubbinsTwuSidecar{6}, KGubbinsTwuSidecar{222,333}, sigma_m, epsilon_over_k, mubar2, Qbar2, multipolar_rhostar_approach::calculate_Gubbins_rhostar}; auto z = (Eigen::ArrayXd(2) << 0.1, 0.9).finished(); auto rhoN = std::complex(300, 1e-100); GTL.eval(300.0, rhoN, rhoN, z);