Burgoyne-Nielsen-Stanko Viscosity

/**
 * Gas viscosity using tuned Lohrenz–Bray–Clark (BNS implementation). Adjusted critical volumes and coefficients are
 * taken from table 11, equation 10 and table 15 in SPE-229932.
 *
 * @param pressure the gas pressure in psia
 * @param temperature the gas temperature in degrees Fahrenheit
 * @param sg gas specific gravity (air = 1)
 * @param co2 mole fraction of carbon dioxide
 * @param h2s mole fraction of hydrogen sulphide
 * @param n2 mole fraction of nitrogen
 * @param h2 mole fraction of hydrogen
 * @return the viscosity of the fluid in centiPoise
 */
public static double gasViscosityBNS(double pressure, double temperature, double sg, double co2, double h2s,
        double n2, double h2) {
    pressure = min(pressure, 15000.0); // limit of validity for this correlation is 15,000 psia
    final int N = 5; // always 5 components in BNS composition
    double t_degr = Amount.valueOf(temperature, FAHRENHEIT).doubleValue(RANKINE);
    double p_psia = pressure;
    double[] mws = BNS_MWS.clone();
    double[] tcs = BNS_TCS.clone();
    double[] pcs = BNS_PCS.clone();
    double[] vcvis = BNS_VCVIS.clone();
    FluidComponent[] fcs = {PureComponent.CO2, PureComponent.H2S, PureComponent.N2, PureComponent.H2,
        PureComponent.CH4}; // technically last component is a pseudo-hydrocarbon component
    double[] z = {co2, h2s, n2, h2, 1.0 - (co2 + h2s + n2 + h2)}; // build composition vector
    double sg_hc = hydrocarbonSpecificGravity(sg, co2, h2s, n2, h2);
    double hc_gas_mw = sg_hc * PureComponent.AIR.mass();

    // Override pseudo-hydrocarbon MW
    int gas_idx = N - 1; // index 4
    mws[gas_idx] = hc_gas_mw;

    // BNS Tc,Pc for a pure pseudo-hydrocarbon gas (no non-hc components)
    double tpc_hc = criticalTemperatureBNS(sg_hc, 0.0, 0.0, 0.0, 0.0, false); // degR
    double ppc_hc = criticalPressureBNS(sg_hc, 0.0, 0.0, 0.0, 0.0, false); // psia
    tcs[gas_idx] = tpc_hc;
    pcs[gas_idx] = ppc_hc;
    vcvis[gas_idx] = 0.0576710 * (hc_gas_mw - PureComponent.CH4.mass()) + 1.44383; // eq 10 in SPE-229932

    // Dilute gas viscosity from Stiel-Thodos
    FastMap<FluidComponent, Amount<DynamicViscosity>> xi_mui = new FastMap<>(N);
    FastMap<FluidComponent, Double> comp_bns = new FastMap<>(N);
    for (int i = 0; i < N; i++) {
        double tri = t_degr / tcs[i];
        double tc_k = Amount.valueOf(tcs[i], RANKINE).doubleValue(KELVIN);
        double pc_atm = Amount.valueOf(pcs[i], POUND_PER_SQUARE_INCH).doubleValue(ATMOSPHERE);
        double xi = PVT.viscosityReducingParameter(tc_k, pc_atm, mws[i]);
        double mui_cp = PVT.gasViscosityLowPressureStielThodos(tri, xi);
        xi_mui.put(fcs[i], Amount.valueOf(mui_cp, CENTIPOISE));
        comp_bns.put(fcs[i], z[i]);
    }

    // Apply Herning and Zipperer (1936) mixing rule to get mixture dilute gas viscosity.
    Amount<DynamicViscosity> mu_star = mixingRuleHerningZipperer(xi_mui, comp_bns);

    // Determine mixture values that will be used for LBC
    double vcvis_mix = 0.0;
    double tc_k_mix = 0.0;
    double mw_mix = 0.0;
    double pc_atm_mix = 0.0;
    for (int i = 0; i < N; i++) {
        double tc_k_i = Amount.valueOf(tcs[i], RANKINE).doubleValue(KELVIN);
        double pc_atm_i = Amount.valueOf(pcs[i], POUND_PER_SQUARE_INCH).doubleValue(ATMOSPHERE);
        vcvis_mix += vcvis[i] * z[i];
        tc_k_mix += z[i] * tc_k_i;
        mw_mix += z[i] * mws[i];
        pc_atm_mix += z[i] * pc_atm_i;
    }

    // Apply modified LBC form to determine viscosity.
    double rhoc = 1.0 / vcvis_mix;
    double zfac_bns = zFactorBNS(pressure, temperature, sg, co2, h2s, n2, h2);
    double rho_r = p_psia / (zfac_bns * Mixture.GAS_CONSTANT * t_degr * rhoc);
    double xi = PVT.viscosityReducingParameter(tc_k_mix, pc_atm_mix, mw_mix);
    double mu_bns = PVT.viscosityBNSModifiedLBC(rho_r, xi, mu_star.doubleValue(CENTIPOISE));
    return mu_bns;
}

