//======================================================
// Thermodynamic Package in java
// Class Lee Kesler Properties of real gases
// Dr. Turhan Coban
// TUBÝTAK UME(National metrology Institute of Turkey)
// email : Turhan.Coban@ume.tubitak.gov.tr
// File Name   : LeeKesler.java
// This file contains the LeeKesler class
// this class sets basic properties of real gases
// required data is read from LeeKesler.txt
// =====================================================
// Description : This file contains the LeeKesler class
//               class LeeKesler calculates thermophysical properties of
//               real fluids
// DATA FILE DEFINATION
// datas are written in the data file "RealGas.txt"
//============================================================
//  VARIABLE IDENTIFICATION
//============================================================
//  GasEquation : String type variable
//                values : "LeeKesler"
//
//
//  Formula     : Fluid chemical formula
//  Name        : Fluid name
//  atomList    : Array of atom class, structral information
//  natom       : number of atoms
//  M           : Molecular weight as calculated from atomic structure
//                through class atom and atomList
//  MolWt       : Molecular Weigth from data file RealGas
//  Tfp         : Normal freezing point (1 atm), K
//  Tb          : Normal Boiling point  (1 atm), K
//  Tc          : Critical temperature, K
//  Pc          : Critical pressure, bar
//  Vc          : Critical volume, cm^3/mole
//  Zc          : Critical compresibility factor PcVc/RTc
//  Omega       : Pitzer's accentric factor
//  Dipm        : Dipole moment, debyes
//  CPC         : constants array to calculate the isobaric heat capacity of
//                the ideal gas Cp kJ/(kmolK)
//  DELHF       : formation enthalpy, 298.2 K, 1.01325 bar, KJ/kmol
//  DELGF       : formation gibbs free energy, 298.2 K, 1.01325 bar,KJ/kmol
//  DELS0       : entropy at 298.2 K, 1.01325 bar, KJ/kmol K
//  Vpeq        : vapor pressure equation number
//  VPC         : vapor pressure equation coefficients (array)
//  Tmin        : minimum twmp. K
//  Tmax        : maximum temp. K
//  LDEN        : liquid density, g/cm^3
//  TDEN        : temperature that liquid density is given, K
//============================================================
// Reference : The Properties of Gases & Liquids,
// Robert C. Reid, John M. Prausnitz, Bruce E. Pouling
// McGraw-Hill Book Company, ISBN 0-07-051799-1
//==============================================

import java.io.*;

//==============================================
// implementation of constructors functions
// for class LeeKesler (Equation of state)

public class LeeKesler
{
   String GasEquation;
   int natom;
   Atom atomList[];
   String Formula,Name;
   double M,MolWt,Tfp,Tb,Tc,Pc,Vc,Zc,Omega,Dipm;
   public double CPC[]=new double[4];
   public double VPC[]=new double[4];
   double DELHF,DELGF,DELS0,DELHO;
   int VPEq;
   double Tmin,Tmax,LDEN,TDEN;
   int ierror;
   BufferedReader fin;
   double R=8314.5; //(J/kmolK)
   //Peng Robinson constants
   double a,b,u,w;
   //
   //================================//
   // Equation of State coefficients //
   //================================//
   //Simple Fluid (omega=0) :
   //------------------------
   double b0[]=new double[4];
   double c0[]=new double[4];
   double d0[]=new double[2];
   double beta0;
   double gama0;
   double OmegaR;
   //Reference Fluid (n-octane):
   //------------------------
   double bR[]=new double[5];
   double cR[]=new double[5];
   double dR[]=new double[3];
   double betaR;
   double gamaR;
   double N;
   //constructor methods for Lee Kesler
   public LeeKesler()
   {
   // isobutane is assumed to be initial fluid
   try{ startup("isobutane",1.0); }
   catch(IOException e){}
   }

   public LeeKesler(String s2)
   {try{ startup(s2,1.0); }
   catch(IOException e){}
   }

   public LeeKesler(String s2,double Ne)
   {
   try{ startup(s2,Ne); }
   catch(IOException e){}
   }

   public LeeKesler(LeeKesler lk)
   {
   try{ startup(lk.Name,lk.N); }
   catch(IOException e){}
   }

   public void startup(String  s1,double Ne) throws IOException
   {
   GasEquation="Lee - Kesler";
   //setting the simple fluid coefficients :
   b0[0]=0.1181193;
   b0[1]=0.265728;
   b0[2]=0.15479;
   b0[3]=0.030323;
   c0[0]=0.0236744;
   c0[1]=0.0186984;
   c0[2]=0.0;
   c0[3]=0.042724;
   d0[0]=0.155488e-4;
   d0[1]=0.623689e-4;
   beta0=0.65392;
   gama0=0.060167;
   OmegaR=0.3978;
   //setting the reference fluid coefficients :
   bR[0]=0.2026579;
   bR[1]=0.331511;
   bR[2]=0.027655;
   bR[3]=0.203488;
   cR[0]=0.0313385;
   cR[1]=0.0503618;
   cR[2]=0.016901;
   cR[3]=0.041577;
   dR[0]=0.48736e-4;
   dR[1]=0.0740336e-4;
   betaR=1.226;
   gamaR=0.03754;
   int i;
   String atomName[];
   double atomN[];
   int anatom;
   ierror=1;
   String tempName="";
   String tempFormula="";
   String aFormula,aName;
   double aMolWt,aTfp,aTb,aTc,aPc,aVc,aZc,aOmega,aDipm ;
   double aCPC[]=new double[4];
   double aDELHF,aDELGF,aDELS0;
   int aVPEq;
   double aVPC[]=new double[4];
   double aTmin,aTmax,aLDEN,aTDEN ;
   try{
       fin=new BufferedReader(new FileReader("RealGas.txt"));
       } catch(IOException e)
         {
         System.err.println("Error Opening File RealGas.dat\n"+e.toString());
         System.exit(1);
         }
   try{
      while(fin!=null)
         {
         aFormula=Text.readString(fin);
         aName=Text.readString(fin);
         anatom=Text.readInt(fin);
         atomName=new String[anatom];
         atomN=new double[anatom];
         for(i=0;i<anatom;i++)
         {
         atomName[i]=Text.readString(fin);
         atomN[i]=Text.readDouble(fin);
         }
         aMolWt=Text.readDouble(fin);
         aTfp=Text.readDouble(fin);
         aTb=Text.readDouble(fin);
         aTc=Text.readDouble(fin);
         aPc=Text.readDouble(fin);
         aVc=Text.readDouble(fin);
         aZc=Text.readDouble(fin);
         aOmega=Text.readDouble(fin);
         aDipm=Text.readDouble(fin);
         for(int ii=0;ii< aCPC.length;ii++)
         { aCPC[ii]=Text.readDouble(fin); }
         aDELHF=Text.readDouble(fin);
         aDELGF=Text.readDouble(fin);
         aDELS0=Text.readDouble(fin);
         aVPEq=Text.readInt(fin);
         for(int ii=0;ii< aVPC.length;ii++)
         { aVPC[ii]=Text.readDouble(fin); }
         aTmin=Text.readDouble(fin);
         aTmax=Text.readDouble(fin);
         aLDEN=Text.readDouble(fin);
         aTDEN=Text.readDouble(fin);
         if (aName.equals(s1) || aFormula.equals(s1) )
         {
         Name=aName;
         Formula=aFormula;
         MolWt=aMolWt;
         Tfp=aTfp;
         Tb=aTb;
         Tc=aTc;
         Pc=aPc;
         Vc=aVc;
         Zc=aZc;
         Omega=aOmega;
         Dipm=aDipm;
         for(int ii=0;ii< CPC.length;ii++)
         { CPC[ii]=aCPC[ii]; }
         DELHF=aDELHF;
         DELGF=aDELGF;
         DELS0=aDELS0;
         VPEq=aVPEq;
         for(int ii=0;ii< VPC.length;ii++)
         { VPC[ii]=aVPC[ii]; }
         Tmin=aTmin;
         Tmax=aTmax;
         LDEN=aLDEN;
         TDEN=aTDEN;
         natom=anatom;
         N=Ne;
         atomList=new Atom[natom];
         M=0;
         for(i=0;i<natom;i++)
         {
         atomList[i]=new Atom(atomName[i],atomN[i]);
         M+=atomList[i].mass;
         }
         ierror=0;
         break;
         }
         }
      } catch(EOFException e_eof)
        {
        System.out.println("error required gas "+s1+" is not found");
        System.exit(1);
        fin.close();
        return;
        }
        fin.close();
   }

