//======================================================
// Thermodynamic Package in java
// Class Steam Properties of water/steam
// Dr. Turhan Coban
// TUBITAK Ulusal Metroloji Enstitüsü
// email : Turhan.Coban@ume.tubitak.gov.tr
//
// File Name   : steam.java
// This file contains the Steam class
// this class sets thermophysical properties of steam
// =====================================================
// Description : This file contains the  steam class
// class steam calculates thermophysical properties of
// water and steam. It utilise Keenan & Keyes equation of state
// DATA FILE DEFINATION
// required datas are written directly inside of the class as array structures
// reference : Steam tables, thermodynamic properties of water including vapor,
// liquid and solid phases (SI unites), Joseph H. Keenan, Frederick G. Keyes,
// Philip G. Hill, Joan G. Moore, Willey-Interscience publication,
// John Willey and Sons, 1969, ISBN 0-471-04210-2
// VERSION 1.0 : 25/9/2001
//============================================================
//  EQUATION OF STATE
//============================================================
// fi=fi0(T)+R*T[ln(ro)+ro*Q(ro,to)]
// helmholts free energy function
// fi0(T)=sum(i=0;i<=6){C[i]/to^i+C[6]*ln(T)+C[7]*lnT/to}
// helmholts free energy function for the reference pressure
// Q=(to-toc)*sum(j=0;j<=6){sum(i=0;i<=7){A[i][j]*(ro-roaj)^i}
// +{exp(-E*ro)*sum(i=8;i<=9){A[i][j]*ro^(i-8)}}
//============================================================
//  VARIABLE IDENTIFICATION
//============================================================
//  A : array coefficient of eqaution of state
//  u : switch for unit system, SI or ENglish : default state : SI
// unit : the assigned value of the unit
//  C : vector coefficient of the equation of state
//  R1: Gas constant bar/cm^3/g K
//  R2: Gas constant J/g K
//  E : constant of the equation of state
//  toc :1000/Tc
//  to  : 1000/T
//  F : coefficients of the saturation state
//  Pc: critical pressure bar
//  Tc: critical temperature K
// toaj: constant of equation of state toaj=toc(for j=0),toaj=2.5(for j>0)
// roaj: constant of equation of state roaj=0.634(for j=0),roaj=1.0(for j>0)
//============================================================
//  METHOD IDENTIFICATION
//============================================================
// deriv_dQdro=dQ/dro : derivative of Q with respect to ro
// setUnit : method to change unit system
// B : second virial coefficient
// Q : sub function of equation of state
// dQdto = dQ/dto : derivative of Q with respect to to
// dQdro = dQ/dro : derivative of Q with respect to ro
// d2Qdro2 = d2Q/dro2 : second derivative of Q with respect to ro
// d2Qdro2 = d2Q/drodto : derivative of dQ/dro with respect to to
// fi0 : helmholts free energy function for the reference pressure
// P : pressure as a function of temperature and specific volume
// Pi : pressure as a function of temperature and specific volume
//      also predict saturation presure if the volume in saturation region
// dPdv : derivative of P with respect to v
// dPdt : derivative of P with respect to t
// property : method to calculate several thermodynamic properties
// u    : internal energy as a function of t and v
// h    : enthalpy as a function of t and v
// s    : entropy as a function of t and v
// Ps   : saturated vapor pressure as a function of T
// vsl_t :saturated liquid specific volume as a function of T
// vsg_t :saturated vapor specific volume as a function of T
// v_tp  : specific volume as a function of T and v
// v_th  : specific volume as a function of T and h
// v_tu  : specific volume as a function of T and u
// v_tx  : specific volume as a function of T and x
// v_ts  : specific volume as a function of T and s
// t_pv  : temperature as a function of pressure and specific volume
// t_ph  : temperature as a function of pressure and enthalpy
// t_pu  : temperature as a function of pressure and internal energy
// t_ps  : temperature as a function of pressure and internal entropy
// phase_EN : phase in Eglish unit
// x_tv : quality as a function of temperature and specific volume
// property_EN : method to calculate several thermodynamic properties in english units
// property_SI : method to calculate several thermodynamic properties in SI units
// toString : result output through a string variable
//v_tpPR : specific volume as a function of temperature and volume calculated by
//          peng robinson equation of state

public class steamT
{
final double A[][]={
{29.492937,-5.1985860,6.8335354,-0.1564104,-6.3972405,-3.9661401,-0.69048554},
{-132.13917,7.7779182,-26.149751,-0.72546108,26.409282,15.453061,2.7407416},
{274.64632,-33.301902,65.326396,-9.2734289,-47.740374,-29.14247,-5.1028070},
{-360.93828,-16.254622,-26.181978,4.312584,56.32313,29.568796,3.9636085},
{342.18431,-177.31074,0,0,0,0,0},
{-244.50042,127.48742,0,0,0,0,0},
{155.18535,137.46153,0,0,0,0,0},
{5.9728487,155.97836,0,0,0,0,0},
{-410.30848,337.3118,-137.46618,6.7874983,136.87317,79.847970,13.041253},
{-416.05860,-209.88866,-733.96848,10.401717,645.8188,399.1757,71.531353}
                  };
final String u[]={"SI","EN"};
String unit=u[0];
double toaj;
double roaj;
final double C[]={1857.065,3229.12,-419.465,36.6649,
-20.5516,4.85233,46.0,-1011.249}; //  J/gr
final double R1=4.6150943;  //bar cm^3/g K
final double R2=0.46150943; // J/g K
final double E=4.8;
final double toc=1.544912;
final double F[]={-741.9242,-29.721,-11.55286,-0.8685635,0.1094098,0.439993,
0.2520658,0.05218684};
final double Pc=220.88;   //bar
final double Tc=647.286;  //K
final double vc=3.155e-3; //m^3/kg

public double deriv_dQdro(double ro, double to)
{
// This method calculates derivatives of a simple function
// accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_x
double x=ro;
double h0=0.0256;
int i,m;
int n=7;
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
T[i][0]=( Q((x + h[i]),to) - Q((( x - h[i])/(2.0*h[i])),to));
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] - h[i+m]*h[i+m]*T[i][m-1])/
          (h[i]*h[i]- h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}