/**
 * Based on method of Jossi, this correlates the residual viscosity to the reduced density but with the last two
 * terms modified by Burgoyne, Nielsen and Stanko (2025). This is specifically used with the BNS viscosity method.
 *
 * @param rho_r reduced density
 * @param xi viscosity reducing parameter (1/cP)
 * @param mu_star dilute gas viscosity at target temperature (cP)
 * @return viscosity at target pressure and temperature in centiPoise
 */
protected static final double viscosityBNSModifiedLBC(double rho_r, double xi, double mu_star) {
    final double[] A_LBC = {0.1023, 0.023364, 0.058533, -0.0392852, 0.00926279};
    double mu = A_LBC[0] + rho_r * (A_LBC[1] + rho_r * (A_LBC[2] + rho_r * (A_LBC[3] + rho_r * A_LBC[4])));
    mu = pow(mu, 4.0);
    mu -= 1.0e-04;
    mu /= xi;
    mu += mu_star;
    return mu;
}

/**
 * Uses BNS tuned Peng Robinson EOS to solve for Z factor. Ported from z_bur at
 * https://github.com/mwburgoyne/pyResToolbox/blob/main/pyrestoolbox/gas/gas.py.
 *
 * @param pressure the gas pressure in psia
 * @param temperature the gas temperature in degrees Fahrenheit
 * @param sg gas specific gravity (air = 1)
 * @param co2 mole fraction of carbon dioxide
 * @param h2s mole fraction of hydrogen sulphide
 * @param n2 mole fraction of nitrogen
 * @param h2 mole fraction of hydrogen
 * @return the Z factor
 */
public static double zFactorBNS(double pressure, double temperature, double sg, double co2, double h2s, double n2,
        double h2) {
    final int N = 5; // always 5 components in BNS composition
    double t_degr = Amount.valueOf(temperature, FAHRENHEIT).doubleValue(RANKINE);
    double p_psia = pressure;
    double[] tcs = BNS_TCS.clone();
    double[] pcs = BNS_PCS.clone();
    double[] acf = BNS_ACF.clone();

    /*
     * Workflow from SPE-229932 for z-factor property calculation is:
     * 
     *  1. Define mole fraction composition.
     *  2. Override pseudo-hydrocarbon Tc, Pc using BNS correlations into arrays of tuned critical temperatures and
     *     pressures.
     *  3. Calculate temperature-dependent pairs of BIPs (kij matrix) at temperature of interest using Tc of 
     *     pseudo-hydrocarbon.
     *  4. Use the PR EOS to solve for Z-factor analytically at the given P, T. For most spreadsheet style
     *     applications, the largest root can be taken as the gas phase Z-factor. In this implementation we apply
     *     the same trick as employed by Burgoyne to determine if there are actually three roots. If this is the
     *     case then we determine the lowest fugacity coefficient to select the correct root.
     */
    // 1. Define mole fraction composition.
    double[] z = {co2, h2s, n2, h2, 1.0 - (co2 + h2s + n2 + h2)}; // build composition vector

    // 2. Override pseudo-hydrocarbon Tc, Pc using BNS correlations into arrays of tuned critical temperatures and
    //    pressures
    double tpc_hc = criticalTemperatureBNS(sg, co2, h2s, n2, h2, false); // degR
    double ppc_hc = criticalPressureBNS(sg, co2, h2s, n2, h2, false); // psia
    int gas_idx = N - 1; // index 4
    tcs[gas_idx] = tpc_hc;
    pcs[gas_idx] = ppc_hc;

    // 3. Calculate temperature-dependent pairs of BIPs (kij matrix) at temperature of interest using Tc of 
    //    pseudo-hydrocarbon.
    double[][] kij = calcBipsBNS(t_degr, tpc_hc);

    // 4. Use the PR EOS to solve for Z-factor analytically at the given P, T. For most spreadsheet style
    //    applications, the largest root can be taken as the gas phase Z-factor.
    // Compute component PR parameters ai, bi.
    double[] ai = new double[N];
    double[] bi = new double[N];
    for (int i = 0; i < N; i++) {
        double tci_degr = tcs[i];
        double tri = t_degr / tci_degr; // component reduced temperature
        double sqrt_tr = Math.sqrt(tri);
        double pci = pcs[i];
        double pri = p_psia / pci; // component reduced pressure
        double omega = acf[i];

        // Standard PR alpha and m(omega).
        double m = 0.37464 + 1.54226 * omega - 0.26992 * omega * omega;
        double alpha = Math.pow(1.0 + m * (1.0 - sqrt_tr), 2.0);
        ai[i] = BNS_OMEGAA[i] * alpha * pri / (tri * tri);
        bi[i] = BNS_OMEGAB[i] * pri / tri;
    }

    // Mixture a_mix (cubic A) and b_mix (cubic B) with BIPs.
    // Mixing rule A: A = dblsum_ij z_i * z_j * sqrt(Ai_i * Ai_j) * (1 - kij)
    // Mixing rule B: B = sum_i z_i * Bi_i
    double a_mix = 0.0;
    double b_mix = 0.0;
    for (int i = 0; i < N; i++) {
        b_mix += z[i] * bi[i];
        for (int j = 0; j < N; j++) {
            double aij = Math.sqrt(ai[i] * ai[j]) * (1.0 - kij[i][j]);
            a_mix += z[i] * z[j] * aij;
        }
    }
    final double A = a_mix;
    final double B = b_mix;

    // Cubic EOS coefficients for PR: Z^3 + c2 Z^2 + c1 Z + c0 = 0
    double c2 = -(1.0 - B);
    double c1 = A - 3.0 * B * B - 2.0 * B;
    double c0 = B * (B * B + B - A);

    // Solve cubic with Halley’s method (single vapour root).
    double z_root = AdvMath.halleyCubic(c2, c1, c0);

    // This is a fugacity-based root selection for sub-critical conditions. When 3 real roots exist, the 
    // thermodynamically stable phase is the one with the lowest fugacity coefficient (Gibbs criterion). First we 
    // find the discriminant of depressed cubic in order to check for 3 real roots.
    double p_d = (3.0 * c1 - c2 * c2) / 3.0;
    double q_d = (2.0 * c2 * c2 * c2 - 9.0 * c2 * c1 + 27.0 * c0) / 27.0;
    double disc = q_d * q_d / 4.0 + p_d * p_d * p_d / 27.0;

    // If there are three real roots then we are subcritical and must do fugacity-based root selection. Discriminant
    // less than zero implies three real roots.
    if (disc < -1.0e-15) {
        double z_max = z_root;

        // Deflate cubic by z_max to get quadratic for remaining roots.
        double b_q = c2 + z_max;
        double c_q = c1 + z_max * b_q;
        double det = max(b_q * b_q - 4.0 * c_q, 0.0);
        double sqrt_det = Math.sqrt(det);
        double z_min = (-b_q - sqrt_det) / 2.0; // smallest root

        // PR ln(phi) for root selection.
        double sqrt2 = Math.sqrt(2.0);
        double s2p1 = 1.0 + sqrt2;
        double s2m1 = sqrt2 - 1.0;

        // Only compare where z_min is physically valid (Z > B)
        if (z_min > B) {
            double lnphi_max = ((z_max - 1.0)
                    - Math.log(z_max - B)
                    - A / (2.0 * sqrt2 * B)
                    * Math.log((z_max + s2p1 * B) / (z_max - s2m1 * B)));

            double lnphi_min = ((z_min - 1.0)
                    - Math.log(z_min - B)
                    - A / (2.0 * sqrt2 * B)
                    * Math.log((z_min + s2p1 * B) / (z_min - s2m1 * B)));
            if (lnphi_min < lnphi_max) {
                z_root = z_min;
            }
        }
    }

    // Apply volume shift for BNS-style shifted Z
    double vshift_mix = 0.0;
    for (int i = 0; i < N; i++) {
        vshift_mix += z[i] * BNS_VSHIFT[i] * bi[i];
    }
    return z_root - vshift_mix;
}