   public String toString()
   {
   //return the chemical symbol of the gas
   String s="";
   for(int i=0;i<natom;i++)
      s=s+atomList[i].toString();
   return s;
   }

   public double Z_TV_LK_V(double T, double V)
   {
   //Lee Kesler vapor compresibility factor
   double Tr=T/Tc;
   double Vr=Pc*1e5*V/(R*Tc);
   double B0,C0,D0;
   double BR,CR,DR;
   double Z0,ZR,ZX;
   B0=b0[0]-b0[1]/Tr-b0[2]/(Tr*Tr)-b0[3]/(Tr*Tr*Tr);
   C0=c0[0]-c0[1]/Tr+c0[2]/(Tr*Tr*Tr);
   D0=d0[0]+d0[1]/Tr;
   Z0=1+B0/Vr+C0/(Vr*Vr)+D0/(Vr*Vr*Vr*Vr*Vr)
   + c0[3]/(Tr*Tr*Tr*Vr*Vr)*(beta0+gama0/(Vr*Vr))*Math.exp(-gama0/(Vr*Vr));
   BR=bR[0]-bR[1]/Tr-bR[2]/(Tr*Tr)-bR[3]/(Tr*Tr*Tr);
   CR=cR[0]-cR[1]/Tr+cR[2]/(Tr*Tr*Tr);
   DR=dR[0]+dR[1]/Tr;
   ZR=1+BR/Vr+CR/(Vr*Vr)+DR/(Vr*Vr*Vr*Vr*Vr)
   + cR[3]/(Tr*Tr*Tr*Vr*Vr)*(betaR+gamaR/(Vr*Vr))*Math.exp(-gamaR/(Vr*Vr));
   ZX=Z0+Omega/OmegaR*(ZR-Z0);
   return ZX;
   }

   public void setconstants(double T)
   {
   //Peng-robinson cubic equation of state coefficients
     u=2;
     w=-1;
     double fw=0.37464+1.54226*Omega-0.26992*Omega*Omega;
     double Tr=T/Tc;
     b=0.07780*R*Tc/(Pc*1e5);
     a=0.45724*R*R*Tc*Tc/(Pc*1e5)*
     Math.pow(1+fw*(1.0-Math.pow(Tr,0.5)),2);
   }

   public double Z_TP_PR_V(double T, double P)
   {
     //This value is calculated by using Peng-Robinson Equation of state
     //vapor phase compresibility factor
     setconstants(T);
     double AY=a*(P*1e5)/(R*R*T*T);
     double BY=b*(P*1e5)/(R*T);
     double A=-(1+BY-u*BY);
     double B=AY+w*BY*BY-u*BY-u*BY*BY;
     double C=-AY*BY-w*BY*BY-w*BY*BY*BY;
     double cof[]=new double[4];
     cof[0]=1;
     cof[1]=A;
     cof[2]=B;
     cof[3]=C;
     complex zc[]=Matrix.poly_rootsC(cof);
     //complex zc[]=Matrix.cubic_rootsC(cof);
     //System.out.println(Matrix.toString(zc));
     return max(zc);
   }

   public double Z_TP_PR_L(double T, double P)
   {
     // This value is calculated by using Peng-Robinson Equation of state
     // liquid phase compressibility factor
     // if state is above critical liquid compressibility factor will be same as
     // vapor compressibility factor.
        setconstants(T);
        double AY=a*(P*1e5)/(R*R*T*T);
        double BY=b*(P*1e5)/(R*T);
        double A=-(1+BY-u*BY);
        double B=AY+w*BY*BY-u*BY-u*BY*BY;
        double C=-AY*BY-w*BY*BY-w*BY*BY*BY;
        double cof[]=new double[4];
        cof[0]=1;
        cof[1]=A;
        cof[2]=B;
        cof[3]=C;
       complex zc[]=Matrix.poly_rootsC(cof);
       //complex zc[]=Matrix.cubic_rootsC(cof);
       double zz=min(zc);
       return zz;
   }


   double max(complex x[])
   {
     //maximum of complex compressibility factor solution (Peng Robinson E.S.)
     // this solution also check if solution is real or complex
     // and ignore complex solution points
     double nmax=-1e99;
     for(int i=0;i<x.length;i++)
     {
     if((x[i].nreal>=0.0 && x[i].nreal<=1.0) &&
       (x[i].nimag>=0.0 && x[i].nimag<1e-4))
       {
       if(x[i].nreal>nmax) nmax=x[i].nreal;
       }
     }
   return nmax;
   }

   double min(complex x[])
   {
   //minimum of complex compressibility factor solution (Peng Robinson E.S.)
   // this solution also chack if solution is real or complex
   // and ignore complex solution points
     double nmin=1e99;
     for(int i=0;i<x.length;i++)
     {
     if((x[i].nreal>=0.0 && x[i].nreal<=1.0) &&
       (x[i].nimag>=0.0 && x[i].nimag<1e-4))
       {
       if(x[i].nreal<nmin) nmin=x[i].nreal;
       }
     }
   return nmin;
   }

