//======================================================
// Thermodynamic Package in java
// Class Steam Properties of amonia
// Dr. Turhan Coban
// TUBITAK Ulusal Metroloji Enstitüsü
// email : Turhan.Coban@ume.tubitak.gov.tr
// 
// File Name   :amonia.java
// This file contains the amonia 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 :1.2333498
//  to  : 500/T
//  F : coefficients of the saturation state
//  Pc: critical pressure kpa
//  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 ammonia
{
ref_CS3 r;
refrigerantData rd;
final double A[][]={
{-0.645302230405E+01,-0.137199267705E+02,-0.810062031570E+01,-0.488009642100E+01,-0.120287756268E+02,0.680634593000E+01},
{0.808009436769E+01,0.143569200060E+02,-0.450529767000E+02,-0.166188998570E+03,0.379089502300E+02,-0.407302083337E+02},
{0.103299488070E+01,0.558439558100E+02,0.492016650818E+03,0.173783600000E+04,-0.308749152638E+02,0.714835304160E+02},
{-0.894826463200E+01,-0.169777744140E+03,-0.123653237167E+04,-0.781216116832E+04,0.177954826914E+01,-0.389746109600E+02},
{-0.669228588200E+02,-0.175394377532E+01,0.208553371335E+03,0.213489466000E+05,0.0,0.0},
{0.247341746000E+03,0.299983915550E+03,0.450908057879E+04,-0.379808500000E+05,0.0,0.0},
{-0.306557885430E+03,0.241165511000E+02,-0.932335680000E+04,0.427240985306E+05,0.0,0.0},
{0.161791003337E+03,-0.507478070464E+03,0.813947039740E+04,-0.274587106300E+05,0.0,0.0},
{-0.278216888000E+02,0.298812917300E+03,-0.277259735200E+04,0.766892867800E+04,0.0,0.0}};
final String u[]={"SI","EN"};
String unit=u[0];
double toaj;
double roaj;
final double C[]={-0.38727270E+01,0.64463724,0.32238759E+01,-0.21376925E-02,0.86890833E-05,-0.24085149E-07,0.36893175E-10,-0.35034664E-13,0.20563027E-16,-0.68534200E-20,0.99392427E-24}; //  J/gr
final double R=4.8180;  // atm-cm^3/g K
//final double R1=4.8818385; // bar cm^3/g K
final double R2=0.488219205; // j/g K
final double Tc=405.4;  //K
final double toc=500.0/Tc;
final double B[]={-7.296510,1.618053,-1.956546,-2.114118};
final double Pc=113.04*1e2;   //kPA

final double vc=4.25531914893; //cm^3/g

public ammonia()
{
r=new ref_CS3("R717");
}


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=6;
//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 Q(double ro, double t)
{
//coefficient of amonia state equation, fi
  double q=0;
  double T=t+273.15;
  double to=500.0/T;
  for(int i=0;i<9;i++)
  {
    for(int j=0;j<6;j++)
    {
	if(j==0)
	{if(i!=0) q+=A[i][j]*Math.pow(ro,i);
	 else q+=A[i][j];
	}
	else 
    {if(i!=0) q+=A[i][j]*Math.pow(ro,i)*Math.pow((to-toc),j);
     else     q+=A[i][j]*Math.pow((to-toc),j);
    }       
    }
  }
  return q;
}

double dQdto(double ro, double t)
{
//coefficient of amonia state equation, fi
  double q=0;
  double T=t+273.15;
  double to=500.0/T;
  for(int i=0;i<9;i++)
  {
    for(int j=1;j<6;j++)
    {
	if(j==1)
	{if(i!=0) q+=A[i][j]*Math.pow(ro,i);
	 else q+=A[i][j];
	}
	else 
    {if(i!=0) q+=A[i][j]*j*Math.pow(ro,i)*Math.pow((to-toc),(j-1));
     else     q+=A[i][j]*j*Math.pow((to-toc),(j-1));
    }       
    }
  }
  return q;
}

double dQdro(double ro, double t)
{
//coefficient of amonia state equation, fi
  double q=0;
  double T=t+273.15;
  double to=500.0/T;
  for(int i=1;i<9;i++)
  {
    for(int j=0;j<6;j++)
    {
	if(j==0)
	{if(i!=1) q+=A[i][j]*i*Math.pow(ro,(i-1));
	 else q+=A[i][j];
	}
	else 
    {if(i!=1) q+=A[i][j]*i*Math.pow(ro,(i-1))*Math.pow((to-toc),j);
     else     q+=A[i][j]*Math.pow((to-toc),j);
    }       
    }
  }
  return q;
}

double dQ2drodt(double ro, double t)
{
//coefficient of amonia state equation, fi
  double q=0;
  double T=t+273.15;
  double to=500.0/T;
  for(int i=1;i<9;i++)
  {
    for(int j=1;j<6;j++)
    {
	if(j==1)
	{if(i!=1) q+=A[i][j]*i*Math.pow(ro,(i-1));
	 else q+=A[i][j];
	}
	else 
    {if(i!=1) q+=A[i][j]*i*j*Math.pow(ro,(i-1))*Math.pow((to-toc),j);
     else     q+=A[i][j]*j*Math.pow((to-toc),(j-1));
    }       
    }
  }
  return q;
}


double dQ2dro2(double ro, double t)
{
//coefficient of amonia state equation, fi
  double q=0;
  double T=t+273.15;
  double to=500.0/T;
  for(int i=2;i<9;i++)
  {
    for(int j=0;j<6;j++)
    {
	if(j==0)
	{if(i!=2) q+=A[i][j]*i*(i-1)*Math.pow(ro,(i-2));
	 else q+=A[i][j]*i;
	}
	else 
    {if(i!=2) q+=A[i][j]*i*(i-1)*Math.pow(ro,(i-2))*Math.pow((to-toc),j);
     else     q+=A[i][j]*i*Math.pow((to-toc),j);
    }       
    }
  }
  return q;
}


double A0(double t)
{
// helmholts free energy function for reference pressure
 double R=4.8180;
 double s=0;
 double T=t+273.15;  //K
 s+=R2*T*C[0]*Math.log(T);
 for(int i=1;i<11;i++)
 {
 s+=R2*T*C[i]*Math.pow(T,(i-2));
 }
 s-=R2*T*(1.0-Math.log(R*T));
 return s;
}

double dA0dt(double t)
{
// helmholts free energy function for reference pressure
 double R=4.8180;
 double s=0;
 double T=t+273.15;  //K
 s+=R2*C[0]*(Math.log(T)+1.0);
 s+=R2*Math.log(R*T);
 for(int i=1;i<11;i++)
 {
 s+=R2*(i-1)*C[i]*Math.pow(T,(i-2));
 }
 return s;
}

double Abar(double t, double v)
{ double R=4.8180;
  double T=t+273.15; //K
  double to=500.0/T;
  double ro=1/v*1e-3;
  return R2*T*(Math.log(ro)+ro*Q(ro,t));
}
	
double A(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=500.0/T;
  double ro=1/v*1e-3;
  return Abar(t,v)+A0(t);
}

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

double Pi(double t, double v)
{
// pressure p(t,v)
// T C
// v m^3/kg
double T=t+273.15; //K
double to=500.0/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);}
else return P(t,v);
}


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

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