/**
 * Component viscosity-reducing parameter ≈ 1 / μci (cP-1).
 *
 * @param tc critical temperature, K
 * @param pc critical pressure, atm
 * @param mw molecular weight, g/mol
 * @return viscosity-reducing parameter ≈ 1 / μci (1/cP)
 */
public static final double viscosityReducingParameter(double tc, double pc, double mw) {
    double xi = pow(tc, 1.0 / 6.0) / (sqrt(mw) * pow(pc, 2.0 / 3.0));
    return xi;
}

/**
 * Calculates the dilute gas viscosity (typically at low pressure 1 to 5 atm) for a component.
 *
 * @param tr reduced temperature
 * @param xi viscosity reducing parameter (1/cP)
 * @return dilute gas viscosity (cP)
 */
public static final double gasViscosityLowPressureStielThodos(double tr, double xi) {
    double mu_star;
    if (tr < 1.5) {
        mu_star = (34.0e-05 * pow(tr, 0.94)) / xi;
    } else {
        mu_star = (17.78e-05 * pow(max(4.58 * tr - 1.67, 1.0e-30), 5.0 / 8.0)) / xi;
    }
    return mu_star;
}

/**
 * In a simplification of the kinetic theory approach, ignoring second order effects, Herning and Zipperer (1936)
 * proposed a viscosity mixing rule based on molecular weights.
 *
 * @param xi_mui map of fluid components to pure component viscosity (allows different pure component viscosities)
 * @param comp map of fluid components to pure component mixture fraction
 * @return viscosity of the mixture
 */
public static Amount<DynamicViscosity> mixingRuleHerningZipperer(
        FastMap<FluidComponent, Amount<DynamicViscosity>> xi_mui, FastMap<FluidComponent, Double> comp) {
    double mui_zi_sqrtmi = 0.0;
    for (FluidComponent fc_i : comp.keySet()) {
        final Double zi = comp.get(fc_i);
        final double m_i = fc_i.mass();
        double zj_phi_ij = 0.0;
        for (FluidComponent fc_j : comp.keySet()) {
            final Double zj = comp.get(fc_j);
            final double m_j = fc_j.mass();
            zj_phi_ij += zj * sqrt(m_j / m_i);
        }
        mui_zi_sqrtmi += (xi_mui.get(fc_i).doubleValue(CENTIPOISE) * zi) / zj_phi_ij;
    }
    return Amount.valueOf(mui_zi_sqrtmi, CENTIPOISE);
}