void setUnit(boolean x)
{
//sets SI or English units
if(x) unit=u[0];else unit=u[1];
}

void setUnit(String x)
{
//sets SI or English units
if(x.equals("EN")) unit=u[1];
else unit=u[0];
}

double B(double t)
{
//second virial coefficient
double T=t+273.15; // K
double to=1000/T;
//second virial coefficient
return 2.0624-2.61204*to*Math.pow(10,(0.1008*to*to/(1+0.0349*to*to)));
}

double Q(double ro, double to)
{
//coefficient of steam state equation, fi
  double q=0;
  for(int j=0;j<7;j++)
  {
  double s1=0;
  double s2=0;
  if(j==0) {toaj=toc;roaj=0.634;}
  else     {toaj=2.5;roaj=1.0;  }
    s1=A[0][j];
    for(int i=1;i<8;i++)
    {
    s1+=A[i][j]*Math.pow((ro-roaj),i);
    }
    s2=A[8][j]+A[9][j]*ro;
  q+=Math.pow((to-toaj),(j-1))*(s1+Math.exp(-E*ro)*s2);
  }
  q*=(to-toc);
  return q;
}

double dQdto(double ro, double to)
{
 // derivative of Q respect to to (temperature function)
  double q1=0;
  double q2=0;
  double toij;
  for(int j=0;j<7;j++)
  {
  double s1=0;
  double s2=0;
  if(j==0) {toaj=toc;roaj=0.634;}
  else     {toaj=2.5;roaj=1.0;  }
    s1=A[0][j];
    for(int i=1;i<8;i++)
    {
    s1+=A[i][j]*Math.pow((ro-roaj),i);
    }
    s2=A[8][j]+A[9][j]*ro;
    q1+=Math.pow((to-toaj),(j-1))*(s1+Math.exp(-E*ro)*s2);
    q2+=(j-1)*Math.pow((to-toaj),(j-2))*(s1+Math.exp(-E*ro)*s2);
  }
  return q1+q2*(to-toc);
}

double dQdro(double ro, double to)
{
 // derivative of Q respect to density (specific volume  function)
  double q=0;
  double toij;
  for(int j=0;j<7;j++)
  {
  double s1=0;
  double s2=0;
  double s3=0;
  if(j==0) {toaj=toc;roaj=0.634;}
  else     {toaj=2.5;roaj=1.0;  }
    for(int i=1;i<8;i++)
    {
    s1+=i*A[i][j]*Math.pow((ro-roaj),(i-1));
    }
    s2+=A[8][j]+A[9][j]*ro;
    s3+=A[9][j];
  q+=Math.pow((to-toaj),(j-1))*(s1-E*Math.exp(-E*ro)*s2+Math.exp(-E*ro)*s3);
  }
  q*=(to-toc);
  return q;
}

double dQ2drodt(double ro, double to)
{
 // derivative of Q respect to density (specific volume  function)
  double q=0;
  double toij;
  for(int j=0;j<7;j++)
  {
  double s1=0;
  double s2=0;
  double s3=0;
  if(j==0) {toaj=toc;roaj=0.634;}
  else     {toaj=2.5;roaj=1.0;  }
    for(int i=1;i<8;i++)
    {
    s1+=i*A[i][j]*Math.pow((ro-roaj),(i-1));
    }
    s2+=A[8][j]+A[9][j]*ro;
    s3+=A[9][j];
  q+=(Math.pow((to-toaj),j-1)+(to-toc)*(j-1)*Math.pow((to-toaj),j-2))*
  (s1-E*Math.exp(-E*ro)*s2+Math.exp(-E*ro)*s3);
  }
  return q;
}


double dQ2dro2(double ro, double to)
{
// second derivative of Q respect to density (specific volume  function)
  double q=0;
  double toij;
  for(int j=0;j<7;j++)
  {
  double s1=0;
  double s2=0;
  double s3=0;
  if(j==0) {toaj=toc;roaj=0.634;}
  else     {toaj=2.5;roaj=1.0;  }
    for(int i=0;i<8;i++)
    {
    s1+=i*(i-1)*A[i][j]*Math.pow((ro-roaj),(i-2));
    }
    for(int i=8;i<10;i++)
    {
    s2+=A[i][j]*Math.pow(ro,(i-8));
    }
    for(int i=8;i<10;i++)
    {
    s3+=(i-8)*(i-9)*A[i][j]*Math.pow(ro,(i-10));
    }
  q+=Math.pow((to-toaj),(j-1))*(s1+E*E*Math.exp(-E*ro)*s2+Math.exp(-E*ro)*s3);
  }
  q*=(to-toc);
  return q;
}


double fi0(double to)
{
// helmholts free energy function for reference pressure

 double ff0=0;
 for(int i=0;i<6;i++)
 {
  ff0+=C[i]/Math.pow(to,i);
 }
 ff0+=C[6]*Math.log(1000.0/to)+C[7]*Math.log(1000.0/to)/to;
 return ff0;
}

double dfi0todto(double to)
{
// derivative helmholts free energy function for reference pressure
  double s=0;
  s+=C[0]-C[2]/(to*to)-2*C[3]/(to*to*to)-3*C[4]/(to*to*to*to)-4*C[5]/(to*to*to*to*to)+
  C[6]*Math.log(1000.0/to)-C[6]-C[7]/to;
  return s;
}

double dfidt(double to)
{
  double s=0;
  s+=-C[1]/(to*to)-2*C[2]/(to*to*to)-3*C[3]/(to*to*to*to)
  -4*C[4]/(to*to*to*to*to)-5*C[5]/(to*to*to*to*to*to)+
  -C[6]/to-C[7]/(to*to)*(Math.log(1000.0/to)+1.0);
  s*=(-to*to/1000);
  return s;
}

double fi(double t, double v)
{
  //helmholts free energy
  //actual equation of state  for steam
  // T C
  // v m^3/kg
  double T=t+273.15; //K
  double to=1000/T;
  double ro=1/v*1e-3;
  return fi0(to)+R2*1000.0/to*(Math.log(ro)+ro*Q(ro,to));
}

