Implement the Luckas J and K integrals

include/teqp/models/saft/correlation_integrals.hpp

+#pragma once
+
+#include <map>
+#include <array>
+
+namespace teqp {
+namespace SAFTpolar{
+
+static const std::map<int, std::array<double, 12>> Luckas_J_coeffs = {
+    {4,  {-1.38410152e00,  -7.05792933e-01, 2.60947023e00,  1.96828333e01, 1.13619510e01, -2.98510490e01,  -3.15686398e01, -2.00943290e01,  5.11029320e01, 1.44194150e01, 9.40061069e00,  -2.36844608e01}},
+    {5,  {-6.89702637e-01, -1.62382602e-01, 1.16302441e00,  1.42067443e01, 4.59642681e00, -1.81421003e01,  -2.45012804e01, -8.42839734e00,  3.25579587e01, 1.16339969e01, 4.00080085e00,  -1.54419815e01}},
+    {6,  {-4.57433648e-01, -4.74592642e-02, 7.36156406e-01, 1.19830560e01, 2.52201906e00, -1.40252699e01,  -2.15154325e01, -4.69572629e00,  2.59035512e01, 1.04390419e01, 2.24160266e00,  -1.24580635e01}},
+    {7,  {-3.46551277e-01, -8.67899015e-03, 5.41729996e-01, 1.07102246e01, 1.58762511e00, -1.19085356e01,  -1.98486869e01, -2.98067526e00,  2.25425086e01, 9.79650295e00, 1.42953263e00,  -1.09749041e01}},
+    {8,  {-2.85932401e-01,  7.44866141e-03, 4.35471488e-01, 9.87890138e00, 1.08489226e00, -1.06280961e01,  -1.88507949e01, -2.04803671e00,  2.06131402e01, 9.44834359e00, 9.88551760e-01, -1.01576036e01}},
+    {9,  {-2.50931500e-01,  1.49604396e-02, 3.71768229e-01, 9.29877256e00, 7.84100755e-01, -9.78235175e00, -1.82609129e01, -1.48667167e00,  1.94537822e01, 9.28527644e00, 7.24545266e-01, -9.70374243e00}},
+    {10, {-2.30585563e-01,  1.86890174e-02, 3.31660880e-01, 8.88275358e00, 5.90996555e-01, -9.19786291e00, -1.79530632e01, -1.12502567e00,  1.87717642e01, 9.25235661e00, 5.55909350e-01, -9.47729033e00}},
+    {11, {-2.19221482e-01,  2.05982865e-02, 3.05823216e-01, 8.58399185e00, 4.60555635e-01, -8.78667564e00, -1.78552884e01, -8.80186074e-01, 1.84146650e01, 9.31691675e00, 4.43027169e-01, -9.40608814e00}},
+    {12, {-2.13591301e-01,  2.15809377e-02, 2.89113218e-01, 8.37291362e00, 3.68864477e-01, -8.49747359e00, -1.79188678e01, -7.07751483e-01, 1.82906357e01, 9.45598216e00, 3.64636215e-01, -9.44579903e00}},
+    {13, {-2.11696363e-01,  2.20767096e-02, 2.78453228e-01, 8.23384704e00, 3.02478832e-01, -8.30238277e00,  1.81183432e01, -5.82685843e-01, 1.83497132e01, 9.65764061e00, 3.08765033e-01, -9.57226243e00}},
+    {14, {-2.12238116e-01,  2.23157164e-02, 2.71883622e-01, 8.15208365e00, 2.53073339e-01, -8.17919021e00, -1.84273649e01, -4.89407317e-01, 1.85502740e01, 9.90906951e00, 2.67965960e-01, -9.76495200e00}},
+    {15, {-2.14335965e-01,  2.24240261e-02, 2.68094773e-01, 8.11899188e00, 2.15506735e-01, -8.11465705e00, -1.88310645e01, -4.18309476e-01, 1.88679367e01, 1.02033085e01, 2.37674032e-01, -1.00120648e01}},
+};
+
+template<int n>
+class LuckasJIntegral{
+public:
+    
+    const std::array<double, 12> a;
+    double a00,a01,a02,a10,a11,a12,a20,a21,a22,a30,a31,a32;
+    
+    LuckasJIntegral() : a(Luckas_J_coeffs.at(n)){
+        a00 = a[0]; a01 = a[1]; a02 = a[2];
+        a10 = a[3]; a11 = a[4]; a12 = a[5];
+        a20 = a[6]; a21 = a[7]; a22 = a[8];
+        a30 = a[9]; a31 = a[10]; a32 = a[11];
+    }
+    
+    template<typename TType, typename RhoType>
+    auto get_J(const TType& Tstar, const RhoType& rhostar){
+        auto Z_1 = 0.3 + 0.05*n;
+        auto Z_2 = 1.0/n;
+        auto A_0 = a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar;
+        auto A_1 = a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar;
+        auto A_2 = a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar;
+        std::common_type_t<TType, RhoType> out = (A_0 + A_1*pow(Tstar, Z_1) + A_2*pow(Tstar, Z_2))*exp(1.0/(Tstar + 4.0/pow(std::abs(log(rhostar/sqrt(2.0))), 3.0)));
+        return out;
+    }
+};
+
+static const std::map<std::tuple<int, int>, std::array<double, 16>> Luckas_K_coeffs = {
+    {{222,333},{  1.23746712e-02,  1.90127007e-04,  5.92861480e-03,  1.88387469e-04,  1.51499090e-01, -4.03661304e-03, -1.71653218e-01,  7.44958473e-02,  4.61939219e-01,  7.61527540e-03, -1.30641988e00,   7.63483402e-01,  7.36061418e-03, -3.14115085e-03, -2.62332295e-01,  2.99061150e-01}},
+    {{233,344},{  3.51281055e-03,  2.12654833e-04,  3.26699985e-03,  1.71420492e-04,  9.78591707e-02, -3.24426472e-03, -9.85906322e-02,  4.31053752e-02,  2.68515490e-01,  6.25909413e-03, -7.46179517e-01,  4.10644114e-01,  1.28223189e-01, -2.54348393e-03, -5.40215100e-01,  4.50896967e-01}},
+    {{334,445},{ -1.77119569e-03, -2.60437070e-04, -2.65044611e-03, -1.96937245e-04, -8.97801454e-02,  2.93120169e-03,  8.77755068e-02, -3.76799716e-02, -2.19892289e-01, -5.73110168e-03,  6.37326406e-01, -3.52780095e-01, -1.57228990e-01,  2.04489194e-03,  5.59739399e-01, -4.42956307e-01}},
+    {{444,555},{  1.07862851e-03,  2.88647260e-04,  2.07369614e-03,  1.79600129e-04,  7.57084957e-02, -2.28761197e-03, -7.33670248e-02,  3.12736172e-02,  1.73518964e-01,  4.54174801e-03, -5.14488122e-01,  2.85483011e-01,  1.71670939e-01, -1.07503182e-03, -5.51880129e-01,  4.16887229e-01}}
+};
+
+template<int n1, int n2>
+class LuckasKIntegral{
+public:
+    
+    const std::array<double, 16> a;
+    double a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31, a32, a33;
+    
+    LuckasKIntegral() : 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];
+    }
+    
+    template<typename TType, typename RhoType>
+    auto get_K(const TType& Tstar, const RhoType& rhostar){
+        auto Z_1 = 2.0;
+        auto Z_2 = 3.0;
+        auto Z_3 = 4.0;
+        auto b_0 = a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar;
+        auto b_1 = a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar;
+        auto b_2 = a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar;
+        auto b_3 = a03 + a13*rhostar + a23*rhostar*rhostar + a33*rhostar*rhostar*rhostar;
+        std::common_type_t<TType, RhoType> 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;
+    }
+};
+
+}
+}

src/tests/catch_test_SAFTpolar.cxx

+#include <catch2/catch_test_macros.hpp>
+#include <catch2/catch_approx.hpp>
+
+using Catch::Approx;
+
+#include "teqp/models/saft/correlation_integrals.hpp"
+
+using namespace teqp::SAFTpolar;
+
+TEST_CASE("Evaluation of J^{(n)}", "[LuckasJn]")
+{
+    LuckasJIntegral<12> J12;
+    auto Jval = J12.get_J(3.0, 1.0);
+}
+
+TEST_CASE("Evaluation of K(xxx, yyy)", "[LuckasKnn]")
+{
+    auto Kval23 = LuckasKIntegral<222, 333>().get_K(1.0, 0.9);
+    CHECK(Kval23 == Approx(0.03332).margin(0.02));
+    
+    auto Kval45 = LuckasKIntegral<444, 555>().get_K(1.0, 0.9);
+    CHECK(Kval45 == Approx(0.01541).margin(0.02));
+}