double u(double t, double v)
{
//internal energy
double u=0;
double R=4.8180;
double T=t+273.15; //K
double to=500.0/T;
double ro=1/v*1e-3;
u+=R2*T*ro*to*(dQdto(ro,t))+A0(t)-T*dA0dt(t);
//System.out.println("R2="+R2+"to="+to+"ro="+ro+"T="+T+"u="+u); 
return u;
}

public double h(double t, double v)
{
//enthalpy
double h=0;
double R=4.8180;
double T=t+273.15; //K
double to=500/T;
double ro=1/v*1e-3;
h+=R2*T*(ro*(Q(ro,t))+ro*ro*(dQdro(ro,t))+ro*to*(dQdto(ro,t))+1.0)+A0(t)-T*dA0dt(t);
//System.out.println("R2="+R2+"to="+to+"ro="+ro+"T="+T+"h="+h); 
return h;

}

public double s(double t, double v)
{
//entropy
double s=0;
double R=4.8180;
double T=t+273.15; //K
double to=500.0/T;
double ro=1/v*1e-3;
s+=-R2*(Math.log(ro)+ro*(Q(ro,t))-ro*to*(dQdto(ro,t)))-dA0dt(t);
//System.out.println("R2="+R2+"to="+to+"ro="+ro+"T="+T+"s="+s); 
return s;
}

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

// Doyma termodinamik ?zellikleri
//doyma gaz ?zg?l hacim
//doyma s?v? ?zg?l hacim
double vf(double t)
{return 1.0/r.rd.rol(t);}
// doyma gaz ?zg?l hacim
double vg(double t)
{return 1.0/r.rd.rov(t);}


double Ps(double t)
{return r.rd.Psb(t);}

double Psb(double t)
{return r.rd.Psb(t);}

double Psd(double t)
{return r.rd.Psd(t);}
double Ts(double ps)
{return r.Tsd(ps);}

double Tsd(double ps)
{return r.Tsd(ps);}

double Tsb(double ps)
{return r.Tsb(ps);}

// doyma s?v? entalpi
double hf(double t)
{return r.rd.h_l(t);}
//doyma gaz entalpi

double hg(double t)
{return r.rd.h_v(t);}

//doyma s?v? i? enerji
double uf(double t)
{return r.u_l(t);}

//doyma gaz i? enerji

double ug(double t)
{return r.u_v(t);}

// doyma s?v? entropi
double sf(double t)
{return r.s_l(t);}
//doyma gaz entropi

double sg(double t)
{return r.rd.s_v(t);}

double vsl_t(double t)
{return r.rd.vl(t);}

double vsg_t(double t)
{return r.rd.vv(t);}

double vg_app(double t)
{return 1.0/r.rov(t);}

double v_tp2(double t,double P)
{
// 
double a,b;
  double Tc=r.Tc;
  double Td=Tsd(P);
  double Tb=Tsb(P);
//double R=rd.R;
double vc,vg,vf;
vc=r.vc;
vg=vg(t);
vf=vf(t);
double b1=(8314.5/r.M)*(t+273.15)/(P*1.0e3);
if(t>Td ) {a=vg;b=b1;} //supercritical or superheated gas
else if(t<Tb) {a=5e-5;b=vf;}  //liquid
else {a=vf;b=vg;} //saturated
double test;
double p=0;
double es,ea;
double f1,f2,f3,fp;
int maxit=500,iter=0;
double tol=1.0e-5;
es=0.0000001;
double x1=a;f1=P(t,x1)-P;
double x2=b;f2=P(t,x2)-P;
double x3=(x1+x2)/2.0;f3=P(t,x3)-P;;
if(f1==0) return x1;
else if(f2==0) return x2;
else if(f3==0) return x3; 
//if(f1*f2>0) System.out.println("No root is existed in the given region");
p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
fp=P(t,p)-P;;
ea=Math.abs(x3-p);
while((ea>es)&&(iter<maxit))
{	if(Math.abs(f3)<tol) return x3; 
    if((p<x1) && (p>x2)) 
    {p=(x1+x2)/2.0;    
     if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
     else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    }
    else
    {        
    if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
    else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
    fp=P(t,p)-P;;
    ea=Math.abs(x3-p);
    }
    iter++;
} 
//if(iter>=maxit) JOptionPane.showMessageDialog(null,"Warning : MAximum number of iteration is exceeded \n"+
//   " result may not be valid","MAXIMUM ITERATION NUMBER WARNING",JOptionPane.WARNING_MESSAGE);
return p;
}

double v_tp3(double t,double P)
{ // P kPa
  // Pc kritik bas?n? r.Pc
  double Tc=r.Tc;
  double Td=Tsd(P);
  double Tb=Tsb(P);
// brent metodu ile k?k de?erini bulur
double vc,vg,vf;
double x1=0,x2=0;
vc=r.vc;
vg=vg(t);
vf=vf(t);
double b=R2*(t+273.15)/P;
double x=v_tpPR("v",t,P)*1e3;
System.out.println("Tc="+Tc+"Td="+Td+"Tb="+Tb+"t="+t+"b="+b+"x="+x+"vg="+vg+"vf="+vf);
if(t>Td ) {x1=vg;x2=x;} //supercritical or superheated gas
else if(t<Tb) {x1=5e-5;x2=vf;}  //liquid
else {x1=vf;x2=vg;} //saturated
double test;
double es,ea;
double f1,f2,f3,fp;
int maxit=500,iter=0;
double tol=1.0e-5;
double p=0;
es=0.0001;
f1=P(t,x1)-P;
f2=P(t,x2)-P;
//System.out.println("f1="+f1+"f2="+f2+"P="+P);
double x3=(x1+x2)/2.0;f3=P(t,x3)-P;
if(f1==0) return x1;
else if(f2==0) return x2;
else if(f3==0) return x3; 
//if(f1*f2>0) System.out.println("verilen b?lgede k?k yok");
p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
fp=P(t,p)-P;
ea=Math.abs(x3-p);
while((ea>es)&&(iter<maxit))
{	if(Math.abs(f3)<tol) return x3; 
    if((p<x1) && (p>x2)) 
    {p=(x1+x2)/2.0;    
     if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
     else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    }
    else
    {        
    if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
    else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
    fp=P(t,p)-P; 
    ea=Math.abs(x3-p);
    }
    iter++;
} 
return p;
}
public double v_tp1(double t,double Pi)
{
//specific volume as a function of temperature and pressure	
// x in f(x)=0 given also df(x)/dx 
// Newton - Raphson type method 
double T=t+273.15;
double x=v_tpPR("v",t,Pi)*1e3;//(T*R2/Pi)/r.M;
int nmax=500;
double tolerance=1.0e-10;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
fx=P(t,x)-Pi;
dfx=dPdv(t,x);
x-=fx/dfx;
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 Td=Tsd(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);
     System.out.println("ts="+ts+"Td="+Td+"Pc="+Pc+"p="+p+"vsl="+vsl+"vsg="+vsg+"x="+x);
     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 amonia system
if( Ps(t) > Pc)
  return "kritik?st? s?v?";
else
 {
 double ts=t;
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 if(v<vsl)
   return "s?v?";
 else if(v==vsl)
   return "doymu? s?v?";
 else if(v==vsg)
   return "doymu? buhar";
 else if(v>vsl && v<vsg)
   return "doymu? s?v? buhar kar???m?";
 else // (v>vsg)
   return "k?zg?n buhar";
 }
}

public String phase(double x)
{
//phase of the amonia system
if( x==3)
  return "kritik ?st? s?v?";
else
 {
 if(x<0)
   return "s?v?";
 else if( x==0)
   return "doymu? s?v?";
 else if(x==1)
   return "doymu? buhar";
 else if( x>0 && x<1)
   return "doymu? s?v? 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 xx=0;
double p=P(t,v);
double T=t+273.15;//K
double Psb=Psb(t);//buble 
double Psd=Psd(t);
double Pc=r.rd.Pc;

double vsl=vsl_t(t);
double vsg=vsg_t(t);

if(t>= r.Tc) xx=3;


else
{
  if(v==vsl)
    xx=0;
  else if(v==vsg)
    xx=1;
  else if(v>vsl && v<vsg)
    xx=(v-vsl)/(vsg-vsl);
  else if(v>vsg )
    xx=2;
  else if(v<vsl )
    xx=-1;  
//System.out.println("p="+p+"pc="+Pc+"psb="+Psb+"psd="+Psd+"v="+v+"vsl="+vsl+"vsg="+vsg+"x="+xx); 
}
   return xx;  
}
public double x_pv(double p, double v)
{
//phase and/or quality 
double xx=0;

double Ts=Ts(p);
double Pc=r.rd.Pc;
double vsl=vsl_t(Ts);
double vsg=vsg_t(Ts);
if(p >=Pc )
{ xx=3; }
else
{
 if(v<vsl) xx=-1;
 else if(v>vsg) xx=2;
 else if(v==vsl) xx=0;
 else if(v==vsg) xx=1;
 else {xx=(v-vsl)/(vsg-vsl);}
}
return xx;  
}

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[11];
char s1=s.charAt(0);
char s2=s.charAt(1);
r[6]=x_tv(v1,v2);
if(s.equals("tv")) //tv se?imi
{
double vsl=vsl_t(v1);
double vsg=vsg_t(v1);	
double qsl=vf(v1);
double qsg=vg(v1);
 r[0]=P(v1,v2);
 r[1]=v1; 
 r[2]=v2; 
 r[3]=h(v1,v2);
 r[4]=u(v1,v2);
 r[5]=s(v1,v2);	
 if(v1<this.r.Tc)
 {r[6]=(v2-qsl)/(qsg-qsl);}
 if(r[6]>=0 && r[6]<=1)
 { r[0]=Ps(v1);
   r[1]=v1; 
   r[2]=v2; 
   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];
 }
 } //tv se?imi sonu