   public double Z_TP_LK_V(double T, double P)
   {
   // Vapor compressibility factor by using Lee-Kesler equation of state
   double Tr=T/Tc;
   double ZPR=Z_TP_PR_V(T,P);
   //System.out.println("ZPR = "+ZPR);
   double Pr=P/Pc;
   double V=ZPR*R*T/(P*1e5);
   //System.out.println("V = "+V);
   // First volume guess assume Z value from Peng Robinson equation of State
   double Vr=Pc*1e5*V/(R*Tc);
   double B0,C0,D0;
   double BR,CR,DR;
   double Z0,ZR;
   double ZX=0;
   B0=b0[0]-b0[1]/Tr-b0[2]/(Tr*Tr)-b0[3]/(Tr*Tr*Tr);
   C0=c0[0]-c0[1]/Tr+c0[2]/(Tr*Tr*Tr);
   D0=d0[0]+d0[1]/Tr;
   BR=bR[0]-bR[1]/Tr-bR[2]/(Tr*Tr)-bR[3]/(Tr*Tr*Tr);
   CR=cR[0]-cR[1]/Tr+cR[2]/(Tr*Tr*Tr);
   DR=dR[0]+dR[1]/Tr;
   double f0,df0,fR,dfR,f,df;
   int nmax=150;
   double tolerance=1.0e-15;
   int i=0;
   while(i++<nmax)
   {
   f0=1+B0/Vr+C0/(Vr*Vr)+D0/(Vr*Vr*Vr*Vr*Vr)
   + c0[3]/(Tr*Tr*Tr*Vr*Vr)*(beta0+gama0/(Vr*Vr))*Math.exp(-gama0/(Vr*Vr));
   fR=1+BR/Vr+CR/(Vr*Vr)+DR/(Vr*Vr*Vr*Vr*Vr)
   + cR[3]/(Tr*Tr*Tr*Vr*Vr)*(betaR+gamaR/(Vr*Vr))*Math.exp(-gamaR/(Vr*Vr));
   f=f0+Omega/OmegaR*(fR-f0)-Pr*Vr/Tr;
   df0=-B0/(Vr*Vr)-2*C0/(Vr*Vr*Vr)-5*D0/(Vr*Vr*Vr*Vr*Vr*Vr)+
   2.0*c0[3]/(Tr*Tr*Tr*Vr*Vr*Vr)*Math.exp(-gama0/(Vr*Vr))*
   (-beta0-2*gama0/(Vr*Vr)+gama0*beta0/(Vr*Vr)+gama0*gama0/(Vr*Vr*Vr*Vr));
   dfR=-BR/(Vr*Vr)-2*CR/(Vr*Vr*Vr)-5*DR/(Vr*Vr*Vr*Vr*Vr*Vr)+
   2.0*cR[3]/(Tr*Tr*Tr*Vr*Vr*Vr)*Math.exp(-gamaR/(Vr*Vr))*
   (-betaR-2*gamaR/(Vr*Vr)+gamaR*betaR/(Vr*Vr)+gamaR*gamaR/(Vr*Vr*Vr*Vr));
   df=df0+Omega/OmegaR*(dfR-df0)-Pr/Tr;
   Vr-=f/df;
   ZX=Pr*Vr/Tr;
   //System.out.println("VR= "+Vr+"ZX="+ZX);
   if(Math.abs(f)<tolerance)
      {break;}
   }
   if(ZX<0) return ZPR;
   else return ZX;
   }

   public double Z_TP_LK_L(double T, double P)
   {
   // Liquid compressibility factor by using Lee-Kesler equation of state
   double Tr=T/Tc;
   double Pr=P/Pc;
   double ZPR=Z_TP_PR_L(T,P);
   double V=ZPR*R*T/(P*1e5);
   {
   // First guess assume Z value from Peng Robinson equation of State
   double Vr=Pc*1e5*V/(R*Tc);
   double B0,C0,D0;
   double BR,CR,DR;
   double Z0,ZR,ZX;
   B0=b0[0]-b0[1]/Tr-b0[2]/(Tr*Tr)-b0[3]/(Tr*Tr*Tr);
   C0=c0[0]-c0[1]/Tr+c0[2]/(Tr*Tr*Tr);
   D0=d0[0]+d0[1]/Tr;
   BR=bR[0]-bR[1]/Tr-bR[2]/(Tr*Tr)-bR[3]/(Tr*Tr*Tr);
   CR=cR[0]-cR[1]/Tr+cR[2]/(Tr*Tr*Tr);
   DR=dR[0]+dR[1]/Tr;
   double f0,df0,fR,dfR,f,df;
   int nmax=150;
   double tolerance=1.0e-15;
   int i=0;
   while(i++<nmax)
   {
   f0=1+B0/Vr+C0/(Vr*Vr)+D0/(Vr*Vr*Vr*Vr*Vr)
   + c0[3]/(Tr*Tr*Tr*Vr*Vr)*(beta0+gama0/(Vr*Vr))*Math.exp(-gama0/(Vr*Vr));
   fR=1+BR/Vr+CR/(Vr*Vr)+DR/(Vr*Vr*Vr*Vr*Vr)
   + cR[3]/(Tr*Tr*Tr*Vr*Vr)*(betaR+gamaR/(Vr*Vr))*Math.exp(-gamaR/(Vr*Vr));
   f=f0+Omega/OmegaR*(fR-f0)-Pr*Vr/Tr;
   df0=-B0/(Vr*Vr)-2*C0/(Vr*Vr*Vr)-5*D0/(Vr*Vr*Vr*Vr*Vr*Vr)+
   2.0*c0[3]/(Tr*Tr*Tr*Vr*Vr*Vr)*Math.exp(-gama0/(Vr*Vr))*
   (-beta0-2*gama0/(Vr*Vr)+gama0*beta0/(Vr*Vr)+gama0*gama0/(Vr*Vr*Vr*Vr));
   dfR=-BR/(Vr*Vr)-2*CR/(Vr*Vr*Vr)-5*DR/(Vr*Vr*Vr*Vr*Vr*Vr)+
   2.0*cR[3]/(Tr*Tr*Tr*Vr*Vr*Vr)*Math.exp(-gamaR/(Vr*Vr))*
   (-betaR-2*gamaR/(Vr*Vr)+gamaR*betaR/(Vr*Vr)+gamaR*gamaR/(Vr*Vr*Vr*Vr));
   df=df0+Omega/OmegaR*(dfR-df0)-Pr/Tr;
   Vr-=f/df;
   if(Math.abs(f0)<tolerance)
      {break;}
   }
   ZX=Pr*Vr/Tr;
   if(ZX>0)
      return ZX;
   else
      return ZPR;
   }
   }

   double dZ_TP_LK_VL(double T,double P)
   {
   return Z_TP_LK_V(T,P)-Z_TP_LK_L(T,P);
   }

   double P_TV_LK_V(double T, double V)
   {
   return Z_TV_LK_V(T,V)*R*T/V*1e-5;
   }

   double V_TP_LK_V(double T, double P)
   {
   return Z_TP_LK_V(T,P)*R*T/(P*1e5);
   }

   double V_TP_LK_L(double T, double P)
   {
   return Z_TP_LK_L(T,P)*R*T/(P*1e5);
   }

   double V_T_LK_SL(double T)
   {
   if(T<Tc)
   {
   return Z_TP_LK_L(T,P_T_SV(T))*8314.5*T/(P_T_SV(T)*1e5);
   }
   else return 0;
   }

   double V_T_LK_SV(double T)
   {
   if(T<Tc)
   {
   return Z_TP_LK_V(T,P_T_SV(T))*8314.5*T/(P_T_SV(T)*1e5);
   }
   else
   return 0;
   }

   double V_T_BE_SL(double T)
   {
   //Volume of saturated liquid as a function of temperature
   // Bhirud equation
   if(T<Tc)
   {
   double Tr=T/Tc;
   double Tr2=Tr*Tr;
   double Tr3=Tr*Tr2;
   double Tr4=Tr2*Tr2;
   double Tr5=Tr3*Tr2;
   double Tr6=Tr3*Tr3;
   double lnV0=1.39644-24.076*Tr+102.615*Tr2-255.719*Tr3+
   355.805*Tr4-256.671*Tr5+75.1088*Tr6;
   double lnV1=13.4412-135.7437*Tr+533.38*Tr2-1091.453*Tr3+
   1231.43*Tr4-728.227*Tr5+176.737*Tr6;
   return (R*T/(Pc*1e5)*Math.exp(lnV0+Omega*lnV1));
   }
   else
   {
   return 0;
   }
   }

   double V_T_RE_SL(double T)
   {
   //Volume of saturated liquid as a function of temperature
   // Rackett equation
   if(T<Tc)
   {
   double deger=0;
   double Tr=T/Tc;
   double ZRA;
      if(0.98>Tr && Tr>=1.0 )
        {ZRA=Zc;}
      else
        {ZRA=0.29506-0.08775*Omega;}
     if(LDEN!=0)
     {
      double VSR=1.0/(LDEN*1e3/M);
      double TrR=TDEN/Tc;
      double SI=Math.pow((1-Tr),(2.0/7.0))-Math.pow((1-TrR),(2.0/7.0));
      ZRA=0.29506-0.08775*Omega;
      deger=VSR*Math.pow(ZRA,SI);
     }
     else
     {
     double SI=1+Math.pow((1-Tr),(2.0/7.0));
     deger=R*Tc/Pc*Math.pow(ZRA,SI);
     }
   return deger;
   }
   else
   return 0;
   }