double P(double t, double v)
{
// pressure p(t,v)
// T C
// v m^3/kg
double T=t+273.15; //K
double to=1000/T;
double ro=1/v*1e-3;
return ro*R1*T*(1+ro*Q(ro,to)+ro*ro*dQdro(ro,to));
}

double Pi(double t, double v)
{
// pressure p(t,v)
// T C
// v m^3/kg
double T=t+273.15; //K
double to=1000/T;
double ro=1/v*1e-3;
double vsl=vsl_t(t);
double vsg=vsg_t(t);
if(vsl<=v && v>=vsg) {return Ps(t);}
return P(t,v);
}


double dPdv(double t, double v)
{
//derivative of pressure dp(t,v)/dv
double T=t+273.15; //K
double to=1000/T;
double ro=1/v*1e-3;
return R1*T*(1.0+2.0*ro*Q(ro,to)+4.0*ro*ro*dQdro(ro,to)+
ro*ro*ro*dQ2dro2(ro,to))*(-1/(v*v)*1e-3);
}

double dPdt(double t, double v)
{
double T=t+273.15; //K
double to=1000/T;
double ro=1/v*1e-3;
   return ro*R1*(1+ro*Q(ro,to)+ro*ro*dQdro(ro,to))+
        ro*R1*(ro*dQdto(ro,to)+ro*ro*dQ2drodt(ro,to))*(-1000/(T));
}

double u(double t, double v)
{
//internal energy
double T=t+273.15; //K
double to=1000/T;
double ro=1/v*1e-3;
return R2*T*ro*to*dQdto(ro,to)+dfi0todto(to);
}

public double h(double t, double v)
{
//enthalpy
double T=t+273.15; //K
double to=1000/T;
double ro=1/v*1e-3;
return R2*T*(ro*to*dQdto(ro,to)+1+ro*Q(ro,to)+ro*ro*dQdro(ro,to))+dfi0todto(to);
}

public double s(double t, double v)
{
//entropy
double T=t+273.15; //K
double to=1000/T;
double ro=1/v*1e-3;
return -R2*(Math.log(ro)+ro*Q(ro,to)-ro*to*dQdto(ro,to))-dfidt(to);
}


// saturation state
//=========================================================

public double Ps(double Ts)
{
  // pressure of saturated steam
  double to=1000.0/(Ts+273.15);
  double ttc=374.136;
  double s=0;
  for(int i=0;i<8;i++)
  {
    s+=F[i]*Math.pow((0.65-0.01*Ts),i);
  }
  return Pc*Math.exp(to*1e-5*(ttc-Ts)*s);}

public double dPs(double Ts)
{
  double to=1000.0/(Ts+273.15);
  double ttc=374.136;
  double s1=0;
  double s2=0;
  for(int i=0;i<8;i++)
  {
    s1+=F[i]*Math.pow((0.65-0.01*Ts),i);
    s2+=F[i]*i*Math.pow((0.65-0.01*Ts),(i-1));
  }
  return (-to*to*1e-8*(ttc-Ts)*s1-
  to*1e-5*s1+to*1e-5*(ttc-Ts)*s2)*Pc*Math.exp(to*1e-5*(ttc-Ts)*s1);
}

public double Ts(double ps)
{
double f0,df0;
int nmax=150;
//you can adjust tolerance if you require better accuracy
double ttc=374.136;
int i=0;
double aa[]=new double[6];
aa[0] = 99.61598078344355;
aa[1] = 27.87371293352203;
aa[2] = 2.4012458375125334;
aa[3] = 0.2298589975523555;
aa[4] = 0.018713728673908747;
aa[5] = 5.666201496623557E-4;
double pl=Math.log(ps);
double pl2=pl*pl;
double pl4=pl2*pl2;
double x=aa[0]+aa[1]*pl+aa[2]*pl2+aa[3]*pl2*pl+aa[4]*pl4+aa[5]*pl4*pl;
double dx=1.0e-5;
double dx1=2.0e-5;
double x1m;
double tolerance=1.0e-10;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
for(i=0;i<nmax;i++)
{
fx=Ps(x)-ps;
fxm=Ps((x-dx))-ps;
fxm1=Ps(x-dx1)-ps;
dfx=(fx-fxm)/dx;
dfxm=(fx-fxm1)/dx1;
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}

public double vsl_t(double t)
{
// specific volume of saturated liquid
// secant method
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
// function should be supplied
if(t>=Tc) return 3.155e-3;
double x=1.001;
int nmax=500;
double dx=1.0e-7;
double dx1=2.0e-7;
double x1m;
double tolerance=1.0e-10;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
if((t+273.15) > Tc) return 3.155e-3;
for(int i=0;i<nmax;i++)
{
fx=P(t,x*1e-3)-Ps(t);
fxm=P(t,(x-dx)*1e-3)-Ps(t);
fxm1=P(t,(x-dx1)*1e-3)-Ps(t);
dfx=(fx-fxm)/dx;
dfxm=(fx-fxm1)/dx1;
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x*1e-3;
}
}
return x*1e-3;
}


public double vsg_t(double t)
{
// specific volume of the saturated vapor
// secant method
// a polynomial curve fitting is used as first guess
if( (t+273.15) >= Tc) return 3.155e-3;
double aa[]=new double[7];
aa[0] = 5.30262888569996;
aa[1] = -0.0637415944379427;
aa[2] = 1.118825833543724E-4;
aa[3] = 1.1808820813419274E-6;
aa[4] = -9.85398318281358E-9;
aa[5] = 2.8443958795501176E-11;
aa[6] = -2.9157288848571675E-14;
double x=Math.exp(aa[0]+aa[1]*t+aa[2]*t*t+aa[3]*t*t*t+aa[4]*t*t*t*t+
aa[5]*t*t*t*t*t+aa[6]*t*t*t*t*t*t);
int nmax=500;
double dx=1.0e-5;
double dx1=2.0e-5;
double x1m;
double tolerance=1.0e-8;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=P(t,x)-Ps(t);
fxm=P(t,(x-dx))-Ps(t);
fxm1=P(t,(x-dx1))-Ps(t);
dfx=(fx-fxm)/dx;
dfxm=(fx-fxm1)/dx1;
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}