else if(s.equals("pv"))  // pv se?imi
{
double ts=Ts(v1);
double vsl=vsl_t(ts);
double vsg=vsg_t(ts);
double qsl=vf(ts);
double qsg=vg(ts);
r[0]=v1;
r[1]=t_pv(v1,v2);
r[2]=v2;
r[3]=h(v1,v2);
r[4]=u(v1,v2);
r[5]=s(v1,v2);
if(v1<this.r.Pc)
{r[6]=(v2-qsl)/(qsg-qsl);}
else r[6]=3;
 if(r[6]>=0 && r[6]<=1)
 {r[0]=v1;
  r[1]=Ts(v1);
  r[2]=vsl*(1-r[6])+vsg*r[6];
  r[3]=h(r[1],vsl)*(1-r[6])+ h(r[1],vsg)*r[6];
  r[4]=u(r[1],vsl)*(1-r[6])+ u(r[1],vsg)*r[6];
  r[5]=s(r[1],vsl)*(1-r[6])+ s(r[1],vsg)*r[6];
 }
} //pv se?imi sonu
else if(s.equals("tp")) //tp se?imi
{
   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]);
} //tp se?imi sonu
else if(s.equals("th")) //th se?imi
{
 //double hsl=h(v1,vsl_t(v1));
 //double hsg=h(v1,vsg_t(v1));
 double vth=v_th(v1,v2);
 double hsl=hf(v1);
 double hsg=hg(v1);
   r[0]=P(v1,vth);
   r[1]=v1;
   r[2]=vth;
   r[3]=v2;
   r[4]=u(r[1],r[2]);
   r[5]=s(r[1],r[2]);
   if(v1<this.r.Tc)
   r[6]=(v2-hsl)/(hsg-hsl);
 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]=u(v1,vsl)*(1-r[6])+ u(v1,vsg)*r[6];
   r[5]=s(v1,vsl)*(1-r[6])+ s(v1,vsg)*r[6];
  }
} //th se?imi sonu
else if(s.equals("tu")) //tu se?imi
{
   double vtu=v_tu(v1,v2);
   double usl=uf(v1);
   double usg=ug(v1);
   r[0]=P(v1,vtu);
   r[1]=v1;
   r[2]=vtu;
   r[3]=h(r[1],r[2]);
   r[4]=v2;
   r[5]=s(r[1],r[2]);
   if(v1<this.r.Tc)
   r[6]=(v2-usl)/(usg-usl);
   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];
  }
} // tu se?imi sonu
else if(s.equals("ts")) //ts se?imi
{
   double ps=Ps(v1);
   double vts=v_ts(v1,v2);
   double ssl=sf(v1);
   double ssg=sg(v1);
   r[0]=P(v1,vts);
   r[1]=v1;
   r[2]=vts;
   r[3]=h(r[1],r[2]);
   r[4]=u(r[1],r[2]);
   r[5]=v2;
   if(v1<this.r.Tc) r[6]=(v2-ssl)/(ssg-ssl);
   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];
   r[5]=v2;
 }
} //ts se?imi sonu
else if(s.equals("tx"))
{
   r[6]=v2;
   r[0]=Ps(v1);
   r[1]=v1;
   double hl=hf(v1);
   double hg=hg(v1);
   double vsl=vsl_t(v1);
   double vsg=vsg_t(v1);
   double ul=hl-r[0]*vsl;
   double ug=hg-r[0]*vsg;
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=hl*(1-r[6])+ hg*r[6];
   r[4]=ul*(1-r[6])+ ug*r[6];
   r[5]=sf(v1)*(1-r[6])+ sg(v1)*r[6];
   //System.out.println("r5=s="+r[5]+"v1=t="+v1+"sf="+sf(v1)+"sg="+sg(v1));
} //tx se?imi sonu
else if(s.equals("vt"))
{
 r=property("tv",v2,v1);
} //vt se?imi sonu
else if(s.equals("vp"))
{
 r=property("pv",v2,v1);
} //vp se?imi sonu
else if(s.equals("pt"))
{
 r=property("tp",v2,v1);
} //pt se?imi sonu
else if(s.equals("ph"))
{
 double ts=Ts(v1);
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 double hsl=hf(ts);
 double hsg=hg(ts);
 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 < this.r.Pc)
   { r[6]=(v2-hsl)/(hsg-hsl);}
   else r[6]=3;
 if(r[6]>=0 && r[6]<=1)
 { r[0]=v1;
   r[1]=Ts(v1);
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=v2;
   r[4]=u(r[1],vsl)*(1-r[6])+ u(r[1],vsg)*r[6];
   r[5]=s(r[1],vsl)*(1-r[6])+ s(r[1],vsg)*r[6];
  }
} //ph se?imi sonu
else if(s.equals("pu"))
{
 double ts=Ts(v1);
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 double usl=uf(ts);
 double usg=ug(ts);
 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 < this.r.Pc)
   { r[6]=(v2-usl)/(usg-usl);}
   else r[6]=3;
 if(r[6]>=0 && r[6]<=1)
 { r[0]=v1;
   r[1]=Ts(v1);
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=h(r[1],vsl)*(1-r[6])+ h(r[1],vsg)*r[6];
   r[4]=v2;
   r[5]=s(r[1],vsl)*(1-r[6])+ s(r[1],vsg)*r[6];
  }
} //pu se?imi sonu 
else if(s.equals("ps"))
{
 double ts=Ts(v1);
 double vsl=vsl_t(ts);
 double vsg=vsg_t(ts);
 double ssl=sf(ts);
 double ssg=sg(ts);
 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 < this.r.Pc)
   { r[6]=(v2-ssl)/(ssg-ssl);}
   else r[6]=3;
 if(r[6]>=0 && r[6]<=1)
 { r[0]=v1;
   r[1]=Ts(v1);
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=h(r[1],vsl)*(1-r[6])+ h(r[1],vsg)*r[6];
   r[4]=u(r[1],vsl)*(1-r[6])+ u(r[1],vsg)*r[6];
   r[5]=v2;
  }
} //ps se?imi sonu
else if(s.equals("px"))
{
   r[6]=v2;
   double ts=Ts(v1);
   r[0]=v1;
   r[1]=Ts(v1);
   double hl=hf(ts);
   double hg=hg(ts);
   double sl=sf(ts);
   double sg=sg(ts);
   double vsl=vsl_t(ts);
   double vsg=vsg_t(ts);
   double ul=hl-r[0]*vsl;
   double ug=hg-r[0]*vsg;
   r[2]=vsl*(1-r[6])+ vsg*r[6];
   r[3]=hl*(1-r[6])+ hg*r[6];
   r[4]=ul*(1-r[6])+ ug*r[6];
   r[5]=sl*(1-r[6])+ sg*r[6];
} //px se?imi sonu
r[7]=1.0/r[2];
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, pressure        ";
   s1[1][0]="T, temperature     ";
   s1[2][0]="v, specific volume ";
   s1[3][0]="h, enthalpy        ";
   s1[4][0]="u, internal energy ";
   s1[5][0]="s, entropy         ";
   s1[6][0]="x, quality         ";
   s1[7][0]="   density         ";
   s1[8][0]="phase              ";
   if(unit.equals("SI"))
   {
   s1[0][2]=" kpa           ";
   s1[1][2]=" deg 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 vap/kg mix ";
   s1[7][2]=" kg/m^3        ";
   s1[8][2]=" ";
   }
   else
   {
   s1[0][2]=" lbf/in^2       ";
   s1[1][2]=" deg 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 vap/lbm mix";
   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]=" kpa           ";
   s2[1]=" deg C         ";
   s2[2]=" m^3/kg        ";
   s2[3]=" KJ/kg         ";
   s2[4]=" KJ/kg         ";
   s2[5]=" KJ/kg         ";
   s2[6]=" kg vap/kg mix ";
   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 vap/lbm mix";
   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]=" kpa           ";
   s2[1]=" deg C         ";
   s2[2]=" m^3/kg        ";
   s2[3]=" KJ/kg         ";
   s2[4]=" KJ/kg         ";
   s2[5]=" KJ/kg K       ";
   s2[6]=" kg vap/kg mix ";
   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 vap/lbm mix";
   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.25;
        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;
   }

}
