//======================================================
// 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  steam class
// class steam calculates thermophysical properties of
// water and steam. It utilise Keenan & Keyes equation of state
// DATA FILE DEFINATION
// datas are written directly inside of the class as array structure
// 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

//============================================================
//  VARIABLE IDENTIFICATION
//============================================================
//

//============================================================
//  METHOD IDENTIFICATION
//============================================================
//


public class steam2
{
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}
              };
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.6151;  //bar cm^3/g K
final double R2=0.46151; // 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

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 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 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 property(String s, double v1, double v2)
{
//several properties of steam
//s values
// v_tp : v(t,p)
// t_pv : t(p,v)
// fi_tv: fi(t,v)

//bisection method to find roots of a function Ps(Ts)-ps
//function should be supplied
// 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;
double xu;
double test;
double xr;
double es,ea;
double fxl;
double fxr;
xl=1.0001e-3;
xu=726.025;
int maxit=500,iter=0;
es=0.0000000001;
ea=1.1*es;
if(s.equals("fi_tv"))
{return fi(v1,v2);}
else if(s.equals("p_tv"))
{return P(v1,v2);}
else if(s.equals("p_tu"))
{return P(v1,property("v_tu",v1,v2));}
else if(s.equals("p_th"))
{return P(v1,property("v_th",v1,v2));}
else if(s.equals("p_ts"))
{return P(v1,property("v_ts",v1,v2));}
else if(s.equals("h_tu"))
{return h(v1,property("v_tu",v1,v2));}
else if(s.equals("h_ts"))
{return h(v1,property("v_ts",v1,v2));}
else if(s.equals("u_tv"))
{return u(v1,v2);}
else if(s.equals("h_tv"))
{return h(v1,v2);}
else if(s.equals("s_tv"))
{return s(v1,v2);}
else if(s.equals("u_tp"))
{return u(v1,property("v_tp",v1,v2));}
else if(s.equals("h_tp"))
{return h(v1,property("v_tp",v1,v2));}
else if(s.equals("s_tp"))
{return s(v1,property("v_tp",v1,v2));}
else if(s.equals("fi_tp"))
{return fi(v1,property("v_tp",v1,v2));}
else if(s.equals("v_pu"))
{return property("v_tu",property("t_pu",v1,v2),v2);}
else if(s.equals("v_ph"))
{return property("v_th",property("t_ph",v1,v2),v2);}
else if(s.equals("v_ps"))
{return property("v_ts",property("t_ps",v1,v2),v2);}
else if(s.equals("v_tx"))
{return property("vsl_t",v1)*(1-v2)+property("vslg_t",v1)*v2;}
else if(s.equals("v_px"))
{return property("vsl_p",v1)*(1-v2)+property("vslg_p",v1)*v2;}
else if(s.equals("u_tx"))
{return property("usl_t",v1)*(1-v2)+property("uslg_t",v1)*v2;}
else if(s.equals("u_px"))
{return property("usl_p",v1)*(1-v2)+property("uslg_p",v1)*v2;}
else if(s.equals("h_tx"))
{return property("hsl_t",v1)*(1-v2)+property("hslg_t",v1)*v2;}
else if(s.equals("h_px"))
{return property("hsl_p",v1)*(1-v2)+property("hslg_p",v1)*v2;}
else if(s.equals("s_tx"))
{return property("ssl_t",v1)*(1-v2)+property("sslg_t",v1)*v2;}
else if(s.equals("s_px"))
{return property("ssl_p",v1)*(1-v2)+property("sslg_p",v1)*v2;}
else if(s.equals("x_pv"))
{return x_pv(v1,v2);}
else if(s.equals("x_tv"))
{return x_tv(v1,v2);}
else if(s.equals("x_tp"))
{return x_tv(v1,property("v_tp",v1,v2));}
else if(s.equals("x_th"))
{return x_tv(v1,property("v_th",v1,v2));}
else if(s.equals("x_tu"))
{return x_tv(v1,property("v_tu",v1,v2));}
else if(s.equals("x_ts"))
{return x_tv(Ps(v1),property("v_ts",v1,v2));}
else if(s.equals("x_ph"))
{return x_pv(v1,property("v_ph",v1,v2));}
else if(s.equals("x_pu"))
{return x_pv(v1,property("v_pu",v1,v2));}
else if(s.equals("x_ps"))
{return x_pv(v1,property("v_ps",v1,v2));}
else if(s.equals("v_tp"))
{ return v_tp(v1,v2); }
else if(s.equals("v_th"))
{ return v_th(v1,v2); }
else if(s.equals("v_ts"))
{ return v_ts(v1,v2); }
else if(s.equals("v_tx"))
{ return v_tx(v1,v2); }
else if(s.charAt(0)=='v')
{ xl=1.0001e-3; xu=726.025; } //end of v
else if(s.charAt(0)=='t')
{
xl=0.0;
xu=1300;
}
else
{
xl=1.0001e-3;
xu=726.025;
}
xr=(xl+xu)/2.0;
fxl= P(v1,xl)-v2;
fxr= P(v1,xr)-v2;
while((ea>es)&&(iter<maxit))
{
xr=(xl+xu)/2.0;
iter++;
if((xl+xu)!=0)
  { ea=Math.abs((xu-xl)/(xl+xu))*100;}
else if(s.equals("v_ts"))
{
fxl= s(v1,xl)-v2;
fxr= s(v1,xr)-v2;
// v1 = t
// v2 = s
}
else if(s.equals("t_pv"))
{
fxl=P(xl,v2)-v1;
fxr=P(xr,v2)-v1;
}
else if(s.equals("t_pu"))
{
double voll=property("v_tp",xl,v1);
double volr=property("v_tp",xr,v1);
fxl=v2-u(xl,voll);
fxr=v2-u(xr,volr);
}
else if(s.equals("t_ph"))
{
double voll=property("v_tp",xl,v1);
double volr=property("v_tp",xr,v1);
fxl=v2-h(xl,voll);
fxr=v2-h(xr,volr);
}
else if(s.equals("t_ps"))
{
double voll=property("v_ts",xl,v2);
double volr=property("v_ts",xr,v2);
fxl=P(xl,voll)-v1;
fxr=P(xr,volr)-v1;
}
else if(s.equals("t_vh"))
{
fxl=h(v1,xl)-v2;
fxr=h(v1,xr)-v2;
}
else if(s.equals("t_vu"))
{
fxl=u(v1,xl)-v2;
fxr=h(v1,xr)-v2;
}
else if(s.equals("t_vs"))
{
fxl=s(v1,xl)-v2;
fxr=s(v1,xr)-v2;
}
else
{
fxl = P(v1,xl)-v2;
fxr = P(v1,xr)-v2;
}

test= fxl*fxr;
if(fxl==0.0) return xl;
if(test==0.0)       ea=0;
else if(test<0.0)  xu=xr;
else               { xl=xr; }
// System.out.println("fxl="+fxl+"xl= "+xl+"fxr="+fxr+"xr="+xr+"xu="+xu);
} //end of while
return xr;
}