public double v_tp(double t, double p)
{
double ps=Ps(t);
double ts=Ts(p);
double x=1.0002;
int nmax=100;
double dx=1.0e-4;
double dx1=2.0e-4;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
//using Peng robinson equation of state as first estimate
  if( p > Pc)
  {
  x=v_tpPR("v",t,p)*1e3;
     if( x*1e-3< vc )
     { double a1[]=new double[6];
       a1[0] = 0.9678434673818401;
       a1[1] = 0.003343621695478305;
       a1[2] = -7.512696317379941E-5;
       a1[3] = 6.649577638309645E-7;
       a1[4] = -2.315472955248364E-9;
       a1[5] = 2.8425996940750715E-12;
       double x1=a1[0]+a1[1]*t+a1[2]*t*t+a1[3]*t*t*t+a1[4]*t*t*t*t+a1[5]*t*t*t*t*t;
       double a[]=new double[4];
       a[0] = 1.0435466274372924;
       a[1] = -5.11086759702184E-5;
       a[2] = 9.36450821986505E-9;
       a[3] = -1.7399029542296548E-12;
       double x2=a[0]+a[1]*p+a[2]*p*p+a[3]*p*p*p;
       x=x2*x1/1.0313;
     }
  }
  else if(p <= Pc)
  {
     double vsl=vsl_t(ts);
     double vsg=vsl_t(ts);
     double a[]=new double[4];
     a[0] = 1.0435466274372924;
     a[1] = -5.11086759702184E-5;
     a[2] = 9.36450821986505E-9;
     a[3] = -1.7399029542296548E-12;
     double xx=a[0]+a[1]*p+a[2]*p*p+a[3]*p*p*p;
     if(t < ts) x=xx/vsl_t(100)*vsl_t(t);
     else if (t>ts)  x=v_tpPR("v",t,p)*1e3;//x=(8314.5/18.02*(t+273.15)/(p*1.0e5)*1e3+vsg_t(t))/2.0;
     if(t==ts) x=v_tpPR("v",t,p)*1e3;
  }
  for(int i=0;i<nmax;i++)
  {
  fx=P(t,x*1e-3)-p;
  fxm=P(t,(x-dx)*1e-3)-p;
  fxm1=P(t,(x-dx1)*1e-3)-p;
  dfx=(fx-fxm)/dx;
  dfxm=(fx-fxm1)/dx1;
  d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
  x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
  if(Math.abs(fx)<tolerance)
    { return x*1e-3; }
  }
  return x*1e-3;
}

public double v_th(double t,double hi)
{
// bisection method
// defination of variables :
// xl : lower guess
// xu : upper guess
// xr : root estimate
// es : stopping criterion
// ea :approximate error
// maxit : maximum iterations
// iter  : number of iteration
double xl=0.94;
double xu=726025;
double test;
double xr=0;
double es,ea;
double fxl,fxr;
double tolerance=1e-20;
int maxit=500,iter=0;
es=0.000001;
ea=1.1*es;
     double ps=Ps(t);
     double vsl=vsl_t(t);
     double vsg=vsg_t(t);
     double hsl=h(t,vsl);
     double hsg=h(t,vsg);
    if( (t+273.15) > Tc)
    { if(hi< hsl) {xl=0.94;xu=vsl*1e3;}
      else  {xl=vsl*1e3;xu=726025;}
    }
    else if((t+273.15) < Tc)
    {
      double xx1=(hi-hsl)/(hsg-hsl);
      double xx=(hsl-hi)/hsl;
      if(hi<=hsl*1.1) {xl=vsl*1e3*0.9;xu=1.1*vsl*1e3;}
      //else if (hi>=hsl && hi<=hsg) {return vsl*(1-xx1)+vsg*xx1;}
      else if (hi>=hsg) {xl=vsg_t(t)*1e3;xu=726025;}
      else {xl=vsl*1e3*0.9;xu=xu=1.1*vsg*1e3;}
    }
while((ea>es)&&(iter<maxit))
{
xr=(xl+xu)/2.0;
iter++;
if((xl+xu)!=0)
  { ea=Math.abs((xu-xl)/(xl+xu))*100;}
fxl=h(t,xl*1e-3)-hi;
fxr=h(t,xr*1e-3)-hi;
test= fxl*fxr;
if(fxl==0) return xl;
if(test==0.0 || (Math.abs(fxl) < tolerance) )       ea=0;
else if(test<0.0)  xu=xr;
else
  {
  xl=xr;
  }
} //end of while
return xr*1e-3;
}

public double v_tu(double t,double ui)
{
//bisection method
// defination of variables :
// xl : lower guess
// xu : upper guess
// xr : root estimate
// es : stopping criterion
// ea :approximate error
// maxit : maximum iterations
// iter  : number of iteration
double xl=0.94;
double xu=726025;
double test;
double xr=0;
double es,ea;
double fxl,fxr;
double tolerance=1e-20;
int maxit=500,iter=0;
es=0.000001;
ea=1.1*es;
     double ps=Ps(t);
     double vsl=vsl_t(t);
     double vsg=vsg_t(t);
     double usl=u(t,vsl);
     double usg=u(t,vsg);
    if( (t+273.15) > Tc)
    { if(ui< usl) {xl=0.94;xu=vsl*1e3;}
      else  {xl=vsl*1e3;xu=726025;}
    }
    else if((t+273.15) < Tc)
    {
      double xx1=(ui-usl)/(usg-usl);
      double xx=(usl-ui)/usl;
      if(ui<=usl*1.1) {xl=vsl*1e3*0.9;xu=1.1*vsl*1e3;}
      //else if (ui>=usl && ui<=usg) {return vsl*(1-xx1)+vsg*xx1;}
      else if (ui>=usg) {xl=vsg_t(t)*1e3;xu=726025;}
      else {xl=vsl*1e3*0.9;xu=xu=1.1*vsg*1e3;}
    }
while((ea>es)&&(iter<maxit))
{
xr=(xl+xu)/2.0;
iter++;
if((xl+xu)!=0)
  { ea=Math.abs((xu-xl)/(xl+xu))*100;}
fxl=u(t,xl*1e-3)-ui;
fxr=u(t,xr*1e-3)-ui;
test= fxl*fxr;
if(fxl==0) return xl;
if(test==0.0 || (Math.abs(fxl) < tolerance) )       ea=0;
else if(test<0.0)  xu=xr;
else
  {
  xl=xr;
  }
} //end of while
return xr*1e-3;
}