   double V_P_RE_SL(double P)
   {
   //Volume of saturated liquid as a function of pressure
   // Rackett equation
   if(P<Pc)
   {
   return V_T_RE_SL(T_P_SV(P));
   }
   else
   return 0;
   }

   double V_P_BE_SL(double P)
   {
   //Volume of saturated liquid as a function of pressure
   if(P<Pc)
   {
   return V_T_BE_SL(T_P_SV(P));
   }
   else
   return 0;
   }

   public double P_T_SV(double T)
   {
   //boiling (saturated vapor) pressure (bara) of the substance
   // at a given temperature
   if(T<Tc)
   {
   double deger=0;
   if(VPEq==1)
      {
       double x=1-T/Tc;
       double x3=x*x*x;
       double x6=x3*x3;
       deger=Pc*Math.exp((VPC[0]*x+VPC[1]*Math.pow(x,1.5)+VPC[2]*x3+VPC[3]*x6)/(1-x));
      }
   else if(VPEq==2)
      {
      //first guess is done by using Gomez-Thodos method
      double VP=P_T_GT_SV(T);
      double fx,dfx;
      int nmax=100;
      double tolerance=1.0e-10;
      int i=0;
      while(i++<nmax)
      {
        fx = -Math.log(VP)+(VPC[0]-VPC[1]/T+VPC[2]*Math.log(T)+VPC[3]/(T*T)*VP);
        dfx=VPC[3]/(T*T)-1.0/VP;
        VP-=fx/dfx;
        if(Math.abs(fx)<tolerance)
          {deger=VP;break;}
      }
      deger=VP;
      }
   else if(VPEq==3)
      {
      deger=Math.exp(VPC[0]-VPC[1]/(T+VPC[2]));
      }
   else if(VPEq==4)
      {
      //assume Lee - Kesler estimation
      deger=P_T_LK_SV(T);
      }
   return deger;
   }
   else
   return 0;
   }//end of PVp

   public double T_P_SV(double P)
   {
   //boiling Temperature (K) of the substance
   // at a given pressure
   {
   double deger=0;
   int nmax=150;
   double tolerance=1.0e-10;
   double Pr=P/Pc;
   //Lee Kesler equation is used as first estimation
   double T=T_P_LK_SV(P);
   if(VPEq==1)
      {
       double x=1-T/Tc;
       double fx,dfx;
      int i=0;
      while(i++<nmax)
      {
      double x3=x*x*x;
      double x6=x3*x3;
        double fx1 = (VPC[0]*x+VPC[1]*Math.pow(x,1.5)+VPC[2]*x3+VPC[3]*x6)/(1-x);
        fx=fx1-Math.log(Pr);
        dfx=(VPC[0]+1.5*VPC[1]*Math.pow(x,0.5)+3*x*x*VPC[2]+6*VPC[3]*x6)/x +
        (VPC[0]*x+VPC[1]*Math.pow(x,1.5)+VPC[2]*x3+VPC[3]*x6)/(x*x);
        x-=fx/dfx;
        T=(1-x)*Tc;
        if(Math.abs(fx)<tolerance)
          {break;}
      }
      deger=T;
      }
   else if(VPEq==2)
      {
      double x=1-T/Tc;
      double fx,dfx;
      int i=0;
      while(i++<nmax)
      {
        fx = VPC[0]-VPC[1]/T+VPC[2]*Math.log(T)+VPC[3]*P/(T*T)-Math.log(P);
        dfx= VPC[1]/(T*T)+VPC[2]/T-2.0*VPC[3]*P/(T*T*T);
        T-=fx/dfx;
         if(Math.abs(fx)<tolerance)
          {break;}
      }
      deger=T;
      }
   else if(VPEq==3)
      {
      deger=VPC[1]/(VPC[0]-Math.log(Pr))-VPC[2];
      }
   else if(VPEq==4 && VPEq==0)
      {
       deger=T_P_LK_SV(P);
      }
   return deger;
   }
   }//end of

   public double P_T_LK_SV(double T)
   {
   //Lee-Kesler saturated vapor pressure estimation
   double Tr=T/Tc;
   double Tr6=Tr*Tr*Tr*Tr*Tr*Tr;
   double f1=5.92714-6.09648/Tr-1.28862*Math.log(Tr)+0.169347*Tr6;
   double f2=15.2518-15.6875/Tr-13.4721*Math.log(Tr)+0.43577*Tr6;
   double teta=Tb/Tc;
   double teta6=teta*teta*teta*teta*teta*teta;
   double alfa=-Math.log(Pc)-5.92714+6.09648/teta+1.28862*Math.log(teta)-
   0.169347*teta6;
   double beta=15.2518-15.6875/teta-13.4721*Math.log(teta)+0.43577*teta6;
   double ww=alfa/beta;
   return Pc*Math.exp(f1+ww*f2);
   }

   public double T_P_LK_SV(double P)
   {
   //Lee-Kesler vapor temperature estimation
   //boiling Temperature (K) of the substance
   // at a given pressure
   double teta=Tb/Tc;
   double teta6=teta*teta*teta*teta*teta*teta;
   double alfa=-Math.log(Pc)-5.92714+6.09648/teta+1.28862*Math.log(teta)-
   0.169347*teta6;
   double beta=15.2518-15.6875/teta-13.4721*Math.log(teta)+0.43577*teta6;
   double ww=alfa/beta;
      //first guess is carried out by using gomez - thodos equation
      double T=T_P_GT_SV(P);
      double Tr=T/Tc;
      double Tr6=Tr*Tr*Tr*Tr*Tr*Tr;
      double Pr=P/Pc;
      double fx,dfx;
      int nmax=500;
      double tolerance=1.0e-15;
      int i=0;
      while(i++<nmax)
      {
      double f1=5.92714-6.09648/Tr-1.28862*Math.log(Tr)+0.169347*Tr6;
      double f2=15.2518-15.6875/Tr-13.4721*Math.log(Tr)+0.43577*Tr6;
      double df1=6.09648/(Tr*Tr)-1.28862/Tr-1.1016082*Tr6/Tr;
      double df2=15.6875/(Tr*Tr)-13.4721/Tr+2.61462*Tr6/Tr;
        fx = f1+ww*f2-Math.log(Pr);
        dfx=df1+ww*df2;
        Tr-=fx/dfx;
        if(Math.abs(fx)<tolerance)
          {break;}
      }
      T=Tr*Tc;
     return T;
   }

   public double P_T_WE_SV(double T)
   {
   // Lee Kesler - Wagner equation for vapor pressure estimation
     if(T<Tc)
     {
     double A[][]=
    {{369.85, 42.4535, 0.153, -6.72219, 1.33236, -2.13868, -1.38551},
     {562.16, 48.9794, 0.212, -6.98273, 1.33213, -2.62863, -3.33399},
     {568.81, 24.8617, 0.398, -7.91211, 1.38007, -3.80435, -4.50132},
     {566.52, 31.2481, 0.415, -8.05688, 1.46673, -3.82439, -2.78724}};
     double deger=0;
     double Tr=T/Tc;
     double x=1-Tr;
     double x3=x*x*x;
     double x6=x3*x3;
     for(int i=0;i<3;i++)
     {
      if(Omega>A[i][2]&&Omega<A[i+1][2])
      {
       double lnPvpR1=(A[i][3]*x+A[i][4]*Math.pow(x,1.5)+A[i][5]*x3+A[i][6]*x6)/(1-x);
       double lnPvpR2=(A[i+1][3]*x+A[i+1][4]*Math.pow(x,1.5)+A[i+1][5]*x3+A[i+1][6]*x6)/(1-x);
       double y1=(Omega-A[i][2])/(A[i+1][2]-A[i][2]);
       double y=lnPvpR1+(lnPvpR2-lnPvpR1)*y1;
       deger=Pc*Math.exp(y);
       break;
      }
     }
     return deger;
     }
     return 0;
   }