protected static double hydrocarbonSpecificGravity(double sg, double co2, double h2s, double n2, double h2) {
    double sg_hc = 0.75; // default value if 100% inerts to avoid numerical problems
    if (co2 + h2s + n2 + h2 < 1.0) {
        final double MW_CO2 = PureComponent.CO2.mass();
        final double MW_H2S = PureComponent.H2S.mass();
        final double MW_N2 = PureComponent.N2.mass();
        final double MW_H2 = PureComponent.H2.mass();
        sg_hc = (sg - (co2 * MW_CO2 + h2s * MW_H2S + n2 * MW_N2 + h2 * MW_H2) / PureComponent.AIR.mass())
                / (1.0 - co2 - h2s - n2 - h2);
        sg_hc = max(PureComponent.CH4.specificGravity(), sg_hc);
    }
    return sg_hc;
}

/**
 * Calculates pseudo-critical temperature for the pseudo-hydrocarbon component. Uses custom linear fits for gas
 * condensates and associated gas (Burgoyne, 2025). These relations are derived from property fitting and may differ
 * from literature.
 *
 * @param sg gas specific gravity (air = 1.0)
 * @param co2 mole fraction of carbon dioxide
 * @param h2s mole fraction of hydrogen sulphide
 * @param n2 mole fraction of nitrogen
 * @param h2 mole fraction of hydrogen
 * @param is_ag flag for if this associated gas (true) or gas-condensate (false)
 * @return the critical temperature for the pseudo-hydrocarbon component in degrees Rankine (degR)
 */
public static double criticalTemperatureBNS(double sg, double co2, double h2s, double n2, double h2,
        boolean is_ag) {
    int idx = is_ag ? 0 : 1;
    double sg_hc = hydrocarbonSpecificGravity(sg, co2, h2s, n2, h2);
    final double[] A = {2695.14765, 1098.10948};
    final double[] B = {274.341701, 101.529237};
    final double[] C = {PureComponent.CH4.t_crit(), PureComponent.CH4.t_crit()};
    double x = max(0.0, PureComponent.AIR.mass() * sg_hc - PureComponent.CH4.mass());
    double tpc_hc = A[idx] * x / (B[idx] + x) + C[idx];
    return tpc_hc;
}

/**
 * Calculates pseudo-critical pressure for the pseudo-hydrocarbon component. Uses custom linear fits for gas
 * condensates and associated gas (Burgoyne, 2025). These relations are derived from property fitting and may differ
 * from literature.
 *
 * @param sg gas specific gravity (air = 1.0)
 * @param co2 mole fraction of carbon dioxide
 * @param h2s mole fraction of hydrogen sulphide
 * @param n2 mole fraction of nitrogen
 * @param h2 mole fraction of hydrogen
 * @param is_ag flag for if this associated gas (true) or gas-condensate (false)
 * @return the critical pressure for the pseudo-hydrocarbon component in psia
 */
public static double criticalPressureBNS(double sg, double co2, double h2s, double n2, double h2, boolean is_ag) {
    int idx = is_ag ? 0 : 1;
    double sg_hc = hydrocarbonSpecificGravity(sg, co2, h2s, n2, h2);
    final double[] VC_SLOPE = {0.177497835, 0.170931432};
    double x = max(0.0, PureComponent.AIR.mass() * sg_hc - PureComponent.CH4.mass());
    double tpc_hc = criticalTemperatureBNS(sg, co2, h2s, n2, h2, is_ag);
    double vc_on_zc = VC_SLOPE[idx] * x + 5.518525872412144;
    double ppc_hc = Mixture.GAS_CONSTANT * tpc_hc / vc_on_zc;
    return ppc_hc;
}

/**
 * Compute 5x5 BIP matrix using precomputed constants. From equation 4 in paper SPE-229932.
 *
 * @param deg_r temperature in Rankine
 * @param tpc_hc pseudo-critical temperature of the gas (same units as BNS_TCS)
 * @return 5x5 kij matrix
 */