public double v_ts(double t,double si)
{
// bisection method
// defination of variables :
// xl : lower guess
// xu : upper guess
// xr : root estimate
// es : stopping criterion
// ea :approximate error
// maxit : maximum iterations
// iter  : number of iteration
double xl=0.94;
double xu=726025;
double test;
double xr=0;
double es,ea;
double fxl,fxr;
double tolerance=1e-20;
int maxit=500,iter=0;
es=0.000001;
ea=1.1*es;
     double ps=Ps(t);
     double vsl=vsl_t(t);
     double vsg=vsg_t(t);
     double ssl=s(t,vsl);
     double ssg=s(t,vsg);
  if( (t+273.15) > Tc)
  { if(si< ssl) {xl=0.94;xu=vsl*1e3;}
    else  {xl=vsg*1e3;xu=726025;}
  }
  else if((t+273.15) < Tc)
  {
     double xx=(ssl-si)/ssl;
     double xi=(si-ssl)/(ssg-ssl);
     if(si==ssl) return vsl;
     else if(si==ssg) return vsg;
     //else if((si>ssl) && (si<ssg) ) {return (1-xi)*vsl+xi*vsg;}
     else if(si<ssl) {xl=0.94;xu=1.1*vsl*1e3;}
     else if (si>ssg) {xl=vsg_t(t)*1e3;xu=726025;}
  }
while((ea>es)&&(iter<maxit))
{
xr=(xl+xu)/2.0;
iter++;
if((xl+xu)!=0)
  { ea=Math.abs((xu-xl)/(xl+xu))*100;}
fxl=s(t,xl*1e-3)-si;
fxr=s(t,xr*1e-3)-si;
test= fxl*fxr;
if(fxl==0) return xl;
if(test==0.0 || (Math.abs(fxl) < tolerance) )       ea=0;
else if(test<0.0)  xu=xr;
else
  {
  xl=xr;
  }
} //end of while
return xr*1e-3;
}

public double t_ph(double pi,double hi)
{
double xl=0.1;
double xu=3000;
if(pi<Pc)
{
double ts=Ts(pi);
double vsl=vsl_t(ts);
double vsg=vsg_t(ts);
double hsl=h(ts,vsl);
double hsg=h(ts,vsg);
if(hi<hsl) {xl=0.1;xu=ts;}
else if(hi>hsg) {xl=ts;xu=3000;}
else if(hi>=hsl && hi<=hsg) { return ts; }
}
int maxit=50,iter=0;
double es=0.000000001;
double dx=Math.abs((xu-xl)/maxit);
double dxp=dx;
double xr=xl;
double fx=hi-h(xr,v_tp(xr,pi));
double fxr=fx;
for(int i=0;i<maxit;i++)
{
xr+=dxp;
if(xr>xu) break;
fx=hi-h(xr,v_tp(xr,pi));
  if(fx*fxr<0)
  {
    if(dxp<es) return xr;
    else
    {
    xr-=dxp;
    dxp/=2.0;
    }
  }
}
return xr;
}

public double v_ph(double pi,double hi)
{
return v_tp(t_ph(pi,hi),pi);
}

public double t_pu(double pi,double ui)
{
double xl=0.1;
double xu=3000;
if(pi<Pc)
{
double ts=Ts(pi);
double vsl=vsl_t(ts);
double vsg=vsg_t(ts);
double usl=h(ts,vsl);
double usg=h(ts,vsg);
if(ui<usl) {xl=0.1;xu=ts;}
else if(ui>usg) {xl=ts;xu=3000;}
else if(ui>=usl && ui<=usg) {return ts;}
}
int maxit=50,iter=0;
double es=0.000000001;
double dx=Math.abs((xu-xl)/maxit);
double dxp=dx;
double xr=xl;
double fx=ui-u(xr,v_tp(xr,pi));
double fxr=fx;
for(int i=0;i<maxit;i++)
{
xr+=dxp;
if(xr>xu) break;
fx=ui-u(xr,v_tp(xr,pi));
  if(fx*fxr<0)
  {
    if(dxp<es) return xr;
    else
    {
    xr-=dxp;
    dxp/=2.0;
    }
  }
}
return xr;
}

public double v_pu(double pi,double ui)
{
return v_tp(t_pu(pi,ui),pi);
}

public double t_ps(double pi,double si)
{
double xl=0.1;
double xu=3000;
if(pi<Pc)
{
double ts=Ts(pi);
double vsl=vsl_t(ts);
double vsg=vsg_t(ts);
double ssl=s(ts,vsl);
double ssg=s(ts,vsg);
if(si<ssl) {xl=0.1;xu=ts;}
else if(si>ssg) {xl=ts;xu=3000;}
else if(si>=ssl && si<=ssg) {return ts;}
}
int maxit=50,iter=0;
double es=0.000000001;
double dx=Math.abs((xu-xl)/maxit);
double dxp=dx;
double xr=xl;
double fx=si-s(xr,v_tp(xr,pi));
double fxr=fx;
for(int i=0;i<maxit;i++)
{
xr+=dxp;
if(xr>xu) break;
fx=si-s(xr,v_tp(xr,pi));
  if(fx*fxr<0)
  {
    if(dxp<es) return xr;
    else
    {
    xr-=dxp;
    dxp/=2.0;
    }
  }
}
return xr;
}

public double v_ps(double pi,double si)
{
return v_tp(t_ps(pi,si),pi);
}

public double v_tx(double t, double x)
{
  double vsl=vsl_t(t);
  double vsg=vsg_t(t);
  return vsl*(1-x)+vsg*x;
}

