Fix return of unevaluated dual expression for CPA

Closes #82

include/teqp/math/finite_derivs.hpp

@@ -1,20 +1,20 @@
 #pragma once
 
+namespace teqp{
+
 /**
 * Routines for finite differentiation, useful for testing derivatives obtained by other methods
-* 
+*
 * From:
 * Bengt Fornberg, 1988, "Coefficients from Generation of Finite Difference Formulas on Arbitrarily Spaced Grids", MATHEMATICS OF COMPUTATION, v. 51, n. 184, pp. 699-706
-* 
+*
 * Higher derivatives should always be done in extended precision mathematics!
-* 
+*
 * Warning: these routines may give entirely erroneous results for double precision arithmetic,
 * especially the higher derivatives
-* 
+*
 * Warning: these routines are optimized for accuracy, not for speed or memory use
 */
-namespace teqp{
-
 template<int Nderiv, int Norder, typename Function, typename Scalar>
 auto centered_diff(const Function &f, const Scalar x, const Scalar h) {
 
@@ -63,4 +63,10 @@ auto centered_diff(const Function &f, const Scalar x, const Scalar h) {
 	auto val = num / pow(h, Nderiv);
 	return val;
 }
-}; // namespace teqp
+
+template<typename Function, typename Scalarx, typename Scalary>
+auto centered_diff_xy(const Function &f, const Scalarx x, const Scalary y, const Scalarx dx, const Scalary dy) {
+    return (f(x+dx, y+dy) - f(x+dx, y-dy) - f(x-dx, y+dy) + f(x-dx, y-dy))/(4*dx*dy);
+}
+
+}; // namespace teqp

include/teqp/models/CPA.hpp

@@ -159,7 +159,7 @@ class CPACubic {
 
     template<typename TType>
     auto get_ai(TType T, int i) const {
-        return a0[i] * POW2(1.0 + c1[i]*(1.0 - sqrt(T / Tc[i])));
+        return forceeval(a0[i] * POW2(1.0 + c1[i]*(1.0 - sqrt(T / Tc[i]))));
     }
 
     template<typename TType, typename VecType>

src/tests/catch_tests.cxx

@@ -430,11 +430,19 @@ TEST_CASE("Test water Clapeyron.jl", "[CPA]") {
     auto R = 8.31446261815324;
     CPA::CPACubic cub(CPA::cubic_flag::SRK, a0, bi, c1, Tc, R);
     double T = 303.15, rhomolar = 1/1.7915123921401366e-5;
+    using my_float_type = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<100U>>;
+    my_float_type Trecip = 1/T, rhof = rhomolar;
 
     auto z = (Eigen::ArrayXd(1) << 1).finished();
 
     using tdx = TDXDerivatives<decltype(cub)>;
     auto alphar = cub.alphar(T, rhomolar, molefrac);
+    auto Ar11_cub = tdx::get_Ar11(cub, T, rhomolar, z);
+    auto fcub = [&cub, &z](const auto& Trecip, const auto& rho) -> my_float_type{
+        return cub.alphar(forceeval(1/Trecip), rho, z);
+    };
+    auto Ar11_cd_cubic = static_cast<double>(centered_diff_xy(fcub, Trecip, rhof, my_float_type(1e-30), my_float_type(1e-30)) * Trecip * rhof);
+    CHECK(Ar11_cub == Approx(Ar11_cd_cubic));
     CHECK(alphar == Approx(-1.2713135319123854)); // Value from Clapeyron.jl
     double p_noassoc = T*rhomolar*R*(1+tdx::get_Ar01(cub, T, rhomolar, z));
     CAPTURE(p_noassoc);
@@ -451,6 +459,15 @@ TEST_CASE("Test water Clapeyron.jl", "[CPA]") {
     CAPTURE(p_withassoc);
 
     REQUIRE(p_withassoc == Approx(100000.000));
+
+    auto f = [&cpa, &z](const auto& Trecip, const auto& rho) -> my_float_type{
+        return cpa.alphar(forceeval(1/Trecip), rho, z);
+    };
+    auto Ar11_cd = static_cast<double>(centered_diff_xy(f, Trecip, rhof, my_float_type(1e-30), my_float_type(1e-30)) * Trecip * rhof);
+
+    auto Ar11 = tdc::get_Ar11(cpa, T, rhomolar, z);
+    REQUIRE(std::isfinite(Ar11));
+    CHECK(Ar11 == Approx(Ar11_cd));
 }
 
 TEST_CASE("Test water", "[CPA]") {
@@ -481,6 +498,31 @@ TEST_CASE("Test water", "[CPA]") {
     REQUIRE(p_withassoc == Approx(312682.0709));
 }
 
+TEST_CASE("Test [C2mim][NTF2]", "[CPA]") {
+    std::valarray<double> a0 = {25.8}, bi = {251.70e-6}, c1 = {0.273}, Tc = {1244.9},
+                          molefrac = {1.0};
+    auto R = 8.3144598;
+    CPA::CPACubic cub(CPA::cubic_flag::SRK, a0, bi, c1, Tc, R);
+    double T = 313, rhomolar = 1503/0.39132; // From https://pubs.acs.org/doi/epdf/10.1021/je700205n
+
+    auto z = (Eigen::ArrayXd(1) << 1).finished();
+
+    using tdx = TDXDerivatives<decltype(cub)>;
+    auto alphar = cub.alphar(T, rhomolar, molefrac);
+
+    std::vector<CPA::association_classes> schemes = { CPA::association_classes::a2B };
+    std::valarray<double> epsAB = { 177388/10.0 }, betaAB = { 0.00266 };
+    CPA::CPAAssociation cpaa(std::move(cub), schemes, CPA::radial_dist::KG, epsAB, betaAB, R);
+
+    CPA::CPAEOS cpa(std::move(cub), std::move(cpaa));
+    using tdc = TDXDerivatives<decltype(cpa)>;
+    auto Ar01 = tdc::get_Ar01(cpa, T, rhomolar, z);
+    double p_withassoc = T*rhomolar*R*(1 + Ar01);
+    CAPTURE(p_withassoc);
+
+//    REQUIRE(p_withassoc == Approx(312682.0709));
+}
+
 TEST_CASE("Test water w/ factory", "") {
    using namespace CPA;
    nlohmann::json water = {