Add method for residual entropy

include/teqp/derivs.hpp

@@ -38,7 +38,7 @@ get_Psir(const Model& model, const TType T, const ContainerType& rhovec) {
 }
 
 /**
-// Calculate the residual pressure from derivatives of alphar
+/// Calculate the residual pressure from derivatives of alphar
 */
 template <typename Model, typename TType, typename ContainerType>
 typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
@@ -52,6 +52,15 @@ get_pr(const Model& model, const TType T, const ContainerType& rhovec) {
     return pr*rhotot_*model.R*T;
 }
 
+/**
+/// Calculate the residual entropy (s^+=-sr/R) from derivatives of alphar
+*/
+template <typename Model, typename TType, typename ContainerType>
+typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
+get_splus(const Model& model, const TType T, const ContainerType& rhovec) {
+    return model.alphar(T, rhovec) + T*derivT([&model](const auto& T, const auto& rhovec) { return model.alphar(T, rhovec); }, T, rhovec);
+}
+
 template<typename Model, typename TType, typename RhoType>
 auto build_Psir_Hessian(const Model& model, const TType T, const RhoType& rho) {
     // Double derivatives in each component's concentration

src/main.cpp

@@ -60,6 +60,7 @@ void test_vdwMix() {
     auto dPsirdrho0 = rhovec[0]*derivrhoi(fPsir, T, rhovec, 0);
     auto dPsirdrho1 = rhovec[1]*derivrhoi(fPsir, T, rhovec, 1);
     auto pfromderiv = rho*R*T - Psir + dPsirdrho0 + dPsirdrho1;
+    auto sr = get_splus(vdW, T, rhovec);
 
     auto t2 = std::chrono::steady_clock::now();
     auto pfromderiv3 = rhotot*R*T + get_pr(vdW, T, rhovec);