public String phase(double t, double v)
{
//phase of the steam-water system
if( Ps(t) > Pc)
  return "Superkritik sıvı";
else
 {
 double ts=t;
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 if(v<vsl)
   return "Su";
 else if(v==vsl)
   return "Doymuş Su";
 else if(v==vsg)
   return "Doymuş Su-Buhar Karışımı";
 else if(v>vsl && v<vsg)
   return "saturated liquid-vapor mixture";
 else // (v>vsg)
   return "Kızgın Buhar";
 }
}

public String phase(double x)
{
//phase of the steam-water system
if( x==3)
  return "Superkritik sıvı";
else
 {
 if(x<0)
   return "Su";
 else if( x==0)
   return "Doymuş Su";
 else if(x==1)
   return "Doymuş Buhar";
 else if( x>0 && x<1)
   return "Doymuş Su-Buhar Karışımı";
 else // (x>1)
   return "Kızgın Buhar";
 }
}

public String phase_EN(double t, double v)
{
t=(t-32)/1.8;
v/=16.01877892;
return phase(t,v);
}

public double t_pv(double p, double v)
{
double ro=1e-3/v;
double ts=Ts(p);
double vsl=vsl_t(ts);
double vsg=vsg_t(ts);
double x=200.0;
int nmax=100;
double dx=1.0e-3;
double dx1=2.0e-3;
double x1m;
double tolerance=1.0e-8;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
  if( p > Pc)
  { x=p*v/(8314.5/18.02); }
  else if(p < Pc)
  {  x=ts*v/vsg; }
  for(int i=0;i<nmax;i++)
  {
  fx=P(x,v)-p;
  fxm=P((x-dx),v)-p;
  fxm1=P((x-dx1),v)-p;
  dfx=(fx-fxm)/dx;
  dfxm=(fx-fxm1)/dx1;
  d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
  x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
  if(Math.abs(fx)<tolerance)
    { return x; }
  }
  return x;
}

public double x_tv(double t, double v)
{
//phase and/or quality of the steam-water system
double p=Ps(t);
if(p> Pc)
  return 3;
else
 {
 double ts=t;
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 if(v<vsl)
   return (v-vsl)/(vsg-vsl);
 else if(v==vsl)
   return 0;
 else if(v==vsg)
   return 1;
 else if(v>vsl && v<vsg)
   return (v-vsl)/(vsg-vsl);
 else // (v>vsg)
   return (v-vsl)/(vsg-vsl);
 }
}

public double x_pv(double p, double v)
{
//phase and/or quality of the steam-water system
double ts=Ts(p);
if(p > Pc)
  return 3;
else
 {
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 if(v<vsl)
   return (v-vsl)/(vsg-vsl);
 else if(v==vsl)
   return 0;
 else if(v==vsg)
   return 1;
 else if(v > vsl && v < vsg)
   return (v-vsl)/(vsg-vsl);
 else // (v>vsg)
   return (v-vsl)/(vsg-vsl);
 }
}

public double[] property_EN(String s, double v1, double v2)
{
if(s.charAt(0)=='t') v1=(v1-32.0)/1.8;
else if(s.charAt(0)=='p') v1/=14.504;
else if(s.charAt(0)=='v') v1/=16.01877892;
if(s.charAt(1)=='t') v2=(v2-32)/1.8;
else if(s.charAt(1)=='p') v2/=14.504;
else if(s.charAt(1)=='v') v2/=16.01877892;
else if(s.charAt(1)=='h') v2/=0.42992;
else if(s.charAt(1)=='u') v2/=0.42992;
else if(s.charAt(1)=='s') v2/=0.2388444444;
double pp[]=property_SI(s,v1,v2);
pp[0]*=14.504;
pp[1]=pp[1]*1.8+32;
pp[2]*=16.01877892;
pp[3]*=0.42992;
pp[4]*=0.42992;
pp[5]*=0.2388444444;
pp[7]=1.0/pp[2];
return pp;
}

public double[] property_SI(String s, double v1, double v2)
{
//cases of String s
//    tv tp th tu ts tx
// pv    pt ph pu ps px
// vp    vt
double[] r=new double[8];
char s1=s.charAt(0);
char s2=s.charAt(1);
r[6]=x_tv(v1,v2);
if(s.equals("tv"))
{
 r[6]=x_tv(v1,v2);
 if(r[6]>=0 && r[6]<=1)
 { r[0]=Ps(v1);
   double vsl=vsl_t(v1);
   double vsg=vsg_t(v1);
   r[3]=h(v1,vsl)*(1-r[6])+ h(v1,vsl)*r[6];
   r[4]=u(v1,vsl)*(1-r[6])+ u(v1,vsl)*r[6];
   r[5]=s(v1,vsl)*(1-r[6])+ s(v1,vsl)*r[6];
 }
 else
 { r[0]=P(v1,v2);
   r[3]=h(v1,v2);
   r[4]=h(v1,v2);
   r[5]=s(v1,v2);
 }
   r[1]=v1;
   r[2]=v2;
} //end of tv case
else if(s.equals("pv"))
{
 r[6]=x_pv(v1,v2);
 if(r[6]>=0 && r[6]<=1)
 { r[1]=Ts(v1);
   double vsl=vsl_t(r[1]);
   double vsg=vsg_t(r[1]);
   r[3]=h(r[1],vsl)*(1-r[6])+ h(r[1],vsl)*r[6];
   r[4]=u(r[1],vsl)*(1-r[6])+ u(r[1],vsl)*r[6];
   r[5]=s(r[1],vsl)*(1-r[6])+ s(r[1],vsl)*r[6];
 }
 else
 {
   r[1]=t_pv(v1,v2);
   r[3]=h(r[1],v2);
   r[4]=u(r[1],v2);
   r[5]=s(r[1],v2);
 }
  r[0]=v1;
  r[2]=v2;
} //end of pv case
else if(s.equals("tp"))
{
   r[0]=v2;
   r[1]=v1;
   r[2]=v_tp(v1,v2);
   r[3]=h(r[1],r[2]);
   r[4]=u(r[1],r[2]);
   r[5]=s(r[1],r[2]);
   r[6]=x_tv(r[1],r[2]);
} //end of tp case
else if(s.equals("th"))
{
 double hsl=h(v1,vsl_t(v1));
 double hsg=h(v1,vsg_t(v1));
 double vth=v_th(v1,v2);
   r[0]=P(v1,vth);
   r[1]=v1;
   r[2]=vth;
   r[3]=v2;
   r[4]=u(v1,vth);
   r[5]=s(v1,vth);
   r[6]=x_tv(v1,r[2]);
 if(r[6]>=0 && r[6]<=1)
 { r[0]=Ps(v1);
   r[1]=v1;
   double vsl=vsl_t(v1);
   double vsg=vsg_t(v1);
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=v2;
   r[4]=h(v1,vsl)*(1-r[6])+ h(v1,vsg)*r[6];
   r[5]=s(v1,vsl)*(1-r[6])+ s(v1,vsg)*r[6];
  }
} //end of th case
else if(s.equals("tu"))
{
   double vtu=v_tu(v1,v2);
   r[0]=P(v1,vtu);
   r[1]=v1;
   r[2]=vtu;
   r[3]=h(v1,vtu);
   r[4]=v2;
   r[5]=s(v1,vtu);
   r[6]=x_tv(v1,r[2]);
   if(r[6]>=0 && r[6]<=1)
   {
      r[0]=Ps(v1);
      r[1]=v1;
      double vsl=vsl_t(v1);
      double vsg=vsg_t(v1);
      r[2]=vsl*(1-r[6])+ vsg*r[6];
      r[3]=h(v1,vsl)*(1-r[6])+ h(v1,vsg)*r[6];
      r[4]=v2;
      r[5]=s(v1,vsl)*(1-r[6])+ s(v1,vsg)*r[6];
  }
} // end of tu case
else if(s.equals("ts"))
{
   double vts=v_ts(v1,v2);
   r[0]=P(v1,vts);
   r[1]=v1;
   r[2]=v_ts(v1,v2);
   r[3]=h(v1,vts);
   r[4]=u(v1,vts);
   r[5]=v2;
   r[6]=x_tv(v1,r[2]);
 if(r[6]>=0 && r[6]<=1)
 { r[0]=Ps(v1);
   double vsl=vsl_t(v1);
   double vsg=vsg_t(v1);
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=h(v1,vsl)*(1-r[6])+ h(v1,vsg)*r[6];
   r[4]=u(v1,vsl)*(1-r[6])+ u(v1,vsg)*r[6];
 }

} //end of ts case
else if(s.equals("tx"))
{
   r[6]=v2;
   r[0]=Ps(v1);
   r[1]=v1;
   double vsl=vsl_t(v1);
   double vsg=vsg_t(v1);
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=h(v1,vsl)*(1-r[6])+ h(v1,vsg)*r[6];
   r[4]=u(v1,vsl)*(1-r[6])+ u(v1,vsg)*r[6];
   r[5]=s(v1,vsl)*(1-r[6])+ s(v1,vsg)*r[6];
} //end of tx case
else if(s.equals("vt"))
{
 r=property_SI("tv",v2,v1);
} //end of vt case
else if(s.equals("vp"))
{
 r=property_SI("pv",v2,v1);
} //end of vp case
else if(s.equals("pt"))
{
 r=property_SI("tp",v2,v1);
} //end of tp case
else if(s.equals("ph"))
{
 double ts=Ts(v1);
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 double hsl=h(ts,vsl);
 double hsg=h(ts,vsg);
 double vph=v_ph(v1,v2);
   r[0]=v1;
   r[1]=t_ph(v1,v2);
   r[2]=vph;
   r[3]=v2;
   r[4]=u(r[1],vph);
   r[5]=s(r[1],vph);
   if( v1 < Pc)
   { r[6]=(v2-hsl)/(hsg-hsl);}
   else r[6]=3;
 if(r[6]>=0 && r[6]<=1)
 { r[0]=v1;
   r[1]=ts;
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=v2;
   r[4]=u(ts,vsl)*(1-r[6])+ u(ts,vsg)*r[6];
   r[5]=s(ts,vsl)*(1-r[6])+ s(ts,vsg)*r[6];
  }
} //end of ph case
else if(s.equals("pu"))
{
 double ts=Ts(v1);
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 double usl=u(ts,vsl);
 double usg=u(ts,vsg);
 double vpu=v_pu(v1,v2);
   r[0]=v1;
   r[1]=t_pu(v1,v2);
   r[2]=vpu;
   r[3]=h(r[1],r[2]);
   r[4]=v2;
   r[5]=s(r[1],r[2]);
   if( v1 < Pc)
   { r[6]=(v2-usl)/(usg-usl);}
   else r[6]=3;
 if(r[6]>=0 && r[6]<=1)
 { r[0]=v1;
   r[1]=ts;
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=h(ts,vsl)*(1-r[6])+ h(ts,vsg)*r[6];
   r[4]=v2;
   r[5]=s(ts,vsl)*(1-r[6])+ s(ts,vsg)*r[6];
  }
} //end of pu case
else if(s.equals("ps"))
{
 double ts=Ts(v1);
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 double ssl=s(ts,vsl);
 double ssg=s(ts,vsg);
 double vps=v_ps(v1,v2);
   r[0]=v1;
   r[1]=t_ps(v1,v2);
   r[2]=vps;
   r[3]=h(r[1],r[2]);
   r[4]=u(r[1],r[2]);
   r[5]=v2;
   if( v1 < Pc)
   { r[6]=(v2-ssl)/(ssg-ssl);}
   else r[6]=3;
 if(r[6]>=0 && r[6]<=1)
 { r[0]=v1;
   r[1]=ts;
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=h(ts,vsl)*(1-r[6])+ h(ts,vsg)*r[6];
   r[4]=u(ts,vsl)*(1-r[6])+ u(ts,vsg)*r[6];
   r[5]=v2;
  }
} //end of ps case
else if(s.equals("px"))
{
   r[6]=v2;
   double ts=Ts(v1);
   r[0]=v1;
   r[1]=ts;
   double vsl=vsl_t(ts);
   double vsg=vsg_t(ts);
   r[2]=vsl*(1-r[6])+ vsg*r[6];;
   r[3]=h(ts,vsl)*(1-r[6])+ h(ts,vsg)*r[6];;
   r[4]=u(ts,vsl)*(1-r[6])+ u(ts,vsg)*r[6];
   r[5]=s(ts,vsl)*(1-r[6])+ s(ts,vsg)*r[6];
} //end of px case
r[7]=1.0/r[2]*1.00006493;
return r;
}