   public double T_P_WE_SV(double P)
   {
   if(P<Pc)
   {
   //Lee Kesler - Wagner equation for vapor temperature estimation
   //boiling Temperature (K) of the substance at a given pressure
     double A[][]=
    {{369.85, 42.4535, 0.153, -6.72219, 1.33236, -2.13868, -1.38551},
     {562.16, 48.9794, 0.212, -6.98273, 1.33213, -2.62863, -3.33399},
     {568.81, 24.8617, 0.398, -7.91211, 1.38007, -3.80435, -4.50132},
     {566.52, 31.2481, 0.415, -8.05688, 1.46673, -3.82439, -2.78724}};
      //first guess
      double T=T_P_SV(P);
      double Tr=T/Tc;
      double x=1.0-Tr;
      double Pr=P/Pc;
      int nmax=200;
      double tolerance=1.0e-20;
      int j=0;
      for(int i=0;i<3;i++)
      {
      if(Omega>A[i][2]&&Omega<A[i+1][2])
      {
      while(j++<nmax)
      {
       double fx,dfx;
       double x2=x*x;
       double x3=x*x*x;
       double x5=x3*x2;
       double x6=x3*x3;
       double lnPvpR1=(A[i][3]*x+A[i][4]*Math.pow(x,1.5)+A[i][5]*x3+A[i][6]*x6)/(1-x);
       double lnPvpR2=(A[i+1][3]*x+A[i+1][4]*Math.pow(x,1.5)+A[i+1][5]*x3+A[i+1][6]*x6)/(1-x);
       double d_lnPvpR1=(A[i][3]+1.5*A[i][4]*Math.pow(x,0.5)+3*A[i][5]*x2+
       6*A[i][6]*x5)/(1-x)+lnPvpR1/(1-x);
       double d_lnPvpR2=(A[i+1][3]+1.5*A[i+1][4]*Math.pow(x,0.5)+3*A[i+1][5]*x2+
       6*A[i+1][6]*x5)/(1-x)+lnPvpR2/(1-x);
       double y1=(Omega-A[i][2])/(A[i+1][2]-A[i][2]);
       fx=lnPvpR1+(lnPvpR2-lnPvpR1)*y1-Math.log(Pr);
       dfx=d_lnPvpR1+(d_lnPvpR2-d_lnPvpR1)*y1;
       //fx=Pc*Math.exp(y);
       x-=fx/dfx;
        if(Math.abs(fx)<tolerance)
          {break;}
      } //end of while
      } //end of if
      }//end of for

      T=Tc*(1-x);
      return T;
      }
      else
        return 0;
   }


   public double P_T_GT_SV(double T)
   {
   if(T<Tc)
   {
   //Gomez-Thodos vapor pressure estimation
   double Tbr=Tb/Tc;
   double Tr=T/Tc;
   double h=Tbr*Math.log(Pc)/(1-Tbr);
   double delta=0;
   double m,gama,beta,aa,bb;
   double Tbr7=Tbr*Tbr*Tbr*Tbr*Tbr*Tbr*Tbr;
   if(Formula=="H2O" &&  Name=="Water")
     {
      delta=0.19904;
     }
   if(Formula=="H2" && Name=="Hydrogen")
     {
     delta=0.41815;
     }
   if(Formula=="Ne")
     {
     delta=0.02319;
     }
   m=0.78425*Math.exp(0.089315*h)-8.5217/Math.exp(0.74826*h);
   beta=-4.267-221.79/(Math.pow(h,2.5)*Math.exp(0.0384*Math.pow(h,2.5)))+
   3.8126/Math.exp(2272.44/(h*h*h))+delta;
   aa=(1.0-1.0/Tbr)/(Tbr7-1);
   bb=(1.0-1/Math.pow(Tbr,m))/(Tbr7-1);
   gama=aa*h+bb*beta;
   if(Formula=="H2" && Name=="Hydrogen")
     {
     m=0.0052*Math.pow(M,0.29)*Math.pow(Tc,0.72);
     gama=2.464/M*Math.exp(9.8e-6*M*Tc);
     beta=gama/bb-aa*h/bb;
     }
   double sonuc=Pc*Math.exp(beta*(1.0/Math.pow(Tr,m)-1)+
   gama*(Tr*Tr*Tr*Tr*Tr*Tr*Tr-1.0));
   return sonuc;
   }
   else
    return 0;
   }

   public double T_P_GT_SV(double P)
   {
   if(P>Pc)
      return 0;
   else
   {
   //Gomez-Thodos vapor Temperature estimation
   //boiling Temperature (K) of the substance
   // at a given pressure
   double Tbr=Tb/Tc;
   double h=Tbr*Math.log(Pc)/(1-Tbr);
   double delta=0;
   double m,gama,beta,aa,bb;
   double Tbr7=Tbr*Tbr*Tbr*Tbr*Tbr*Tbr*Tbr;
   if(Formula=="H2O" &&  Name=="Water")
     {
      delta=0.19904;
     }
   if(Formula=="H2" && Name=="Hydrogen")
     {
     delta=0.41815;
     }
   if(Formula=="Ne")
     {
     delta=0.02319;
     }
   m=0.78425*Math.exp(0.089315*h)-8.5217/Math.exp(0.74826*h);
   beta=-4.267-221.79/(Math.pow(h,2.5)*Math.exp(0.0384*Math.pow(h,2.5)))+
   3.8126/Math.exp(2272.44/(h*h*h))+delta;
   aa=(1.0-1.0/Tbr)/(Tbr7-1);
   bb=(1.0-1/Math.pow(Tbr,m))/(Tbr7-1);
   gama=aa*h+bb*beta;
   if(Formula=="H2" && Name=="Hydrogen")
     {
     m=0.0052*Math.pow(M,0.29)*Math.pow(Tc,0.72);
     gama=2.464/M*Math.exp(9.8e-6*M*Tc);
     beta=gama/bb-aa*h/bb;
     }
      //first guess is carried out by using clapeyron equation
      double T=T_P_CE_SV(P);
      double Tr=T/Tc;
      double Pr=P/Pc;
      double fx,dfx;
      int nmax=500;
      double tolerance=1.0e-15;
      int i=0;
      while(i++<nmax)
      {
        fx = -Math.log(Pr)+beta*(1.0/Math.pow(Tr,m)-1)+
        gama*(Tr*Tr*Tr*Tr*Tr*Tr*Tr-1.0);
        dfx=beta*(-m/Math.pow(Tr,(m+1))-1)+
        7.0*gama*Tr*Tr*Tr*Tr*Tr*Tr;
        Tr-=fx/dfx;
        if(Math.abs(fx)<tolerance)
          {break;}
      }
      T=Tr*Tc;
     return T;
   }
   }