double property(String s, double v1)
{
if(s.equals("ps_t"))
{return Ps(v1);}
else if(s.equals("ts_p"))
{ return Ts(v1); }
else if(s.equals("vsl_t"))
{return vsl_t(v1); }
else if(s.equals("vsl_p"))
{return vsl_t(Ts(v1)); }
else if(s.equals("vsg_t"))
{return vsg_t(v1); }
else if(s.equals("vslg_t"))
{return vsg_t(v1)-vsl_t(v1);}
else if(s.equals("vsg_p"))
{return vsl_t(Ts(v1)); }
else if(s.equals("vslg_p"))
{return vsg_t(Ts(v1))-vsl_t(Ts(v1)); }
else if (s.equals("hsl_t"))
{return h(v1,vsl_t(v1));}
else if (s.equals("hsl_p"))
{double ts=Ts(v1);return h(ts,vsl_t(ts));}
else if (s.equals("hsg_t"))
{return h(v1,vsg_t(v1));}
else if (s.equals("hsg_p"))
{double ts=Ts(v1);return h(ts,vsg_t(ts));}
else if (s.equals("hslg_t"))
{return h(v1,vsg_t(v1))-h(v1,vsl_t(v1));}
else if (s.equals("hslg_p"))
{double ts=Ts(v1);return h(ts,vsg_t(ts))-h(ts,vsl_t(ts));}
else if (s.equals("usl_t"))
{return u(v1,vsl_t(v1));}
else if (s.equals("usl_p"))
{double ts=Ts(v1);return u(ts,vsl_t(ts));}
else if (s.equals("usg_t"))
{return u(v1,vsg_t(v1));}
else if (s.equals("usg_p"))
{double ts=Ts(v1);return u(ts,vsg_t(ts));}
else if (s.equals("uslg_t"))
{return u(v1,vsg_t(v1))-u(v1,vsl_t(v1));}
else if (s.equals("uslg_p"))
{double ts=Ts(v1);return u(ts,vsg_t(ts))-u(v1,vsl_t(ts));}
else if (s.equals("ssl_t"))
{return s(v1,vsl_t(v1));}
else if (s.equals("ssl_p"))
{double ts=Ts(v1);return s(ts,vsl_t(ts));}
else if (s.equals("ssg_t"))
{return s(v1,vsg_t(v1));}
else if (s.equals("ssg_p"))
{double ts=Ts(v1);return s(ts,vsg_t(ts));}
else if (s.equals("sslg_t"))
{return s(v1,vsg_t(v1))-s(v1,vsl_t(v1));}
else if (s.equals("sslg_p"))
{double ts=Ts(v1);return s(ts,vsg_t(ts))-s(ts,vsl_t(ts));}
return 0;
}

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

public double T(double p,double v)
{
return property("t_pv",p,v);
}

// 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;
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>=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 ts=Ts(p);
double x=1.0002;
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=8314.5/18.01*t/(p*1.0e5)*1e3; }
  else if(p < Pc)
  {  if(t <= ts) x=vsl_t(t)*1e3;
     else if (t>=ts) x=vsg_t(t)*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)
{
double ts=t;
double vsl=vsl_t(t);
double vsg=vsg_t(t);
double hsl=h(ts,vsl);
double hsg=h(ts,vsg);
double x=1.0002;
int nmax=100;
double dx=1.0e-3;
double dx1=2.0e-3;
double x1m;
double tolerance=1.0e-10;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
  if( t >= Tc)
  { x=5.0; }
  else if(t < Tc)
  {  if(hi < hsl) x=vsl*1e3;
     else if (hi >hsg) x=vsg*1e3;
     //else if(hi >= hsl && hi<=hsg) {double xx=(hi-hsl)/(hsg-hsl); return vsl*(1-xx)+vsg;}
  }
  for(int i=0;i<nmax;i++)
  {
  fx=h(t,x*1e-3)-hi;
  fxm=h(t,(x-dx)*1e-3)-hi;
  fxm1=h(t,(x-dx1)*1e-3)-hi;
  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_tu(double t, double ui)
{
double ts=t;
double vsl=vsl_t(t);
double vsg=vsg_t(t);
double usl=u(ts,vsl);
double usg=u(ts,vsg);
double x=1.0002;
int nmax=100;
double dx=1.0e-3;
double dx1=2.0e-3;
double x1m;
double tolerance=1.0e-10;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
  if( t >= Tc)
  { x=5.0; }
  else if(t < Tc)
  {  if(ui < usl) x=vsl*1e3;
     else if (ui >usg) x=vsg*1e3;
     //else if(ui >= usl && ui<=usg) {double xx=(ui-usl)/(usg-usl); return vsl*(1-xx)+vsg;}
  }
  for(int i=0;i<nmax;i++)
  {
  fx=u(t,x*1e-3)-ui;
  fxm=u(t,(x-dx)*1e-3)-ui;
  fxm1=u(t,(x-dx1)*1e-3)-ui;
  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_ts(double t, double si)
{
double ts=t;
double vsl=vsl_t(t);
double vsg=vsg_t(t);
double ssl=s(ts,vsl);
double ssg=s(ts,vsg);
double x=1.0002;
int nmax=100;
double dx=1.0e-3;
double dx1=2.0e-3;
double x1m;
double tolerance=1.0e-10;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
  if( t >= Tc)
  { x=5.0; }
  else if(t < Tc)
  {  if(si < ssl) x=vsl*1e3;
     else if (si >ssg) x=vsg*1e3;
     // else if(si >= ssl && si<=ssg) {double xx=(si-ssl)/(ssg-ssl); return vsl*(1-xx)+vsg;}
  }
  for(int i=0;i<nmax;i++)
  {
  fx=s(t,x*1e-3)-si;
  fxm=s(t,(x-dx)*1e-3)-si;
  fxm1=s(t,(x-dx1)*1e-3)-si;
  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_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 "Supercritical fluid";
else
 {
 double ts=t;
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 if(v<vsl)
   return "liquid";
 else if(v==vsl)
   return "saturated liquid";
 else if(v==vsg)
   return "saturated vapor";
 else if(v>vsl && v<vsg)
   return "saturated liquid-vapor mixture";
 else // (v>vsg)
   return "superheated vapor";
 }
}

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.01)/(1.0+B(ro)); }
  else if(p < Pc)
  {  x=ts; }
  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 v_ph(double p, double hi)
{
double ts=Ts(p);
double vsl=vsl_t(ts);
double vsg=vsg_t(ts);
double hsl=h(ts,vsl);
double hsg=h(ts,vsg);
double x=5.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=5.0; }
  else if(p < Pc)
  {
  if( hi <= hsl) x=vsl*1e3;
  else if( hi >= hsg) x=vsg*1e-3;
  }
  for(int i=0;i<nmax;i++)
  {
  fx=  u(x      , v_tu(x     , p))-hi;
  fxm= u((x-dx) , v_tu((x-dx), p))-hi;
  fxm1=u((x-dx1), v_tu((x-dx1),p))-hi;
  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 t_pu(double p, double ui)
{
double ts=Ts(p);
double usl=vsl_t(ts);
double usg=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=200.0; }
  else if(p < Pc)
  {  x=ts; }
  for(int i=0;i<nmax;i++)
  {
  fx=  u(x      , v_tu(x     , p))-ui;
  fxm= u((x-dx) , v_tu((x-dx), p))-ui;
  fxm1=u((x-dx1), v_tu((x-dx1),p))-ui;
  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 -1;
 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 2;
 }
}
public double x_pv(double p, double v)
{
 return x_tv(Ts(p),v);
}

public double[] proper(String s, double v1, double v2)
{
//cases of s
//    tv tp th tu ts tx
// pv       ph pu ps px
//          vh vu vs vx
double[] r=new double[7];
char s1=s.charAt(0);
char s2=s.charAt(1);
r[6]=property(("x_"+s),v1,v2);
if(s.equals("tv"))
{
 r[6]=property("x_"+s,v1,v2);
 if(r[6]>=0 && r[6]<=1)
 { r[0]=Ps(v1);
   r[3]=property("hsl_t",v1)*(1-r[6])+ property("hsg_t",v1)*r[6];
   r[4]=property("usl_t",v1)*(1-r[6])+ property("usg_t",v1)*r[6];
   r[5]=property("ssl_t",v1)*(1-r[6])+ property("ssg_t",v1)*r[6];
 }
 else
 { r[0]=property("p_tv",v1,v2);
   r[3]=property("h_tv",v1,v2);
   r[4]=property("u_tv",v1,v2);
   r[5]=property("s_tv",v1,v2);
 }
   r[1]=v1;
   r[2]=v2;

} //end of tv case
else if(s.equals("pv"))
{
 r[6]=property("x_"+s,v1,v2);
 if(r[6]>=0 && r[6]<=1)
 { r[1]=Ts(v1);
   r[3]=property("hsl_t",r[1])*(1-r[6])+ property("hsg_t",r[1])*r[6];
   r[4]=property("usl_t",r[1])*(1-r[6])+ property("usg_t",r[1])*r[6];
   r[5]=property("ssl_t",r[1])*(1-r[6])+ property("ssg_t",r[1])*r[6];
 }
 else
 {
   r[1]=property("t_pv",v1,v2);
   r[3]=property("h_tv",r[1],v2);
   r[4]=property("u_tv",r[1],v2);
   r[5]=property("s_tv",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]=property("v_tp",v1,v2);
   r[3]=property("h_tv",r[1],r[2]);
   r[4]=property("u_tv",r[1],r[2]);
   r[5]=property("s_tv",r[1],r[2]);
   r[6]=property("x_tv",r[1],r[2]);
} //end of tp case
else if(s.equals("th"))
{
 double hsl=property("hsl_t",v1);
 double hsg=property("hsl_t",v1);
 r[6]=(v2-hsl)/(hsg-hsl);
 if(r[6]>=0 && r[6]<=1)
 { r[0]=Ps(v1);
   r[2]=property("vsl_t",v1)*(1-r[6])+ property("vsg_t",v1)*r[6];
   r[4]=property("usl_t",v1)*(1-r[6])+ property("usg_t",v1)*r[6];
   r[5]=property("ssl_t",v1)*(1-r[6])+ property("ssg_t",v1)*r[6];
  }
 else
 { r[0]=property("p_th",v1,v2);
   r[2]=property("v_th",v1,v2);
   r[4]=property("u_tv",v1,r[2]);
   r[5]=property("s_tv",v1,r[2]);
 }
   r[1]=v1;
   r[3]=v2;
   r[6]=property("x_tv",v1,r[2]);
} //end of th case
else if(s.equals("tu"))
{
 double usl=property("usl_t",v1);
 double usg=property("usl_t",v1);
 r[6]=(v2-usl)/(usg-usl);
 if(r[6]>=0 && r[6]<=1)
 { r[0]=Ps(v1);
   r[2]=property("vsl_t",v1)*(1-r[6])+ property("vsg_t",v1)*r[6];
   r[3]=property("hsl_t",v1)*(1-r[6])+ property("hsg_t",v1)*r[6];
   r[5]=property("ssl_t",v1)*(1-r[6])+ property("ssg_t",v1)*r[6];
  }
 else
 { double vtu=v_tu(v1,v2);
   r[0]=P(v1,vtu);
   r[2]=vtu;
   r[3]=h(v1,vtu);
   r[5]=s(v1,vtu);
 }
   r[1]=v1;
   r[4]=v2;
   r[6]=x_tv(v1,r[2]);
} // end of tu case
else if(s.equals("ts"))
{
 double ssl=property("ssl_t",v1);
 double ssg=property("ssl_t",v1);
 r[6]=(v2-ssl)/(ssg-ssl);
 if(r[6]>=0 && r[6]<=1)
 { r[0]=Ps(v1);
   r[2]=property("vsl_t",v1)*(1-r[6])+ property("vsg_t",v1)*r[6];
   r[3]=property("hsl_t",v1)*(1-r[6])+ property("hsg_t",v1)*r[6];
   r[4]=property("usl_t",v1)*(1-r[6])+ property("usg_t",v1)*r[6];
 }
 else
 { double vts=v_ts(v1,v2);
   r[0]=P(v1,vts);
   r[2]=v_ts(v1,v2);
   r[3]=h(v1,vts);
   r[4]=u(v1,vts);
 }
   r[1]=v1;
   r[5]=v2;
   r[6]=property("x_tv",v1,r[2]);
} //end of ts case
else if(s.equals("tx"))
{
   r[6]=v2;
   r[0]=Ps(v1);
   r[1]=v1;
   r[2]=property("vsl_t",v1)*(1-r[6])+ property("vsg_t",v1)*r[6];
   r[3]=property("hsl_t",v1)*(1-r[6])+ property("hsg_t",v1)*r[6];
   r[4]=property("usl_t",v1)*(1-r[6])+ property("usg_t",v1)*r[6];
   r[5]=property("ssl_t",v1)*(1-r[6])+ property("ssg_t",v1)*r[6];
} //end of tx case
else if(s.equals("vt"))
{
 r=proper("tv",v2,v1);
} //end of vt case
else if(s.equals("vp"))
{
 r=proper("pv",v2,v1);
} //end of vp case
else if(s.equals("pt"))
{
 r=proper("tp",v2,v1);
} //end of tp case
else if(s.equals("ph"))
{
 double ts=Ts(v1);
 double hsl=property("hsl_t",ts);
 double hsg=property("hsl_t",ts);
 r[6]=(v2-hsl)/(hsg-hsl);
 if(r[6]>=0 && r[6]<=1)
 { r[0]=v1;
   r[1]=Ts(v1);
   r[2]=property("vsl_t",ts)*(1-r[6])+ property("vsg_t",ts)*r[6];
   r[4]=property("usl_t",ts)*(1-r[6])+ property("usg_t",ts)*r[6];
   r[5]=property("ssl_t",ts)*(1-r[6])+ property("ssg_t",ts)*r[6];
  }
 else
 { r[1]=property("t_ph",v1,v2);
   r[2]=property("v_ph",v1,v2);
   r[4]=property("u_ph",v1,v2);
   r[5]=property("s_ph",v1,v2);
 }
   r[0]=v1;
   r[3]=v2;
   r[6]=property("x_tv",r[1],r[2]);
} //end of th case
else if(s.equals("pu"))
{
 double ts=Ts(v1);
 double usl=property("usl_t",ts);
 double usg=property("usl_t",ts);
 r[6]=(v2-usl)/(usg-usl);
 if(r[6]>=0 && r[6]<=1)
 { r[1]=ts;
   r[2]=property("vsl_t",ts)*(1-r[6])+ property("vsg_t",ts)*r[6];;
   r[3]=property("hsl_t",ts)*(1-r[6])+ property("hsg_t",ts)*r[6];
   r[5]=property("ssl_t",ts)*(1-r[6])+ property("ssg_t",ts)*r[6];
 }
 else
 { r[1]=property("t_ph",v1,v2);
   r[2]=property("v_ph",v1,v2);
   r[3]=property("h_ph",v1,v2);
   r[5]=property("s_ph",v1,v2);
 }
   r[0]=v1;
   r[4]=v2;
   r[6]=property("x_tv",r[1],r[2]);
} //end of tu case
else if(s.equals("ps"))
{
 double ts=Ts(v1);
 double ssl=property("ssl_t",ts);
 double ssg=property("ssl_t",ts);
 r[6]=(v2-ssl)/(ssg-ssl);
 if(r[6]>=0 && r[6]<=1)
 { r[1]=ts;
   r[2]=property("vsl_t",ts)*(1-r[6])+ property("vsg_t",ts)*r[6];;
   r[3]=property("hsl_t",ts)*(1-r[6])+ property("hsg_t",ts)*r[6];;
   r[4]=property("usl_t",ts)*(1-r[6])+ property("usg_t",ts)*r[6];
 }
 else
 { r[1]=property("t_ps",v1,v2);
   r[2]=property("v_ps",v1,v2);
   r[3]=v2;
   r[4]=property("u_ps",v1,v2);
 }
   r[1]=v1;
   r[5]=v2;
   r[6]=property("x_tv",r[1],r[2]);
} //end of ts case
else if(s.equals("px"))
{
   r[6]=v2;
   double ts=Ts(v1);
   r[0]=v1;
   r[1]=ts;
   r[2]=property("vsl_t",ts)*(1-r[6])+ property("vsg_t",ts)*r[6];;
   r[3]=property("hsl_t",ts)*(1-r[6])+ property("hsg_t",ts)*r[6];;
   r[4]=property("usl_t",ts)*(1-r[6])+ property("usg_t",ts)*r[6];
   r[5]=property("ssl_t",ts)*(1-r[6])+ property("ssg_t",ts)*r[6];
} //end of px case
return r;
}
}