public double[] property(String s, double v1, double v2)
{
if(unit.equals(u[0])) return property_SI(s,v1,v2);
else                 return property_EN(s,v1,v2);
}


public String[][] toString(String s, double v1, double v2)
{
  String s1[][]=new String[9][3];
   s1[0][0]="P, Basınç          ";
   s1[1][0]="T, Sıcaklık        ";
   s1[2][0]="v, Özgül hacim     ";
   s1[3][0]="h, entalpi         ";
   s1[4][0]="u, İç Enerji       ";
   s1[5][0]="s, entropi         ";
   s1[6][0]="x, buhar doygunluk derecesi";
   s1[7][0]="   yoğunluk        ";
   s1[8][0]="Faz                ";
   if(unit.equals("SI"))
   {
   s1[0][2]=" bar          ";
   s1[1][2]=" derece C     ";
   s1[2][2]=" m^3/kg       ";
   s1[3][2]=" KJ/kg        ";
   s1[4][2]=" KJ/kg        ";
   s1[5][2]=" KJ/kg        ";
   s1[6][2]=" kg buhar/kg  ";
   s1[7][2]=" kg/m^3       ";
   s1[8][2]=" ";
   }
   else
   {
   s1[0][2]=" lbf/in^2       ";
   s1[1][2]=" derece F       ";
   s1[2][2]=" ft^3/lbm       ";
   s1[3][2]=" BTU/lbm        ";
   s1[4][2]=" BTU/lbm        ";
   s1[5][2]=" BTU/lbm        ";
   s1[6][2]=" lbm buhar/lbm  ";
   s1[7][2]=" lbm/ft^3       ";
   s1[8][2]=" ";
   }
double pp[]=property(s,v1,v2);
for(int i=0;i<7;i++)
   {s1[i][1]=""+pp[i];}
    s1[7][1]=""+pp[7];
    s1[8][1]=phase(pp[6]);
return s1;
}

public String toString1(String s, double v1, double v2)
{
   String[] s1=new String[7];
   String[] s2=new String[7];
   String s3="";
   double pp[];
   s1[0]="P :";
   s1[1]="T :";
   s1[2]="v :";
   s1[3]="h :";
   s1[4]="u :";
   s1[5]="s :";
   s1[6]="x :";
   s1[7]="  :";
   if(unit.equals("SI"))
   {
   s2[0]=" bar           ";
   s2[1]=" derece C      ";
   s2[2]=" m^3/kg        ";
   s2[3]=" KJ/kg         ";
   s2[4]=" KJ/kg         ";
   s2[5]=" KJ/kg         ";
   s2[6]=" kg buh/kg     ";
   s2[7]=" kg/m^3        ";
   }
   else
   {
   s2[0]=" lbf/in^2       ";
   s2[1]=" deg F          ";
   s2[2]=" ft^3/lbm       ";
   s2[3]=" BTU/lbm        ";
   s2[4]=" BTU/lbm        ";
   s2[5]=" BTU/lbm        ";
   s2[6]=" lbm buhar/lbm  ";
   s2[7]=" lbm/ft^3       ";
   }
pp=property(s,v1,v2);
s3+=phase(pp[6])+"\n";

for(int i=1;i<9;i++)
{
s3+=s1[i-1]+s2[i-1]+pp[i-1]+"\n";
}
return s3;
}

public String toString()
{
   String[] s1=new String[7];
   String[] s2=new String[7];
   String s3="";
   double pp[];
   s1[0]="P :";
   s1[1]="T :";
   s1[2]="v :";
   s1[3]="h :";
   s1[4]="u :";
   s1[5]="s :";
   s1[6]="x :";
   s1[7]="  :";
   if(unit.equals("SI"))
   {
   s2[0]=" bar          ";
   s2[1]=" derece C         ";
   s2[2]=" m^3/kg        ";
   s2[3]=" KJ/kg         ";
   s2[4]=" KJ/kg         ";
   s2[5]=" KJ/kg K       ";
   s2[6]=" kg buhar/kg   ";
   s2[7]=" kg/m^3        ";
   }
   else
   {
   s2[0]=" lbf/in^2       ";
   s2[1]=" derece F          ";
   s2[2]=" ft^3/lbm       ";
   s2[3]=" BTU/lbm        ";
   s2[4]=" BTU/lbm        ";
   s2[5]=" BTU/lbm        ";
   s2[6]=" lbm buhar/lbm  ";
   s2[7]=" lbm/ft^3       ";
   }
for(int i=0;i<8;i++)
{
s3+=s1[i]+s2[i]+"\n";
}
return s3;
}
   public double v_tpPR(String s,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.
        double T=t+273.15;
        double Omega=0.334;
        double R=8314.5;
        double u=2;
        double w=-1;
        double fw=0.37464+1.54226*Omega-0.26992*Omega*Omega;
        double Tr=T/Tc;
        double b=0.07780*R*Tc/(Pc*1e5);
        double a=0.45724*R*R*Tc*Tc/(Pc*1e5)*
        Math.pow(1+fw*(1.0-Math.pow(Tr,0.5)),2);
        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;
        double zc[][]=Matrix.poly_roots(cof);
        double nmin=1e99;
        double nmax=-1e99;
        for(int i=0;i<3;i++)
        {
        if((zc[0][i]>=0.0 && zc[0][i]<=1.0) && (zc[1][i]>=0.0 && zc[1][i]<1e-10))
        { if(zc[0][i]<nmin) nmin=zc[0][i];
          if(zc[0][i]>nmax) nmax=zc[0][i];
        }
        }
     double y,zz;
     if(s.equals("l")) zz=nmin;
     else zz=nmax;;
     y=zz*R/18.015*T/(P*1e5);
     return y;
   }

}