   public double P_T_CE_SV(double T)
   {
   //Clapeyron equation vapor pressure
   //boiling Temperature (K) of the substance
   // at a given pressure
   if(T>Tc) return 0;
   else
   {
   double Tbr=Tb/Tc;
   double h=Tbr*Math.log(Pc)/(1-Tbr);
   double Tr=T/Tc;
   return Pc*Math.exp(h*(1-1/Tr));
   }
   }

   public double T_P_CE_SV(double P)
   {
   //Clapeyron equation vapor temperature
   if(P>Pc) return 0;
   else
   {
   double Tbr=Tb/Tc;
   double h=Tbr*Math.log(Pc)/(1-Tbr);
   double Pr=P/Pc;
   return Tc/(1-1/h*Math.log(Pr));
   }
   }

   public double dH_T_TS_SVL(double T)
   {
   //Torquato and Stell equation
   // as a function of temperature
   // Torquato S.and G. R. Stell; Ind. Eng. Chem. Fundam., 21:202(1982)
    //Reference fluids benzene and carbazole
   if(T>Tc) return 0;
   else
   {
   double To=1-T/Tc;
   double dHR1=6.537*Math.pow(To,(1.0/3.0))-2.467*Math.pow(To,(5.0/6.0))-
   77.521*Math.pow(To,1.208)+59.634*To+36.009*To*To-14.606*To*To*To;
   double dHR21=-0.133*Math.pow(To,(1.0/3.0))-28.215*Math.pow(To,(5.0/6.0))-
   82.958*Math.pow(To,1.208)+99.0*To+19.105*To*To-2.796*To*To*To;
   double dH=dHR1+4*(Omega-0.21)*dHR21;
   double dH1=dH*R*Tc*1e-3;
   //unit of result KJ/kmol
   return dH1;
   }
   }

    public double dH_P_TS_SVL(double P)
    {
    //Torquato and Stell equation
    // as a function of Pressure
    // Torquato S.and G. R. Stell; Ind. Eng. Chem. Fundam., 21:202(1982)
    //Reference fluids benzene and carbazole
    if(P<Pc)
      return dH_T_TS_SVL(T_P_GT_SV(P));
    else
      return 0;
    }

    public double dH_P_VE_SVL(double P)
    {
    if(P<Pc)
      return dH_T_VE_SVL(T_P_GT_SV(P));
    else
      return 0;
    }

    public double dH_P_LK_SVL(double P)
    {
    if(P<Pc)
      return dH_T_LK_SVL(T_P_GT_SV(P));
    else
      return 0;
    }

    public double dH_P_WE_SVL(double P)
    {
    if(P<Pc)
      return dH_T_WE_SVL(T_P_GT_SV(P));
    else
      return 0;
    }

   public double dH_T_LK_SVL(double T)
   {
    //Lee-kesler equation for enthalpy of evaporation
    if(T>Tc) return 0;
    else
    {
    double Tr=T/Tc;
    double deger=6.09648-1.28862*Tr+1.016*Tr*Tr+
    Omega*(15.6875-13.4721*Tr+2.615*Tr*Tr*Tr*Tr*Tr*Tr*Tr);
    deger*=R*Tc*dZ_TP_LK_VL(T,P_T_SV(T))*1e-3;
    return deger;
    }
   }

    public double dH_P_GT_SVL(double P)
    {
    if(P<Pc)
      return dH_T_GT_SVL(T_P_GT_SV(P));
    else
      return 0;
    }
    
   public double dH_T_GT_SVL(double T)
   {
    //Ghomez-thodos equation for enthalpy of evaporation
   if(T>Tc) return 0;
   else
   {
   double Tbr=Tb/Tc;
   double Tr=T/Tc;
   double h=Tbr*Math.log(Pc)/(1-Tbr);
   double delta=0;
   double m,gama,beta,aa,bb;
   double Tbr7=Tbr*Tbr*Tbr*Tbr*Tbr*Tbr*Tbr;
   if(Formula=="H2O" &&  Name=="Water")
     {
      delta=0.19904;
     }
   if(Formula=="H2" && Name=="Hydrogen")
     {
     delta=0.41815;
     }
   if(Formula=="Ne")
     {
     delta=0.02319;
     }
   m=0.78425*Math.exp(0.089315*h)-8.5217/Math.exp(0.74826*h);
   beta=-4.267-221.79/(Math.pow(h,2.5)*Math.exp(0.0384*Math.pow(h,2.5)))+
   3.8126/Math.exp(2272.44/(h*h*h))+delta;
   aa=(1.0-1.0/Tbr)/(Tbr7-1);
   bb=(1.0-1/Math.pow(Tbr,m))/(Tbr7-1);
   gama=aa*h+bb*beta;
   if(Formula=="H2" && Name=="Hydrogen")
     {
     m=0.0052*Math.pow(M,0.29)*Math.pow(Tc,0.72);
     gama=2.464/M*Math.exp(9.8e-6*M*Tc);
     beta=gama/bb-aa*h/bb;
     }
   double sonuc=7.0*gama*Tr*Tr*Tr*Tr*Tr*Tr-beta*m/Math.pow(Tr,(m-1));
   sonuc*=R*Tc*dZ_TP_LK_VL(T,P_T_SV(T))*1e-3;
   return sonuc;
   }
   }

   public double dH_T_WE_SVL(double T)
   {
   // Lee Kesler - Wagner equation for enthapy of vaporisation
     double A[][]=
    {{369.85, 42.4535, 0.153, -6.72219, 1.33236, -2.13868, -1.38551},
     {562.16, 48.9794, 0.212, -6.98273, 1.33213, -2.62863, -3.33399},
     {568.81, 24.8617, 0.398, -7.91211, 1.38007, -3.80435, -4.50132},
     {566.52, 31.2481, 0.415, -8.05688, 1.46673, -3.82439, -2.78724}};
     double deger=0;
     double Tr=T/Tc;
     double x=1-Tr;
     double x2=x*x;
     double x5=x2*x2*x;
     for(int i=0;i<3;i++)
     {
      if(Omega>A[i][2]&&Omega<A[i+1][2])
      {
      double f1=-A[i][3]+A[i][4]*Math.pow(x,0.5)*(0.5*x-1.5)+
       A[i][5]*x2*(2.0*x-3.0)+A[i][6]*x5*(5.0*x-6.0);
       double f2=-A[i+1][3]+A[i+1][4]*Math.pow(x,0.5)*(0.5*x-1.5)+
       A[i+1][5]*x2*(2.0*x-3.0)+A[i+1][6]*x5*(5.0*x-6.0);
       double y1=(Omega-A[i][2])/(A[i+1][2]-A[i][2]);
       double y=f1+(f2-f1)*y1;
       deger=y*R*Tc*dZ_TP_LK_VL(T,P_T_SV(T))*1e-3;
       break;
      }
     }
     return deger;
   }

   public double dHb_T_VE_SVL()
   {
     double Tbr=Tb/Tc;
     double dH=R*1e-3*Tc*Tbr*(0.4343*Math.log(Pc)-0.69431+0.89584*Tbr)/
     (0.37691-0.37306*Tbr+0.15075/Pc/(Tbr*Tbr));
     return dH;
   }

