From 94bda214209e5aa7a1e97ff2894cd357d6897f8a Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.h.bell@gmail.com>
Date: Wed, 3 Jun 2015 21:32:20 -0600
Subject: [PATCH] Analytic derivatives for M1_star implemented

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
---
 src/Backends/Helmholtz/MixtureDerivatives.h | 44 +++++++++++++++++----
 1 file changed, 37 insertions(+), 7 deletions(-)

diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h
index 4d4adb44..1b61fbab 100644
--- a/src/Backends/Helmholtz/MixtureDerivatives.h
+++ b/src/Backends/Helmholtz/MixtureDerivatives.h
@@ -14,6 +14,7 @@
 #ifndef MIXTURE_DERIVATIVES_H
 #define MIXTURE_DERIVATIVES_H
 
+#include <Eigen/Core>
 #include "HelmholtzEOSMixtureBackend.h"
 
 namespace CoolProp{
@@ -177,18 +178,47 @@ class MixtureDerivatives{
         // Fill in the symmetric elements
         for (std::size_t i = 0; i < N; ++i){
             for (std::size_t j = 0; j < i; ++j){
-                L(i, j) = nAij(HEOS, i, j, xN_flag);
+                L(i, j) = L(j, i);
             }
         }
         return L.determinant();
     }
     static CoolPropDbl M1_star(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){
-        Eigen::Matrix2d M1;
-        M1(0, 0) = nAij(HEOS, 0, 0, xN_flag);
-        M1(0, 1) = nAij(HEOS, 0, 1, xN_flag);
-        M1(1, 0) = nAij(HEOS, 0, 0, xN_flag)*n2Aijk(HEOS, 0, 1, 1, xN_flag) + nAij(HEOS, 1, 1, xN_flag)*n2Aijk(HEOS, 0, 0, 0, xN_flag) - 2*nAij(HEOS, 0, 1, xN_flag)*n2Aijk(HEOS, 0, 0, 1, xN_flag);;
-        M1(1, 1) = nAij(HEOS, 0, 0, xN_flag)*n2Aijk(HEOS, 1, 1, 1, xN_flag) + nAij(HEOS, 1, 1, xN_flag)*n2Aijk(HEOS, 0, 0, 1, xN_flag) - 2*nAij(HEOS, 0, 1, xN_flag)*n2Aijk(HEOS, 0, 1, 1, xN_flag);
-        return M1.determinant();
+
+        std::size_t N = HEOS.mole_fractions.size();
+        Eigen::MatrixXd L;
+        L.resize(N, N);
+        for (std::size_t i = 0; i < N; ++i){
+            for (std::size_t j = i; j < N; ++j){
+                L(i, j) = nAij(HEOS, i, j, xN_flag);
+            }
+        }
+        // Fill in the symmetric elements
+        for (std::size_t i = 0; i < N; ++i){
+            for (std::size_t j = 0; j < i; ++j){
+                L(i, j) = L(j, i);
+            }
+        }
+
+        Eigen::MatrixXd M = L;
+
+        // Last row
+        for (std::size_t i = 0; i < N; ++i){
+            Eigen::MatrixXd n2dLdni(N, N);
+            for (std::size_t j = 0; j < N; ++j){
+                for (std::size_t k = j; k < N; ++k){
+                    n2dLdni(j, k) = n2Aijk(HEOS, j, k, i, xN_flag);
+                }
+            }
+            // Fill in the symmetric elements
+            for (std::size_t j = 0; j < N; ++j){
+                for (std::size_t k = 0; k < j; ++k){
+                    n2dLdni(j, k) = n2dLdni(k, j);
+                }
+            }
+            M(N-1, i) = (adjugate(L, N)*n2dLdni).trace();
+        }
+        return M.determinant();
     }
     
     /** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions