Z-factor: Dranchuk-Purvis-Robinson correlation C# code

Code to calculate Z factor with Dranchuk-Purvis-Robinson Method in C#.

Dranchuk et. al. (1974) developed a correlation based on the Benedict-Webb-Rubin type of equation of state. Fittig the equation to 1500 data points from the Standing and Katz Z-factor chart optimized the eight coefficient of the proposed equations. The equation has the following form:

The reduced density is defined by the next equation:

And the coefficients have the following values:

The first equation is solved using the Newton-Rhapson iteration technique, getting:

and the derivative is obteined by evaluating the next equation:

The process is repeated n times until the error becomes smaller than a preset tolerance.

       

            static void Main(string[] args)
        {
            //Dranchuk purvis robinson
            //Coefficients
            double A1 = 0.31506237;
            double A2 = -1.0467099;
            double A3 = -0.57832729;
            double A4 = 0.53530771;
            double A5 = -0.61232032;
            double A6 = -0.10488813;
            double A7 = 0.68157001;
            double A8 = 0.68446549;
            //Result,zequation,deritive, reduced density, pseudocritical pressure, pseudocritical temperature, actual z, next x value
            double DPRc, fz, dfz, dr, Ppr, Tpr, zt, ztdt;
            double dfz1, dfz2, dfz3, dfz4;
            double z1, z2, z3, z4,error;

            zt = 0.6;
            Ppr = 5.60;
            Tpr = 1.35;
            error = 1;

            while (error > 0.0000001)
            {
                dr = (0.27 * Ppr) / (zt * Tpr);
                fz = zt - (1 + (A1 + A2 / Tpr + A3 / Math.Pow(Tpr, 3)) * dr + (A4 + A5 / Tpr) * Math.Pow(dr, 2) + (A5 * A6 * Math.Pow(dr, 5) / Tpr) + A7 * (1 + A8 * Math.Pow(dr, 2)) * (Math.Pow(dr, 2) / Math.Pow(Tpr, 3)) * Math.Exp(-A8 * Math.Pow(dr, 2)));

                dfz1 = (A1 + (A2 / Tpr) + (A3 / Math.Pow(Tpr, 3))) * (dr / zt);
                dfz2 = 2 * (A4 + (A5 / Tpr)) * (Math.Pow(Tpr, 2) / zt);
                dfz3 = (5 * A5 * A6 * Math.Pow(dr, 5)) / (zt * Tpr);
                dfz4 = ((2 * A7 * Math.Pow(dr, 2)) / (zt * Math.Pow(Tpr, 3))) * (1 + A8 * Math.Pow(dr, 2) - Math.Pow(A8 * Math.Pow(dr, 2), 2)) * Math.Exp(-A8 * Math.Pow(dr, 2));
                dfz = 1 + dfz1 + dfz2 + dfz3 + dfz4;

                ztdt = zt-(fz/dfz);

                error = Math.Abs(zt-ztdt);
                zt = ztdt;
            }

            DPRc = zt;
            Console.WriteLine(DPRc);

       }
       

References

Ahmed, Tarek. Equations of state and PVT analysis Applications for improved reservoir modeling
Bánzer, Carlos. Correlaciones numéricas PVT