From a02a614948401ef6ea3bdbbd2e984e2caf9bf4ca Mon Sep 17 00:00:00 2001
From: Trine Mykkeltvedt <tmyk@norceresearch.no>
Date: Thu, 23 Jun 2022 14:47:50 +0200
Subject: [PATCH] more cleaning in the PengRobinson files

---
 opm/material/eos/PengRobinson.hpp        |  2 +-
 opm/material/eos/PengRobinsonMixture.hpp | 28 ------------------------
 2 files changed, 1 insertion(+), 29 deletions(-)

diff --git a/opm/material/eos/PengRobinson.hpp b/opm/material/eos/PengRobinson.hpp
index 3c57743ec..8fedd2f4e 100644
--- a/opm/material/eos/PengRobinson.hpp
+++ b/opm/material/eos/PengRobinson.hpp
@@ -186,7 +186,7 @@ public:
         Valgrind::CheckDefined(a2);
         Valgrind::CheckDefined(a3);
         Valgrind::CheckDefined(a4);
-       // int numSol = invertCubicPolynomial(Z, a1, a2, a3, a4);
+
         int numSol = cubicRoots(Z, a1, a2, a3, a4);
         if (numSol == 3) {
             // the EOS has three intersections with the pressure,
diff --git a/opm/material/eos/PengRobinsonMixture.hpp b/opm/material/eos/PengRobinsonMixture.hpp
index 94cd8e724..eb9844c26 100644
--- a/opm/material/eos/PengRobinsonMixture.hpp
+++ b/opm/material/eos/PengRobinsonMixture.hpp
@@ -122,26 +122,8 @@ public:
         LhsEval Astar = params.a(phaseIdx)*p/(RT*RT);
         LhsEval Bstar = params.b(phaseIdx)*p/(RT);
 
-        //MAYBE REMOVE THIS COMPLETELY ??
-        // calculate delta_i (see: Reid, p. 145)
-        LhsEval sumMoleFractions = 0.0;
-        for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
-            sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx);
-        }
-        LhsEval deltai = 2*sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx);
-        LhsEval tmp = 0;
-        for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
-           tmp +=
-                fs.moleFraction(phaseIdx, compJIdx)
-                / sumMoleFractions
-                * sqrt(params.aPure(phaseIdx, compJIdx));
-         };
-         deltai *= tmp;
-
-        // Calculate A_s LIKE IN JULIA CODE
         LhsEval A_s = 0.0;
         for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
-            //param
             A_s += params.aCache(phaseIdx, compIdx, compJIdx) * fs.moleFraction(phaseIdx, compJIdx) * p / (RT * RT);
         }
 
@@ -157,16 +139,6 @@ public:
         m1 = 0.5*(u + std::sqrt(u*u - 4*w));
         m2 = 0.5*(u - std::sqrt(u*u - 4*w));
 
-
-        LhsEval base =
-            (Z + Bstar*m1) /
-            (Z + Bstar*m2);
-        LhsEval expo =  Astar/(Bstar*std::sqrt(u*u - 4*w))*(bi_b - deltai);
-
-        //α = -log(Z - B) + (B_i[i]/B)*(Z - 1)
-        //β = log((Z + m2*B)/(Z + m1*B))*A/(Δm*B)
-        //γ = (2/A)*A_s - B_i[i]/B
-        //ln_phi = α + β*γ
         alpha = -log(Z - Bstar) + bi_b * (Z - 1);
         betta = log((Z + m2 * Bstar) / (Z + m1 * Bstar)) * Astar / ((m1 - m2) * Bstar);
         gamma = (2 / Astar ) * A_s - bi_b;