//======================================================
// Thermodynamic Package in java
// Class Steam Properties of water/steam
// Dr. Turhan Coban
// University of Nebrasca
// Mechanical Engineering, N104 WSEC, Lincoln NE 68588-0656
// phone :
// email : MCOBAN2@netscape.net
//         MCOBAN2@unl.edu
// File Name   : steam.java
// This file contains the Steam class
// this class sets thermophysical properties of steam
// =====================================================
// Description  : This file contains the CubicEquationofState class
//                This class calculates real fluid properties by
//                using one of the following cubic equations
//                van der Waals
//                Redlich-Kwong
//                Soave
//                Peng Robinson
//
// DATA FILE DEFINATION
// datas are written in the data file "RealGas.txt"
//============================================================
//  VARIABLE IDENTIFICATION
//============================================================
//  GasEquation : String type variable
//                values : "van der Waals"
//                         "Redlich Kwong"
//                         "Soave"
//                         "Peng Robinson"
//
//  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
//  CPVAPA,CPVAPB,CPVAPC,CPVAPD  : constants 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
//  Vpeq        : vapor pressure equation number
//  VPA,VPB,VPC,VPD : vapor pressure equation coefficients
//  Tmin        : minimum temp. 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, 4th edition
//             The Properties of Gases & Liquids,
//             Bruce E. Pouling, John M. Prausnitz, John P. O'Connell
////           McGraw-Hill Book Company, ISBN 0-07-011682-2, 5th edition
//============================================================
//
import java.io.*;

//==============================================
// implementation of constructors functions
// for class cubicEquationofState

public class cubicEquationofState
{
   String GasEquation;
   int natom;
   Atom atomList[];
   int No;
   String Formula,Name;
   double M,MolWt,Tfp,Tb,Tc,Pc,Vc,Zc,Omega,Dipm;
   double CPVAPA,CPVAPB,CPVAPC,CPVAPD,DELHF,DELGF;
   int VPEq;
   double VPA,VPB,VPC,VPD,Tmin,Tmax,LDEN,TDEN;
   int ierror;
   BufferedReader fin;
   double R=8314.5; // (J/kmolK)
   //================================//
   // Equation of State coefficients //
   //================================//
   double a,b,u,w;
   public cubicEquationofState(String  s1,double le,String ge) throws IOException
   {
   if(ge.equals("van der Waals") || ge.equals("Redlich Kwong") ||
   ge.equals("Soave") || ge.equals("Peng Robinson"))
   {GasEquation=ge;}
   else
   {
   GasEquation="Peng Robinson";
   }
   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 aCPVAPA,aCPVAPB,aCPVAPC,aCPVAPD,aDELHF,aDELGF;
   int aVPEq;
   double aVPA,aVPB,aVPC,aVPD,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);
         aCPVAPA=Text.readDouble(fin);
         aCPVAPB=Text.readDouble(fin);
         aCPVAPC=Text.readDouble(fin);
         aCPVAPD=Text.readDouble(fin);
         aDELHF=Text.readDouble(fin);
         aDELGF=Text.readDouble(fin);
         aVPEq=Text.readInt(fin);
         aVPA=Text.readDouble(fin);
         aVPB=Text.readDouble(fin);
         aVPC=Text.readDouble(fin);
         aVPD=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;
         CPVAPA=aCPVAPA;
         CPVAPB=aCPVAPB;
         CPVAPC=aCPVAPC;
         CPVAPD=aCPVAPD;
         DELHF=aDELHF;
         DELGF=aDELGF;
         VPEq=aVPEq;
         VPA=aVPA;
         VPB=aVPB;
         VPC=aVPC;
         VPD=aVPD;
         Tmin=aTmin;
         Tmax=aTmax;
         LDEN=aLDEN;
         TDEN=aTDEN;
         natom=anatom;
         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");
        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 void setconstants(double T)
   {
     if(GasEquation.equals("van der Waals"))
     {
     u=0;
     w=0;
     b=8314.5*Tc/(8.0*Pc*1e5);
     a=27.0/64.0*(R*R*Tc*Tc)/(Pc*1e5);
     }
     else if(GasEquation.equals("Redlich Kwong"))
     {
     u=1;
     w=0;
     b=0.08664*R*Tc/(Pc*1e5);
     a=0.42748*R*R*Math.pow(Tc,2.5)/(Pc*1e5*Math.sqrt(T));
     }
     else if(GasEquation.equals("Soave"))
     {
     u=1;
     w=0;
     double fw=0.48+1.57*Omega-0.176*Omega*Omega;
     double Tr=T/Tc;
     b=0.08664*R*Tc/(Pc*1e5);
     double xy=1.0+fw*(1.0-Math.sqrt(Tr));
     a=0.42748*R*R*Tc*Tc/(Pc*1e5)*xy*xy;
     }
     else if(GasEquation.equals("Peng Robinson"))
     {
     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);
     double xy=1.0+fw*(1.0-Math.sqrt(Tr));
     a=0.45724*R*R*Tc*Tc/(Pc*1e5)*xy*xy;
     }
    }

   double P(double T, double V)
   {
   setconstants(T);
   double v=V;
   return (R*T/(v-b)-a/(v*v+u*b*v+w*b*b))*1e-5;
   }

   public double Z(double T, double P)
   {
   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=0;
   for(int i=0;i<3;i++)
    {
    if((zc[i].nreal>=0.0 && zc[i].nreal<=1.0) &&
       (zc[i].nimag>=0.0 && zc[i].nimag<1e-4))
       {
       zz=zc[i].nreal;
       break;
       }
    }
   return zz;
   }

   public double V(double T, double P)
   {
   return Z(T,P)*R*T/(P*1e5);
   }

   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;
   }

}