   public double dH_T_VE_SVL(double T)
   {
     if(T>Tc) return 0;
     else
     {
     double Tr=T/Tc;
     double n=(0.00264*dHb_T_VE_SVL()/(R*Tb)+0.8794);
     n=n*n*n*n*n*n*n*n*n*n;
     double dH=dHb_T_VE_SVL()*Math.pow((1-Tr)/(1-Tb/Tc),n);
     return dH;
     }
   }


//===================================================================
// Ideal gas equations
public double Cp_T_IG_V(double T)
{
// Specific heat under constant pressure, Cp of ideal gas
  double Cpi=0;
  Cpi=CPC[CPC.length-1];
  for(int i=CPC.length-2;i>=0;i--)
   { Cpi=Cpi*T+CPC[i]; }
  return Cpi;
}

public double h_T_IG_V(double T)
{
//enthalpy of ideal gas
  double T2=T*T,T02=298.15*298.15;
  double T3=T2*T,T03=298.15*T02;
  double T4=T2*T2,T04=T02*T02;
  double h=CPC[0]*(T-298.15)+CPC[1]*(T2-T02)/2.0+CPC[2]*(T3-T03)/3.0+
  CPC[3]*(T4-T04)/4.0;
  return h;
}

public double hf_T_IG_V(double T)
{
    return DELHF+h_T_IG_V(T);
}


public double s_TP_IG_V(double T,double P)
{
double s=DELS0+CPC[0]*Math.log(T/298.15)+CPC[1]*(T-298.15)+
CPC[2]/2.0*(T*T-298.15*298.15)+
CPC[3]/3.0*(T*T*T-298.15*298.15*298.15)+
-R*Math.log(P/1.01325);
return s;
}

//===================================================================
//departure functions of the real fluids
//
public double dCp_T_BR_L(double T)
{
  //liquid specific heat difference compare to IG vapor specific heat
  if(T<Tc)
  {
  double Tr=T/Tc;
  double x=1-Tr;
  double sonuc=1.45+0.45/x+0.25*Omega*(17.11+25.2*Math.pow(x,(1.0/3.0))/Tr+1.742/x);
  sonuc*=R*1e-3;
  return sonuc;
  }
  else
   return 0;
}

public double dCs_T_BR_L(double T)
{
// saturated liquid specific heat difference compare to IG vapor specific heat

  if(T<Tc)
  {
  double Tr=T/Tc;
  double x=1-Tr;
  double sonuc=dCp_T_BR_L(T)-Math.exp(20.1*Tr-17.9)*R*1e-3;
  return sonuc;
  }
  else
    return 0;
}

public double Cp_T_BR_L(double T)
{
//liquid specific heats
return dCp_T_BR_L(T)+Cp_T_IG_V(T);
}

public double Cp_T_BR_SL(double T)
{
//saturated liquid specific heats
return dCs_T_BR_L(T)+Cp_T_IG_V(T);
}

public double Cp_TP_BR_L(double T,double P)
{
double dT=0.01;	
// liquid enthalpy
return (h_TP_BR_L((T+dT),P)-h_TP_BR_L(T,P))/dT;
}

public double h_T_BR_SL(double T)
{
//enthalpy of the saturated liquid
return Cp_T_BR_SL(T)*T;
}



public double h_P_BR_SL(double P)
{
//enthapy of the saturated liquid
double T=T_P_LK_SV(P);
return Cp_T_BR_SL(T)*T;
}

public double dh_T_BR_L(double T)
{

  double x=1-T/Tc;
  double x0=1-298.15/Tc;
  double sonuc1=-R*Tc*1e-3*(1.45*x+0.45*Math.log(x)+0.25*Omega*(17.11*x+
  25.2*(Math.pow(x,(-2/3))/3.0*(1-x)+Math.pow(x,(1/3)))/((1-x)*(1-x))+1.74*Math.log(x)));
  double sonuc2=-R*Tc*1e-3*(1.45*x0+0.45*Math.log(x0)+0.25*Omega*(17.11*x0+
  25.2*(Math.pow(x0,(-2/3))/3.0*(1-x0)+Math.pow(x0,(1/3)))/((1-x0)*(1-x0))+1.74*Math.log(x0)));
  return (sonuc1-sonuc2);
}


public double dh_TV_LK_V(double T, double V)
{
//This method calculate difference between ideal gas enthalpy
//and real gas enthalpy
   double dH0,dHR;
   double Tr=T/Tc;
   double Vr=Pc*1e5*V/(R*Tc);
   double B0,C0,D0;
   double BR,CR,DR;
   double Z0,ZR,ZX;
   B0=b0[0]-b0[1]/Tr-b0[2]/(Tr*Tr)-b0[3]/(Tr*Tr*Tr);
   C0=c0[0]-c0[1]/Tr+c0[2]/(Tr*Tr*Tr);
   D0=d0[0]+d0[1]/Tr;
   Z0=1+B0/Vr+C0/(Vr*Vr)+D0/(Vr*Vr*Vr*Vr*Vr)
   + c0[3]/(Tr*Tr*Tr*Vr*Vr)*(beta0+gama0/(Vr*Vr))*Math.exp(-gama0/(Vr*Vr));
   BR=bR[0]-bR[1]/Tr-bR[2]/(Tr*Tr)-bR[3]/(Tr*Tr*Tr);
   CR=cR[0]-cR[1]/Tr+cR[2]/(Tr*Tr*Tr);
   DR=dR[0]+dR[1]/Tr;
   ZR=1+BR/Vr+CR/(Vr*Vr)+DR/(Vr*Vr*Vr*Vr*Vr)
   + cR[3]/(Tr*Tr*Tr*Vr*Vr)*(betaR+gamaR/(Vr*Vr))*Math.exp(-gamaR/(Vr*Vr));
   ZX=Z0+Omega/OmegaR*(ZR-Z0);
   double E0=c0[3]/(2.0*Tr*Tr*Tr*gama0)*
    (beta0+1-(beta0+1+gama0/Vr/Vr)*Math.exp(-gama0/Vr/Vr));
   double dh0=-Tr*(Z0-1-(b0[1]+2.0*b0[2]/Tr+3.0*b0[3]/(Tr*Tr))/(Tr*Vr)
    -(c0[1]-3.0*c0[2]/Tr/Tr)/(2.0*Tr*Vr*Vr)
    +d0[1]/(5*Tr*Vr*Vr*Vr*Vr*Vr)
    +3.0*E0);
    double ER=cR[3]/(2*Tr*Tr*Tr*gamaR)*
    (betaR+1-(betaR+1+gamaR/Vr/Vr)*Math.exp(-gamaR/Vr/Vr));
    double dhR=-Tr*(ZR-1-(bR[1]+2*bR[2]/Tr+3*bR[3]/(Tr*Tr))/(Tr*Vr)
    -(cR[1]-3*cR[2]/Tr/Tr)/(2*Tr*Vr*Vr)
    +dR[1]/(5*Tr*Vr*Vr*Vr*Vr*Vr)
    +3*ER);
    double dH1=1.0/OmegaR*(dhR-dh0);
    double dH=dh0+Omega/OmegaR*(dhR-dh0);
    dH*=R*Tc*1e-3;
    return dH;
}

public double ds_TV_LK_V(double T, double V)
{
// difference between ideal gas entropy and real gas entropy
   double dS0,dSR;
   double Tr=T/Tc;
   double Vr=Pc*1e5*V/(R*Tc);
   double B0,C0,D0,E0;
   double BR,CR,DR,ER;
   double Z0,ZR,ZX;
   B0=b0[0]-b0[1]/Tr-b0[2]/(Tr*Tr)-b0[3]/(Tr*Tr*Tr);
   C0=c0[0]-c0[1]/Tr+c0[2]/(Tr*Tr*Tr);
   D0=d0[0]+d0[1]/Tr;
   E0=c0[3]/(2.0*Tr*Tr*Tr*gama0)*
   (beta0+1.0-(beta0+1.0+gama0/Vr/Vr)*Math.exp(-gama0/Vr/Vr));
   Z0=1+B0/Vr+C0/(Vr*Vr)+D0/(Vr*Vr*Vr*Vr*Vr)
   + c0[3]/(Tr*Tr*Tr*Vr*Vr)*(beta0+gama0/(Vr*Vr))*Math.exp(-gama0/(Vr*Vr));
   BR=bR[0]-bR[1]/Tr-bR[2]/(Tr*Tr)-bR[3]/(Tr*Tr*Tr);
   CR=cR[0]-cR[1]/Tr+cR[2]/(Tr*Tr*Tr);
   DR=dR[0]+dR[1]/Tr;
   ER=cR[3]/(2.0*Tr*Tr*Tr*gamaR)*
   (betaR+1.0-(betaR+1.0+gamaR/Vr/Vr)*Math.exp(-gamaR/Vr/Vr));
   ZR=1+BR/Vr+CR/(Vr*Vr)+DR/(Vr*Vr*Vr*Vr*Vr)
   + cR[3]/(Tr*Tr*Tr*Vr*Vr)*(betaR+gamaR/(Vr*Vr))*Math.exp(-gamaR/(Vr*Vr));
   ZX=Z0+Omega/OmegaR*(ZR-Z0);
   double P=P_TV_LK_V(T,V);
   double ds0=-Math.log(1.01325/P)-Math.log(Z0)+(b0[0]+b0[2]/Tr/Tr+2.0*b0[3]/(Tr*Tr*Tr))/Vr
    +(c0[0]-2.0*c0[2]/Tr/Tr/Tr)/(2.0*Vr*Vr)
    +d0[0]/(5.0*Vr*Vr*Vr*Vr*Vr)-2.0*E0;
   double dsR=-Math.log(1.01325/P)-Math.log(ZR)+(bR[0]+bR[2]/Tr/Tr+2.0*bR[3]/(Tr*Tr*Tr))/Vr
    +(cR[0]-2.0*cR[2]/Tr/Tr/Tr)/(2.0*Vr*Vr)
    +dR[0]/(5.0*Vr*Vr*Vr*Vr*Vr)-2.0*ER;
   double dS=(ds0+Omega/OmegaR*(dsR-ds0))*R*1e-3;
   return dS;
}

double dhsh_TP_V(double T, double P)
{
  return h_TP_LK_V(T,P)-h_TP_LK_V(T,P_T_LK_SV(T));
}

double dh_TP_LK_V(double T,double P)
{
return dh_TV_LK_V(T,V_TP_LK_V(T,P));
}

double dh_TP_LK_L(double T,double P)
{
return dh_TV_LK_V(T,V_TP_LK_L(T,P));
}

double ds_TP_LK_V(double T,double P)
{
return ds_TV_LK_V(T,V_TP_LK_V(T,P));
}

//===================================================================
//termodynamic properties of real fluids

public double h_TP_BR_L(double T,double P)
{
// liquid enthalpy
return h_T_BR_SL(T)+V_T_BE_SL(T)*(P-P_T_LK_SV(P));
}

double h_TV_LK_V(double T, double V)
{
//enthalpy of the real fluid
// as a function of T and V
return h_T_IG_V(T)-h_T_IG_V(273.15)+dh_TV_LK_V(T,V)+50839.212;
}

double h_TP_LK_V(double T, double P)
{
//enthalpy of the real fluid
// as a function of T and P
return h_T_IG_V(T)-h_T_IG_V(273.15)+dh_TP_LK_V(T,P)+50839.212;
}

double h_T_LK_SV(double T)
{
 return h_TP_LK_V(T,P_T_LK_SV(T));
}

double h_T_LK_SL(double T)
{
 return h_TP_LK_V(T,P_T_LK_SV(T))-dH_T_LK_SVL(T);
}

double h_P_LK_SL(double P)
{
 double T=T_P_LK_SV(P);	
 return h_TP_LK_V(T,P_T_LK_SV(T))-dH_T_LK_SVL(T);
}

double h_P_LK_SV(double P)
{
 double T=T_P_LK_SV(P);	
 return h_TP_LK_V(T,P_T_LK_SV(T));
}

double h_TP_LK_L(double T, double P)
{
  // liquid enthalpy
return h_T_LK_SL(T)+V_T_LK_SL(T)*(P-P_T_LK_SV(P));
}



double s_TV_LK_V(double T, double V)
{
//entropy of the real fluid
// as a function of T and V
return s_TP_IG_V(T,P_TV_LK_V(T,V))+ds_TV_LK_V(T,V);
}

double s_TP_LK_V(double T, double P)
{
//entropy of the real fluid
// as a function of T and P
return s_TP_IG_V(T,P)+ds_TP_LK_V(T,P);
}

int zone(double T,double P)
{
/*
1 : supercritical fluid
2 : critical point
3 : superheated vapour
4 : saturated mixture
5 : liquid
6 : saturated vapour
7 : saturated liquid
*/
int z=0;
double Tvapour=T_P_SV(P);
if(T>=Tc) z=1;
else if(T==Tc) z=2;
else if(Tc>T && T> Tvapour) z=3;
else if (T==Tvapour && V_TP_LK_V(T,P)==V_T_LK_SV(Tvapour)) z=6;
else if (T==Tvapour && V_TP_LK_L(T,P)==V_T_LK_SL(Tvapour)) z=7;
else if(T==Tvapour) z=4;
else if(T < Tvapour) z=5;
return z;
}

String toString(double T,double P)
{
String s1="";
double Tvapour=T_P_SV(P);
if(T>=Tc) s1="supercritical fluid";
else if(T==Tc) s1="critical point";
else if(Tc>T && T> Tvapour) s1="superheated vapour";
else if (T==Tvapour && V_TP_LK_V(T,P)==V_T_LK_SV(Tvapour)) s1="saturated vapour";
else if (T==Tvapour && V_TP_LK_L(T,P)==V_T_LK_SL(Tvapour)) s1="saturated liquid";
else if(T==Tvapour) s1="saturated mixture";
else if(T < Tvapour) s1="liquid";
return s1;
}

String properties()
{  String s="isim   = "+GasEquation+"\n";
   s+="isim   = "+Name+"\n";
   s+="formul = "+Formula+"\n";
   s+="formul(atom sinifi) = "+toString()+"\n";
   s+="N = "+N+" M = "+M+"\n";
   s+="MolWt  = "+MolWt+"\n";
   s+="M      = "+M+"\n";
   s+="Tfp    = "+Tfp+"\n";
   s+="Tb     = "+Tb+"\n";
   s+="Tc     = "+Tc+"\n";
   s+="Pc     = "+Pc+"\n";
   s+="Vc     = "+Vc+"\n";
   s+="Zc     = "+Zc+"\n";
   s+="Omega  = "+Omega+"\n";
   s+="Dipm   = "+Dipm+"\n";
   s+="CPVAPA = "+CPC[0]+"\n";
   s+="CPVAPB = "+CPC[1]+"\n";
   s+="CPVAPC = "+CPC[2]+"\n";
   s+="CPVAPD = "+CPC[3]+"\n";
   s+="DELHF  = "+DELHF+"\n";
   s+="DELGF  = "+DELGF+"\n";
   s+="DELS0  = "+DELS0+"\n";
   s+="VPEq   = "+VPEq+"\n";
   s+="VPA    = "+VPC[0]+"\n";
   s+="VPB    = "+VPC[1]+"\n";
   s+="VPC    = "+VPC[2]+"\n";
   s+="VPD    = "+VPC[3]+"\n";
   s+="Tmin   = "+Tmin+"\n";
   s+="Tmax   = "+Tmax+"\n";
   s+="LDEN   = "+LDEN+"\n";
   s+="TDEN   = "+TDEN+"\n";
   return s;
}

}