private static double[][] calcBipsBNS(double deg_r, double tpc_hc) {
    final int N = 5;
    double[][] slope_tc = new double[N][N];
    for (int i = 0; i < N; i++) {
        System.arraycopy(BIP_SLOPE_TC[i], 0, slope_tc[i], 0, N);
    }
    for (int j = 0; j < 4; j++) {
        slope_tc[4][j] = BIP_GAS_SLOPES[j] * tpc_hc;
    }
    for (int i = 0; i < 4; i++) {
        slope_tc[i][4] = slope_tc[4][i];
    }
    double[][] kij = new double[N][N];
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            kij[i][j] = BIP_CONST[i][j] + slope_tc[i][j] / deg_r;
        }
    }
    return kij;
}

/**
 * Halley solver for cubic: x^3 + c2 * x^2 + c1 * x + c0 = 0. Returns the largest real root. Falls back to 
 * cubicSolver() if Halley does not converge.
 */
public static double halleyCubic(double c2, double c1, double c0) {
    final int MAXITER = 50;
    final double TOL = 1.0e-12;
    double x = -c2 / 3.0; // Initial guess: Z = -c2/3
    double f0 = x * x * x + c2 * x * x + c1 * x + c0;
    if (f0 < 0.0) {
        x += 1.0;
    }

    // Halley iterations
    for (int iter = 0; iter < MAXITER; iter++) {
        double f = (x + c2) * x * x + c1 * x + c0; // x^3 + c2 * x^2 + c1 * x + c0
        double fp = (3.0 * x + 2.0 * c2) * x + c1; // derivative
        double fpp = 6.0 * x + 2.0 * c2; // second derivative

        // Protect against zero derivative
        double safe_fp = (Math.abs(fp) < 1.0e-30) ? 1.0e-30 : fp;
        double dx = f / safe_fp;
        double denom = safe_fp - 0.5 * dx * fpp;
        dx = f / ((Math.abs(denom) < 1.0e-30) ? 1.0e-30 : denom);
        x -= dx;
        if (Math.abs(dx) < TOL) {
            break;
        }
    }

    // Check residual and fall back to Cardano for any bad element
    double f = x * x * x + c2 * x * x + c1 * x + c0;
    if (Math.abs(f) > 1.0e-6) {
        double[] roots = cubicSolver(1.0, c2, c1, c0);

        // Return the largest real root (vapor root)
        double xmax = roots[0];
        for (double r : roots) {
            if (r > xmax) {
                xmax = r;
            }
        }
        return xmax;
    }
    return x;
}

/**
 * Solves cubic polynomial of form a0 * x^3 + b0 * x^2 + c0 * x + d0 = 0. Uses an analytical scheme to determine the
 * roots of the cubic equation.
 *
 * @param a0 first cubic co-efficient
 * @param b0 second cubic co-efficient
 * @param c0 third cubic co-efficient
 * @param d0 fourth cubic co-efficient
 * @return array containing all the real roots for this cubic which can be 1, 2 or 3 in number
 */
public static double[] cubicSolver(double a0, double b0, double c0, double d0) {
    double a = 0.0, b, c, p, q, q2, p3, r, aa, bb;
    double[] y = new double[3];
    int n;

    // Check if we have invalid cubic equation
    if (a0 == 0.0) {
        n = 0;
    } else {
        a = b0 / a0;
        b = c0 / a0;
        c = d0 / a0;
        p = (3.0 * b - a * a) / 3.0;
        q = (27.0 * c - 9.0 * a * b + 2.0 * a * a * a) / 27.0;
        p3 = p * p * p / 27.0;
        q2 = q * q / 4.0;
        r = p3 + q2;
        if (r < 0.0) {
            n = 3;
            findThreeRoots(p, q, q2, p3, y);
        } else {
            n = 1;
            aa = cubeRoot(Math.sqrt(r) - q / 2.0);
            bb = cubeRoot(-Math.sqrt(r) - q / 2.0);
            y[0] = aa + bb;
            if (r == 0.0) {
                n = 2;
                y[1] = -aa;
            }
        }
    }

    // Create the output array to contain the real roots
    double[] x = new double[n];
    for (int i = 0; i < n; i++) {
        x[i] = y[i] - a / 3;
    }
    return x;
}

private static double cubeRoot(double x) {
    if (x < 0.0) {
        return -Math.exp((Math.log(-x)) / 3.0);
    } else if (x > 0.0) {
        return Math.exp((Math.log(x)) / 3.0);
    }

    // Case x == 0.0
    return 0.0;
}

private static void findThreeRoots(double p, double q, double q2, double p3, double[] y) {
    double phi, cos2phi, tanphi;
    if (q == 0.0) {
        phi = Math.PI / 2.0;
    } else {
        cos2phi = -q2 / p3;
        tanphi = Math.sqrt(1 / cos2phi - 1.0);
        phi = Math.atan(tanphi);
    }
    for (int k = 0; k < 3; k++) {
        y[k] = 2 * Math.sqrt(-p / 3.0) * Math.cos((phi + 2.0 * Math.PI * (k - 1.0)) / 3.0);
        if (q > 0.0) {
            y[k] = -y[k];
        }
    }
}

/*
 * EOS parameters for BNS Peng-Robinson model (shared by Z-factor and viscosity functions)
 */
private static final double[] BNS_MWS = {44.01, 34.082, 28.014, 2.016, 0.0};
private static final double[] BNS_TCS = {547.416, 672.120, 227.160, 47.430, 1.0}; // H2 Tc has been modified
private static final double[] BNS_PCS = {1069.51, 1299.97, 492.84, 187.5300, 1.0};
private static final double[] BNS_ACF = {0.12253, 0.04909, 0.037, -0.21700, -0.03899};
private static final double[] BNS_VSHIFT = {-0.27607, -0.22901, -0.21066, -0.36270, -0.19076};
private static final double[] BNS_OMEGAA = {0.427671, 0.436725, CubicEOS.OMEGA_A_PR, CubicEOS.OMEGA_A_PR,
    CubicEOS.OMEGA_A_PR};
private static final double[] BNS_OMEGAB = {0.0696397, 0.0724345, CubicEOS.OMEGA_B_PR, CubicEOS.OMEGA_B_PR,
    CubicEOS.OMEGA_B_PR};
private static final double[] BNS_VCVIS = {1.46352, 1.46808, 1.35526, 0.68473, 0.0}; // cuft/lbm-mol
/*
 * BIP precomputed matrices: kij[i, j] = BIP_CONST[i, j] + BIP_SLOPE_TC[i, j] / degR
 * Component index order: CO2 = 0, H2S = 1, N2 = 2, H2 = 3, Gas = 4
 * Gas column/row uses tpc_hc at runtime; stored slopes in BIP_GAS_SLOPES
 */
private static final double[][] BIP_CONST = {
    {0.0, 0.248638, -0.25, -0.247153, -0.145561},
    {0.248638, 0.0, -0.204414, 0.0, 0.16852},
    {-0.25, -0.204414, 0.0, -0.166253, -0.108},
    {-0.247153, 0.0, -0.166253, 0.0, -0.0620119},
    {-0.145561, 0.16852, -0.108, -0.0620119, 0.0}
};
private static final double[][] BIP_SLOPE_TC = {
    {0.0, -75.64467996, 63.51120432, 89.65031832, 0.0},
    {-75.64467996, 0.0, 157.55635404, 0.0, 0.0},
    {63.51120432, 157.55635404, 0.0, 17.90313836, 0.0},
    {89.65031832, 0.0, 17.90313836, 0.0, 0.0},
    {0.0, 0.0, 0.0, 0.0, 0.0}
};
private static final double[] BIP_GAS_SLOPES = {0.276572, -0.122378, 0.0605506, 0.0427873};