From 3e5362de33603b9ea9283f59d300f80ed262e8fc Mon Sep 17 00:00:00 2001
From: Even Solbraa <41290109+EvenSol@users.noreply.github.com>
Date: Tue, 30 Apr 2024 22:39:14 +0000
Subject: [PATCH 1/4] initial test rachford rice

---
 .../flashOps/RachfordRice.java                | 205 ++++++++++++++++++
 .../flashOps/TPflash.java                     |   4 +-
 2 files changed, 208 insertions(+), 1 deletion(-)
 create mode 100644 src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java

diff --git a/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java b/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java
new file mode 100644
index 000000000..72838ee25
--- /dev/null
+++ b/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java
@@ -0,0 +1,205 @@
+/*
+ * RachfordRice.java
+ *
+ * Created on 1.april 2024
+ */
+
+package neqsim.thermodynamicOperations.flashOps;
+
+import neqsim.thermo.component.ComponentInterface;
+import neqsim.thermo.system.SystemInterface;
+
+/**
+ * RachfordRice classes.
+ *
+ * @author Even Solbraa
+ */
+public class RachfordRice {
+  private static final long serialVersionUID = 1000;
+  double[] beta = new double[2];
+  SystemInterface fluid = null;
+
+  public RachfordRice(SystemInterface fluid) {
+    this.fluid = fluid;
+  }
+
+  public final double calcBeta() throws neqsim.util.exception.IsNaNException,
+      neqsim.util.exception.TooManyIterationsException {
+    ComponentInterface[] compArray = fluid.getPhase(0).getComponents();
+
+    int i;
+    double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
+    double deriv = 0.0;
+    double gbeta = 0.0;
+    double betal = 0;
+    double nybeta = 0;
+
+    double midler = 0;
+    double minBeta = tolerance;
+    double maxBeta = 1.0 - tolerance;
+    double g0 = -1.0;
+    double g1 = 1.0;
+
+
+    beta[0] = 0.5;
+    beta[1] = 0.5;
+    nybeta = fluid.getBeta(0);
+    betal = 1.0 - nybeta;
+
+    for (i = 0; i < fluid.getNumberOfComponents(); i++) {
+      midler = (compArray[i].getK() * compArray[i].getz() - 1.0) / (compArray[i].getK() - 1.0);
+      if ((midler > minBeta) && (compArray[i].getK() > 1.0)) {
+        minBeta = midler;
+      }
+      midler = (1.0 - compArray[i].getz()) / (1.0 - compArray[i].getK());
+      if ((midler < maxBeta) && (compArray[i].getK() < 1.0)) {
+        maxBeta = midler;
+      }
+      g0 += compArray[i].getz() * compArray[i].getK();
+      g1 += -compArray[i].getz() / compArray[i].getK();
+    }
+
+    if (g0 < 0) {
+      this.beta[1] = 1.0 - tolerance;
+      this.beta[0] = tolerance;
+      return this.beta[0];
+    }
+    if (g1 > 0) {
+      this.beta[1] = tolerance;
+      this.beta[0] = 1.0 - tolerance;
+      return this.beta[0];
+    }
+
+    nybeta = (minBeta + maxBeta) / 2.0;
+    // System.out.println("guessed beta: " + nybeta + " maxbeta: " +maxBeta + "
+    // minbeta: " +minBeta );
+    betal = 1.0 - nybeta;
+
+    // ' *l = 1.0-nybeta;
+    double gtest = 0.0;
+    for (i = 0; i < fluid.getNumberOfComponents(); i++) {
+      gtest += compArray[i].getz() * (compArray[i].getK() - 1.0)
+          / (1.0 - nybeta + nybeta * compArray[i].getK()); // beta
+                                                           // =
+                                                           // nybeta
+    }
+
+    if (gtest >= 0) {
+      minBeta = nybeta;
+    } else {
+      maxBeta = nybeta;
+    }
+
+    if (gtest < 0) {
+      double minold = minBeta;
+      minBeta = 1.0 - maxBeta;
+      maxBeta = 1.0 - minold;
+    }
+
+    int iterations = 0;
+    int maxIterations = 300;
+    // System.out.println("gtest: " + gtest);
+    double step = 1.0;
+    do {
+      iterations++;
+      if (gtest >= 0) {
+        // oldbeta = nybeta;
+        deriv = 0.0;
+        gbeta = 0.0;
+
+        for (i = 0; i < fluid.getNumberOfComponents(); i++) {
+          double temp1 = (compArray[i].getK() - 1.0);
+          double temp2 = 1.0 + temp1 * nybeta;
+          deriv += -(compArray[i].getz() * temp1 * temp1) / (temp2 * temp2);
+          gbeta += compArray[i].getz() * (compArray[i].getK() - 1.0)
+              / (1.0 + (compArray[i].getK() - 1.0) * nybeta);
+        }
+
+        if (gbeta >= 0) {
+          minBeta = nybeta;
+        } else {
+          maxBeta = nybeta;
+        }
+        nybeta -= (gbeta / deriv);
+
+        // System.out.println("beta: " + maxBeta);
+        if (nybeta > maxBeta) {
+          nybeta = maxBeta;
+        }
+        if (nybeta < minBeta) {
+          nybeta = minBeta;
+        }
+
+        /*
+         * if ((nybeta > maxBeta) || (nybeta < minBeta)) { // nybeta = 0.5 * (maxBeta + minBeta);
+         * gbeta = 1.0; }
+         */
+      } else {
+        // oldbeta = betal;
+        deriv = 0.0;
+        gbeta = 0.0;
+
+        for (i = 0; i < fluid.getNumberOfComponents(); i++) {
+          deriv -= (compArray[i].getz() * (compArray[i].getK() - 1.0) * (1.0 - compArray[i].getK()))
+              / Math.pow((betal + (1 - betal) * compArray[i].getK()), 2);
+          gbeta += compArray[i].getz() * (compArray[i].getK() - 1.0)
+              / (betal + (-betal + 1.0) * compArray[i].getK());
+        }
+
+        if (gbeta < 0) {
+          minBeta = betal;
+        } else {
+          maxBeta = betal;
+        }
+
+        betal -= (gbeta / deriv);
+
+        if (betal > maxBeta) {
+          betal = maxBeta;
+        }
+        if (betal < minBeta) {
+          betal = minBeta;
+        }
+
+        /*
+         * if ((betal > maxBeta) || (betal < minBeta)) { gbeta = 1.0; { betal = 0.5 * (maxBeta +
+         * minBeta); } }
+         */
+        nybeta = 1.0 - betal;
+      }
+      step = gbeta / deriv;
+      // System.out.println("step : " + step);
+    } while (Math.abs(step) >= 1.0e-10 && iterations < maxIterations); // &&
+    // (Math.abs(nybeta)-Math.abs(maxBeta))>0.1);
+
+    // System.out.println("beta: " + nybeta + " iterations: " + iterations);
+    if (nybeta <= tolerance) {
+      // this.phase = 1;
+      nybeta = tolerance;
+    } else if (nybeta >= 1.0 - tolerance) {
+      // this.phase = 0;
+      nybeta = 1.0 - tolerance;
+      // superheated vapour
+    } else {
+      // this.phase = 2;
+    } // two-phase liquid-gas
+
+    this.beta[0] = nybeta;
+    this.beta[1] = 1.0 - nybeta;
+    fluid.setBeta(fluid.getPhaseIndex(0), nybeta);
+    fluid.setBeta(fluid.getPhaseIndex(1), 1.0 - nybeta);
+
+    if (iterations >= maxIterations) {
+      throw new neqsim.util.exception.TooManyIterationsException(this, "calcBeta", maxIterations);
+    }
+    if (Double.isNaN(beta[1])) {
+      /*
+       * for (i = 0; i < numberOfComponents; i++) { System.out.println("K " + compArray[i].getK());
+       * System.out.println("z " + compArray[i].getz()); }
+       */
+      throw new neqsim.util.exception.IsNaNException(this, "calcBeta", "beta");
+    }
+    return this.beta[0];
+  }
+
+}
diff --git a/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java b/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java
index 78a856f41..3871ebfe8 100644
--- a/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java
+++ b/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java
@@ -91,7 +91,9 @@ public void sucsSubs() {
 
     double oldBeta = system.getBeta();
     try {
-      system.calcBeta();
+      // system.calcBeta();
+      RachfordRice racrice = new RachfordRice(system);
+      racrice.calcBeta();
     } catch (IsNaNException ex) {
       logger.warn("Not able to calculate beta. Value is NaN");
       system.setBeta(oldBeta);

From 22f52d3ade77099e648b57c001c56a0a86df3a29 Mon Sep 17 00:00:00 2001
From: Even Solbraa <41290109+EvenSol@users.noreply.github.com>
Date: Wed, 1 May 2024 18:05:40 +0000
Subject: [PATCH 2/4] update rachford rice - moved to separate class

---
 .../neqsim/thermo/system/SystemInterface.java |  18 +-
 .../neqsim/thermo/system/SystemThermo.java    | 175 ------------------
 .../flashOps/Flash.java                       |   2 +-
 .../flashOps/RachfordRice.java                | 100 ++++------
 .../flashOps/TPflash.java                     |  10 +-
 5 files changed, 40 insertions(+), 265 deletions(-)

diff --git a/src/main/java/neqsim/thermo/system/SystemInterface.java b/src/main/java/neqsim/thermo/system/SystemInterface.java
index ec518e153..d786f68d1 100644
--- a/src/main/java/neqsim/thermo/system/SystemInterface.java
+++ b/src/main/java/neqsim/thermo/system/SystemInterface.java
@@ -385,18 +385,6 @@ public void addTBPfraction(String componentName, double numberOfMoles, double mo
    */
   public void calc_x_y_nonorm();
 
-  /**
-   * <p>
-   * calcBeta. For simple gas liquid systems.
-   * </p>
-   *
-   * @return Beta Mole fraction contained in the heaviest phase, i.e., liquid phase.
-   * @throws neqsim.util.exception.IsNaNException if any.
-   * @throws neqsim.util.exception.TooManyIterationsException if any.
-   */
-  public double calcBeta()
-      throws neqsim.util.exception.IsNaNException, neqsim.util.exception.TooManyIterationsException;
-
   /**
    * <p>
    * calcHenrysConstant.
@@ -2431,14 +2419,14 @@ public void setImplementedTemperatureDeriativesofFugacity(
    */
   public void setPhaseType(int phaseToChange, PhaseType pt);
 
-    /**
+  /**
    * Change the phase type of a given phase.
    *
    * @param phaseToChange the phase number of the phase to set phase type
-   * @param phaseName     String to set
+   * @param phaseName String to set
    */
   public void setPhaseType(int phaseToChange, String phaseName);
-  
+
   /**
    * Set the physical property model type for each phase of the System.
    *
diff --git a/src/main/java/neqsim/thermo/system/SystemThermo.java b/src/main/java/neqsim/thermo/system/SystemThermo.java
index 2f252704c..7cfedd94b 100644
--- a/src/main/java/neqsim/thermo/system/SystemThermo.java
+++ b/src/main/java/neqsim/thermo/system/SystemThermo.java
@@ -1241,181 +1241,6 @@ public final void calc_x_y_nonorm() {
     }
   }
 
-  /** {@inheritDoc} */
-  @Override
-  public final double calcBeta() throws neqsim.util.exception.IsNaNException,
-      neqsim.util.exception.TooManyIterationsException {
-    ComponentInterface[] compArray = getPhase(0).getComponents();
-
-    int i;
-    double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
-    double deriv = 0.0;
-    double gbeta = 0.0;
-    double betal = 0;
-    double nybeta = 0;
-
-    double midler = 0;
-    double minBeta = tolerance;
-    double maxBeta = 1.0 - tolerance;
-    double g0 = -1.0;
-    double g1 = 1.0;
-    nybeta = this.beta[0];
-    betal = 1.0 - nybeta;
-
-    for (i = 0; i < numberOfComponents; i++) {
-      midler = (compArray[i].getK() * compArray[i].getz() - 1.0) / (compArray[i].getK() - 1.0);
-      if ((midler > minBeta) && (compArray[i].getK() > 1.0)) {
-        minBeta = midler;
-      }
-      midler = (1.0 - compArray[i].getz()) / (1.0 - compArray[i].getK());
-      if ((midler < maxBeta) && (compArray[i].getK() < 1.0)) {
-        maxBeta = midler;
-      }
-      g0 += compArray[i].getz() * compArray[i].getK();
-      g1 += -compArray[i].getz() / compArray[i].getK();
-    }
-
-    if (g0 < 0) {
-      this.beta[1] = 1.0 - tolerance;
-      this.beta[0] = tolerance;
-      return this.beta[0];
-    }
-    if (g1 > 0) {
-      this.beta[1] = tolerance;
-      this.beta[0] = 1.0 - tolerance;
-      return this.beta[0];
-    }
-
-    nybeta = (minBeta + maxBeta) / 2.0;
-    // System.out.println("guessed beta: " + nybeta + " maxbeta: " +maxBeta + "
-    // minbeta: " +minBeta );
-    betal = 1.0 - nybeta;
-
-    // ' *l = 1.0-nybeta;
-    double gtest = 0.0;
-    for (i = 0; i < numberOfComponents; i++) {
-      gtest += compArray[i].getz() * (compArray[i].getK() - 1.0)
-          / (1.0 - nybeta + nybeta * compArray[i].getK()); // beta
-                                                           // =
-                                                           // nybeta
-    }
-
-    if (gtest >= 0) {
-      minBeta = nybeta;
-    } else {
-      maxBeta = nybeta;
-    }
-
-    if (gtest < 0) {
-      double minold = minBeta;
-      minBeta = 1.0 - maxBeta;
-      maxBeta = 1.0 - minold;
-    }
-
-    int iterations = 0;
-    int maxIterations = 300;
-    // System.out.println("gtest: " + gtest);
-    double step = 1.0;
-    do {
-      iterations++;
-      if (gtest >= 0) {
-        // oldbeta = nybeta;
-        deriv = 0.0;
-        gbeta = 0.0;
-
-        for (i = 0; i < numberOfComponents; i++) {
-          double temp1 = (compArray[i].getK() - 1.0);
-          double temp2 = 1.0 + temp1 * nybeta;
-          deriv += -(compArray[i].getz() * temp1 * temp1) / (temp2 * temp2);
-          gbeta += compArray[i].getz() * (compArray[i].getK() - 1.0)
-              / (1.0 + (compArray[i].getK() - 1.0) * nybeta);
-        }
-
-        if (gbeta >= 0) {
-          minBeta = nybeta;
-        } else {
-          maxBeta = nybeta;
-        }
-        nybeta -= (gbeta / deriv);
-
-        // System.out.println("beta: " + maxBeta);
-        if (nybeta > maxBeta) {
-          nybeta = maxBeta;
-        }
-        if (nybeta < minBeta) {
-          nybeta = minBeta;
-        }
-
-        /*
-         * if ((nybeta > maxBeta) || (nybeta < minBeta)) { // nybeta = 0.5 * (maxBeta + minBeta);
-         * gbeta = 1.0; }
-         */
-      } else {
-        // oldbeta = betal;
-        deriv = 0.0;
-        gbeta = 0.0;
-
-        for (i = 0; i < numberOfComponents; i++) {
-          deriv -= (compArray[i].getz() * (compArray[i].getK() - 1.0) * (1.0 - compArray[i].getK()))
-              / Math.pow((betal + (1 - betal) * compArray[i].getK()), 2);
-          gbeta += compArray[i].getz() * (compArray[i].getK() - 1.0)
-              / (betal + (-betal + 1.0) * compArray[i].getK());
-        }
-
-        if (gbeta < 0) {
-          minBeta = betal;
-        } else {
-          maxBeta = betal;
-        }
-
-        betal -= (gbeta / deriv);
-
-        if (betal > maxBeta) {
-          betal = maxBeta;
-        }
-        if (betal < minBeta) {
-          betal = minBeta;
-        }
-
-        /*
-         * if ((betal > maxBeta) || (betal < minBeta)) { gbeta = 1.0; { betal = 0.5 * (maxBeta +
-         * minBeta); } }
-         */
-        nybeta = 1.0 - betal;
-      }
-      step = gbeta / deriv;
-      // System.out.println("step : " + step);
-    } while (Math.abs(step) >= 1.0e-10 && iterations < maxIterations); // &&
-    // (Math.abs(nybeta)-Math.abs(maxBeta))>0.1);
-
-    // System.out.println("beta: " + nybeta + " iterations: " + iterations);
-    if (nybeta <= tolerance) {
-      // this.phase = 1;
-      nybeta = tolerance;
-    } else if (nybeta >= 1.0 - tolerance) {
-      // this.phase = 0;
-      nybeta = 1.0 - tolerance;
-      // superheated vapour
-    } else {
-      // this.phase = 2;
-    } // two-phase liquid-gas
-
-    this.beta[0] = nybeta;
-    this.beta[1] = 1.0 - nybeta;
-
-    if (iterations >= maxIterations) {
-      throw new neqsim.util.exception.TooManyIterationsException(this, "calcBeta", maxIterations);
-    }
-    if (Double.isNaN(beta[1])) {
-      /*
-       * for (i = 0; i < numberOfComponents; i++) { System.out.println("K " + compArray[i].getK());
-       * System.out.println("z " + compArray[i].getz()); }
-       */
-      throw new neqsim.util.exception.IsNaNException(this, "calcBeta", "beta");
-    }
-    return this.beta[0];
-  }
-
   /** {@inheritDoc} */
   @Override
   public double calcHenrysConstant(String component) {
diff --git a/src/main/java/neqsim/thermodynamicOperations/flashOps/Flash.java b/src/main/java/neqsim/thermodynamicOperations/flashOps/Flash.java
index 5a78671c3..d62664e42 100644
--- a/src/main/java/neqsim/thermodynamicOperations/flashOps/Flash.java
+++ b/src/main/java/neqsim/thermodynamicOperations/flashOps/Flash.java
@@ -353,7 +353,7 @@ public boolean stabilityCheck() {
       }
     } else {
       try {
-        system.calcBeta();
+        RachfordRice.calcBeta(system);
       } catch (Exception ex) {
         logger.error(ex.getMessage(), ex);
       }
diff --git a/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java b/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java
index 72838ee25..86e92be30 100644
--- a/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java
+++ b/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java
@@ -6,8 +6,10 @@
 
 package neqsim.thermodynamicOperations.flashOps;
 
+import neqsim.thermo.ThermodynamicModelSettings;
 import neqsim.thermo.component.ComponentInterface;
 import neqsim.thermo.system.SystemInterface;
+import neqsim.util.exception.IsNaNException;
 
 /**
  * RachfordRice classes.
@@ -16,36 +18,34 @@
  */
 public class RachfordRice {
   private static final long serialVersionUID = 1000;
-  double[] beta = new double[2];
-  SystemInterface fluid = null;
 
-  public RachfordRice(SystemInterface fluid) {
-    this.fluid = fluid;
-  }
-
-  public final double calcBeta() throws neqsim.util.exception.IsNaNException,
+  /**
+   * <p>
+   * calcBeta. For gas liquid systems.
+   * </p>
+   *
+   * @return Beta Mole fraction of gas phase
+   * @throws neqsim.util.exception.IsNaNException if any.
+   * @throws neqsim.util.exception.TooManyIterationsException if any.
+   */
+  public static double calcBeta(SystemInterface fluidinp)
+      throws neqsim.util.exception.IsNaNException,
       neqsim.util.exception.TooManyIterationsException {
+    SystemInterface fluid = fluidinp;
+
     ComponentInterface[] compArray = fluid.getPhase(0).getComponents();
 
     int i;
     double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
-    double deriv = 0.0;
-    double gbeta = 0.0;
-    double betal = 0;
-    double nybeta = 0;
 
+
+    double nybeta = fluid.getBeta(0);
     double midler = 0;
     double minBeta = tolerance;
     double maxBeta = 1.0 - tolerance;
     double g0 = -1.0;
     double g1 = 1.0;
 
-
-    beta[0] = 0.5;
-    beta[1] = 0.5;
-    nybeta = fluid.getBeta(0);
-    betal = 1.0 - nybeta;
-
     for (i = 0; i < fluid.getNumberOfComponents(); i++) {
       midler = (compArray[i].getK() * compArray[i].getz() - 1.0) / (compArray[i].getK() - 1.0);
       if ((midler > minBeta) && (compArray[i].getK() > 1.0)) {
@@ -60,28 +60,19 @@ public final double calcBeta() throws neqsim.util.exception.IsNaNException,
     }
 
     if (g0 < 0) {
-      this.beta[1] = 1.0 - tolerance;
-      this.beta[0] = tolerance;
-      return this.beta[0];
+      fluid.setBeta(tolerance);
+      return tolerance;
     }
     if (g1 > 0) {
-      this.beta[1] = tolerance;
-      this.beta[0] = 1.0 - tolerance;
-      return this.beta[0];
+      fluid.setBeta(1.0 - tolerance);
+      return 1.0 - tolerance;
     }
 
     nybeta = (minBeta + maxBeta) / 2.0;
-    // System.out.println("guessed beta: " + nybeta + " maxbeta: " +maxBeta + "
-    // minbeta: " +minBeta );
-    betal = 1.0 - nybeta;
-
-    // ' *l = 1.0-nybeta;
     double gtest = 0.0;
     for (i = 0; i < fluid.getNumberOfComponents(); i++) {
       gtest += compArray[i].getz() * (compArray[i].getK() - 1.0)
-          / (1.0 - nybeta + nybeta * compArray[i].getK()); // beta
-                                                           // =
-                                                           // nybeta
+          / (1.0 - nybeta + nybeta * compArray[i].getK());
     }
 
     if (gtest >= 0) {
@@ -98,12 +89,13 @@ public final double calcBeta() throws neqsim.util.exception.IsNaNException,
 
     int iterations = 0;
     int maxIterations = 300;
-    // System.out.println("gtest: " + gtest);
     double step = 1.0;
+    double gbeta = 0.0;
+    double deriv = 0.0;
+    double betal = 1.0 - nybeta;
     do {
       iterations++;
       if (gtest >= 0) {
-        // oldbeta = nybeta;
         deriv = 0.0;
         gbeta = 0.0;
 
@@ -122,20 +114,13 @@ public final double calcBeta() throws neqsim.util.exception.IsNaNException,
         }
         nybeta -= (gbeta / deriv);
 
-        // System.out.println("beta: " + maxBeta);
         if (nybeta > maxBeta) {
           nybeta = maxBeta;
         }
         if (nybeta < minBeta) {
           nybeta = minBeta;
         }
-
-        /*
-         * if ((nybeta > maxBeta) || (nybeta < minBeta)) { // nybeta = 0.5 * (maxBeta + minBeta);
-         * gbeta = 1.0; }
-         */
       } else {
-        // oldbeta = betal;
         deriv = 0.0;
         gbeta = 0.0;
 
@@ -160,46 +145,25 @@ public final double calcBeta() throws neqsim.util.exception.IsNaNException,
         if (betal < minBeta) {
           betal = minBeta;
         }
-
-        /*
-         * if ((betal > maxBeta) || (betal < minBeta)) { gbeta = 1.0; { betal = 0.5 * (maxBeta +
-         * minBeta); } }
-         */
         nybeta = 1.0 - betal;
       }
       step = gbeta / deriv;
-      // System.out.println("step : " + step);
-    } while (Math.abs(step) >= 1.0e-10 && iterations < maxIterations); // &&
-    // (Math.abs(nybeta)-Math.abs(maxBeta))>0.1);
-
-    // System.out.println("beta: " + nybeta + " iterations: " + iterations);
+    } while (Math.abs(step) >= 1.0e-10 && iterations < maxIterations);
     if (nybeta <= tolerance) {
-      // this.phase = 1;
       nybeta = tolerance;
     } else if (nybeta >= 1.0 - tolerance) {
-      // this.phase = 0;
       nybeta = 1.0 - tolerance;
-      // superheated vapour
-    } else {
-      // this.phase = 2;
-    } // two-phase liquid-gas
+    }
 
-    this.beta[0] = nybeta;
-    this.beta[1] = 1.0 - nybeta;
-    fluid.setBeta(fluid.getPhaseIndex(0), nybeta);
-    fluid.setBeta(fluid.getPhaseIndex(1), 1.0 - nybeta);
+    fluid.setBeta(nybeta);
 
     if (iterations >= maxIterations) {
-      throw new neqsim.util.exception.TooManyIterationsException(this, "calcBeta", maxIterations);
+      throw new neqsim.util.exception.TooManyIterationsException(fluid, "calcBeta", maxIterations);
     }
-    if (Double.isNaN(beta[1])) {
-      /*
-       * for (i = 0; i < numberOfComponents; i++) { System.out.println("K " + compArray[i].getK());
-       * System.out.println("z " + compArray[i].getz()); }
-       */
-      throw new neqsim.util.exception.IsNaNException(this, "calcBeta", "beta");
+    if (Double.isNaN(nybeta)) {
+      throw new neqsim.util.exception.IsNaNException(fluid, "calcBeta", "beta");
     }
-    return this.beta[0];
+    return nybeta;
   }
 
 }
diff --git a/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java b/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java
index 3871ebfe8..b7b648180 100644
--- a/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java
+++ b/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java
@@ -91,9 +91,7 @@ public void sucsSubs() {
 
     double oldBeta = system.getBeta();
     try {
-      // system.calcBeta();
-      RachfordRice racrice = new RachfordRice(system);
-      racrice.calcBeta();
+      RachfordRice.calcBeta(system);
     } catch (IsNaNException ex) {
       logger.warn("Not able to calculate beta. Value is NaN");
       system.setBeta(oldBeta);
@@ -135,7 +133,7 @@ public void accselerateSucsSubs() {
     }
     double oldBeta = system.getBeta();
     try {
-      system.calcBeta();
+      RachfordRice.calcBeta(system);
     } catch (Exception ex) {
       if (system.getBeta() > 1.0 - betaTolerance || system.getBeta() < betaTolerance) {
         system.setBeta(oldBeta);
@@ -180,7 +178,7 @@ public void resetK() {
       system.getPhase(1).getComponents()[i].setK(Math.exp(lnK[i]));
     }
     try {
-      system.calcBeta();
+      RachfordRice.calcBeta(system);
       system.calc_x_y();
       system.init(1);
     } catch (Exception ex) {
@@ -251,7 +249,7 @@ public void run() {
 
     // Calculates phase fractions and initial composition based on Wilson K-factors
     try {
-      system.calcBeta();
+      RachfordRice.calcBeta(system);
     } catch (Exception ex) {
       logger.error(ex.getMessage(), ex);
     }

From b81e091ab089f06b59c062cc979538fad07987d9 Mon Sep 17 00:00:00 2001
From: Even Solbraa <41290109+EvenSol@users.noreply.github.com>
Date: Wed, 1 May 2024 18:33:00 +0000
Subject: [PATCH 3/4] modified rachford rice implementation

---
 .../neqsim/thermo/system/SystemInterface.java | 15 +++++
 .../neqsim/thermo/system/SystemThermo.java    | 20 +++++++
 .../flashOps/Flash.java                       |  2 +-
 .../flashOps/RachfordRice.java                | 56 +++++++------------
 .../flashOps/TPflash.java                     |  8 +--
 5 files changed, 61 insertions(+), 40 deletions(-)

diff --git a/src/main/java/neqsim/thermo/system/SystemInterface.java b/src/main/java/neqsim/thermo/system/SystemInterface.java
index d786f68d1..3002a2c1e 100644
--- a/src/main/java/neqsim/thermo/system/SystemInterface.java
+++ b/src/main/java/neqsim/thermo/system/SystemInterface.java
@@ -2571,4 +2571,19 @@ public void setImplementedTemperatureDeriativesofFugacity(
    * @param newfile a boolean
    */
   public void write(String name, String filename, boolean newfile);
+
+  /**
+   * <p>
+   * getKvector - return vector of equilibrium constants
+   * </p>
+   */
+  public double[] getKvector();
+
+
+  /**
+   * <p>
+   * getzvector - return vector of total molar composition
+   * </p>
+   */
+  public double[] getzvector();
 }
diff --git a/src/main/java/neqsim/thermo/system/SystemThermo.java b/src/main/java/neqsim/thermo/system/SystemThermo.java
index 7cfedd94b..59f412e4b 100644
--- a/src/main/java/neqsim/thermo/system/SystemThermo.java
+++ b/src/main/java/neqsim/thermo/system/SystemThermo.java
@@ -5021,4 +5021,24 @@ public void write(String name, String filename, boolean newfile) {
     file.setValues(table);
     file.createFile();
   }
+
+  /** {@inheritDoc} */
+  @Override
+  public double[] getKvector() {
+    double[] K = new double[this.getNumberOfComponents()];
+    for (int i = 0; i < this.getNumberOfComponents(); i++) {
+      K[i] = this.getPhase(0).getComponent(i).getK();
+    }
+    return K;
+  }
+
+  /** {@inheritDoc} */
+  @Override
+  public double[] getzvector() {
+    double[] z = new double[this.getNumberOfComponents()];
+    for (int i = 0; i < this.getNumberOfComponents(); i++) {
+      z[i] = this.getPhase(0).getComponent(i).getz();
+    }
+    return z;
+  }
 }
diff --git a/src/main/java/neqsim/thermodynamicOperations/flashOps/Flash.java b/src/main/java/neqsim/thermodynamicOperations/flashOps/Flash.java
index d62664e42..c8f47b3f4 100644
--- a/src/main/java/neqsim/thermodynamicOperations/flashOps/Flash.java
+++ b/src/main/java/neqsim/thermodynamicOperations/flashOps/Flash.java
@@ -353,7 +353,7 @@ public boolean stabilityCheck() {
       }
     } else {
       try {
-        RachfordRice.calcBeta(system);
+        system.setBeta(RachfordRice.calcBeta(system.getKvector(), system.getzvector()));
       } catch (Exception ex) {
         logger.error(ex.getMessage(), ex);
       }
diff --git a/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java b/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java
index 86e92be30..89728f8f4 100644
--- a/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java
+++ b/src/main/java/neqsim/thermodynamicOperations/flashOps/RachfordRice.java
@@ -28,51 +28,41 @@ public class RachfordRice {
    * @throws neqsim.util.exception.IsNaNException if any.
    * @throws neqsim.util.exception.TooManyIterationsException if any.
    */
-  public static double calcBeta(SystemInterface fluidinp)
-      throws neqsim.util.exception.IsNaNException,
+  public static double calcBeta(double[] K, double[] z) throws neqsim.util.exception.IsNaNException,
       neqsim.util.exception.TooManyIterationsException {
-    SystemInterface fluid = fluidinp;
-
-    ComponentInterface[] compArray = fluid.getPhase(0).getComponents();
 
     int i;
     double tolerance = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
-
-
-    double nybeta = fluid.getBeta(0);
     double midler = 0;
     double minBeta = tolerance;
     double maxBeta = 1.0 - tolerance;
     double g0 = -1.0;
     double g1 = 1.0;
 
-    for (i = 0; i < fluid.getNumberOfComponents(); i++) {
-      midler = (compArray[i].getK() * compArray[i].getz() - 1.0) / (compArray[i].getK() - 1.0);
-      if ((midler > minBeta) && (compArray[i].getK() > 1.0)) {
+    for (i = 0; i < K.length; i++) {
+      midler = (K[i] * z[i] - 1.0) / (K[i] - 1.0);
+      if ((midler > minBeta) && (K[i] > 1.0)) {
         minBeta = midler;
       }
-      midler = (1.0 - compArray[i].getz()) / (1.0 - compArray[i].getK());
-      if ((midler < maxBeta) && (compArray[i].getK() < 1.0)) {
+      midler = (1.0 - z[i]) / (1.0 - K[i]);
+      if ((midler < maxBeta) && (K[i] < 1.0)) {
         maxBeta = midler;
       }
-      g0 += compArray[i].getz() * compArray[i].getK();
-      g1 += -compArray[i].getz() / compArray[i].getK();
+      g0 += z[i] * K[i];
+      g1 += -z[i] / K[i];
     }
 
     if (g0 < 0) {
-      fluid.setBeta(tolerance);
       return tolerance;
     }
     if (g1 > 0) {
-      fluid.setBeta(1.0 - tolerance);
       return 1.0 - tolerance;
     }
 
-    nybeta = (minBeta + maxBeta) / 2.0;
+    double nybeta = (minBeta + maxBeta) / 2.0;
     double gtest = 0.0;
-    for (i = 0; i < fluid.getNumberOfComponents(); i++) {
-      gtest += compArray[i].getz() * (compArray[i].getK() - 1.0)
-          / (1.0 - nybeta + nybeta * compArray[i].getK());
+    for (i = 0; i < K.length; i++) {
+      gtest += z[i] * (K[i] - 1.0) / (1.0 - nybeta + nybeta * K[i]);
     }
 
     if (gtest >= 0) {
@@ -99,12 +89,11 @@ public static double calcBeta(SystemInterface fluidinp)
         deriv = 0.0;
         gbeta = 0.0;
 
-        for (i = 0; i < fluid.getNumberOfComponents(); i++) {
-          double temp1 = (compArray[i].getK() - 1.0);
+        for (i = 0; i < K.length; i++) {
+          double temp1 = (K[i] - 1.0);
           double temp2 = 1.0 + temp1 * nybeta;
-          deriv += -(compArray[i].getz() * temp1 * temp1) / (temp2 * temp2);
-          gbeta += compArray[i].getz() * (compArray[i].getK() - 1.0)
-              / (1.0 + (compArray[i].getK() - 1.0) * nybeta);
+          deriv += -(z[i] * temp1 * temp1) / (temp2 * temp2);
+          gbeta += z[i] * (K[i] - 1.0) / (1.0 + (K[i] - 1.0) * nybeta);
         }
 
         if (gbeta >= 0) {
@@ -124,11 +113,9 @@ public static double calcBeta(SystemInterface fluidinp)
         deriv = 0.0;
         gbeta = 0.0;
 
-        for (i = 0; i < fluid.getNumberOfComponents(); i++) {
-          deriv -= (compArray[i].getz() * (compArray[i].getK() - 1.0) * (1.0 - compArray[i].getK()))
-              / Math.pow((betal + (1 - betal) * compArray[i].getK()), 2);
-          gbeta += compArray[i].getz() * (compArray[i].getK() - 1.0)
-              / (betal + (-betal + 1.0) * compArray[i].getK());
+        for (i = 0; i < K.length; i++) {
+          deriv -= (z[i] * (K[i] - 1.0) * (1.0 - K[i])) / Math.pow((betal + (1 - betal) * K[i]), 2);
+          gbeta += z[i] * (K[i] - 1.0) / (betal + (-betal + 1.0) * K[i]);
         }
 
         if (gbeta < 0) {
@@ -155,13 +142,12 @@ public static double calcBeta(SystemInterface fluidinp)
       nybeta = 1.0 - tolerance;
     }
 
-    fluid.setBeta(nybeta);
-
     if (iterations >= maxIterations) {
-      throw new neqsim.util.exception.TooManyIterationsException(fluid, "calcBeta", maxIterations);
+      throw new neqsim.util.exception.TooManyIterationsException(new RachfordRice(), "calcBeta",
+          maxIterations);
     }
     if (Double.isNaN(nybeta)) {
-      throw new neqsim.util.exception.IsNaNException(fluid, "calcBeta", "beta");
+      throw new neqsim.util.exception.IsNaNException(new RachfordRice(), "calcBeta", "beta");
     }
     return nybeta;
   }
diff --git a/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java b/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java
index b7b648180..8698b6923 100644
--- a/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java
+++ b/src/main/java/neqsim/thermodynamicOperations/flashOps/TPflash.java
@@ -91,7 +91,7 @@ public void sucsSubs() {
 
     double oldBeta = system.getBeta();
     try {
-      RachfordRice.calcBeta(system);
+      system.setBeta(RachfordRice.calcBeta(system.getKvector(), system.getzvector()));
     } catch (IsNaNException ex) {
       logger.warn("Not able to calculate beta. Value is NaN");
       system.setBeta(oldBeta);
@@ -133,7 +133,7 @@ public void accselerateSucsSubs() {
     }
     double oldBeta = system.getBeta();
     try {
-      RachfordRice.calcBeta(system);
+      system.setBeta(RachfordRice.calcBeta(system.getKvector(), system.getzvector()));
     } catch (Exception ex) {
       if (system.getBeta() > 1.0 - betaTolerance || system.getBeta() < betaTolerance) {
         system.setBeta(oldBeta);
@@ -178,7 +178,7 @@ public void resetK() {
       system.getPhase(1).getComponents()[i].setK(Math.exp(lnK[i]));
     }
     try {
-      RachfordRice.calcBeta(system);
+      system.setBeta(RachfordRice.calcBeta(system.getKvector(), system.getzvector()));
       system.calc_x_y();
       system.init(1);
     } catch (Exception ex) {
@@ -249,7 +249,7 @@ public void run() {
 
     // Calculates phase fractions and initial composition based on Wilson K-factors
     try {
-      RachfordRice.calcBeta(system);
+      system.setBeta(RachfordRice.calcBeta(system.getKvector(), system.getzvector()));
     } catch (Exception ex) {
       logger.error(ex.getMessage(), ex);
     }

From 4577e979d2a2ffc74d897225f0eb7cf710239985 Mon Sep 17 00:00:00 2001
From: Even Solbraa <41290109+EvenSol@users.noreply.github.com>
Date: Wed, 1 May 2024 18:36:38 +0000
Subject: [PATCH 4/4] added Rachford rice tests

---
 .../flashOps/RachfordRiceTest.java            | 20 +++++++++++++++++++
 1 file changed, 20 insertions(+)
 create mode 100644 src/test/java/neqsim/thermodynamicOperations/flashOps/RachfordRiceTest.java

diff --git a/src/test/java/neqsim/thermodynamicOperations/flashOps/RachfordRiceTest.java b/src/test/java/neqsim/thermodynamicOperations/flashOps/RachfordRiceTest.java
new file mode 100644
index 000000000..81149e8f9
--- /dev/null
+++ b/src/test/java/neqsim/thermodynamicOperations/flashOps/RachfordRiceTest.java
@@ -0,0 +1,20 @@
+package neqsim.thermodynamicOperations.flashOps;
+
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+
+public class RachfordRiceTest {
+  @Test
+  void testCalcBeta() {
+
+    double[] z = new double[] {0.1, 0.9};
+    double[] K = new double[] {20.0, 0.1};
+
+    try {
+      Assertions.assertEquals(0.06374269005, RachfordRice.calcBeta(K, z), 1e-6);
+    } catch (Exception e) {
+      e.printStackTrace();
+    }
+
+  }
+}