//======================================================
// java Termodinamik paketi
// Class sogut soğutma akışkanların özellikleri
// Dr. Turhan Coban
// Ege Üniversitesi Mühendislik fakultesi
// Makina Bölümü
// 
//======================================================
// Dosya adı   : sogutMBWR.java
// Bu dosya soğutucu akışkanların termodinamik özelliklerini 
// Değiştirilmiş Benedict-Webb-Rubin hal denklemiyle hesaplar
// 
// =====================================================

import java.io.*;
import java.net.*;
import java.util.*;

public class sogutMBWR
{
boolean SI=true;   //switch for SI/EN units
public sogutMBWRData rd;
// a reference state that is an ideal gas state for all refrigerants is taken
public double Tr; //K
public double vr;
final String u[]={"SI","EN"};
String unit=u[0];

public sogutMBWR(String rName)
{
  readref(rName);
  Tr=rd.Tr;
  vr=vsl_t(rd.Tr);
}

public String readGasNames()
{ String s1=rd.refrigerantName;String s2=readref("null");readref(s1);return s2;}
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/=145.04;
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]*=145.04;
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)
{
//---- String s aldığı değerler --------
//    tv tp th tu ts tx
// pv    pt ph pu ps px
// vp    vt

//---- Çıktı değerleri ---------
// r[0] P Basınç kPa
// r[1] t Sıcaklık derece C
// r[2] v Özgül hacim m^3/kg
// r[3] h KJ/kg 
// r[4] u KJ/kg
// r[5] s KJ/kgK
// r[6] x kg buhar/kg toplam
// r[7] ro yoğunluk kg/m^3

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")) //tv case
{
 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;
} //tv sonu
else if(s.equals("pv"))  // pv seçimi
{
 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;
} //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);
   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];
  }
} //th seçimi sonu
else if(s.equals("tu")) //tu seçimi
{
   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];
  }
} // tu seçimi sonu
else if(s.equals("ts")) //ts seçimi
{
   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];
 }

} //ts seçimi sonu
else if(s.equals("tx"))
{
   r[6]=v2;
   r[0]=Ps(v1);
   r[1]=v1;
   double hl=hsl_t(v1);
   double hg=hsg_t(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]=ssl_t(v1)*(1-r[6])+ ssg_t(v1)*r[6];
} //tx seçimi sonu
else if(s.equals("vt"))
{
 r=property_SI("tv",v2,v1);
} //vt seçimi sonu
else if(s.equals("vp"))
{
 r=property_SI("pv",v2,v1);
} //vp seçimi sonu
else if(s.equals("pt"))
{
 r=property_SI("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=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 < rd.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];
  }
} //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=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 < rd.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];
  }
} //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=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 < rd.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;
  }
} //ps seçimi sonu
else if(s.equals("px"))
{
   r[6]=v2;
   double ts=Ts(v1);
   r[0]=v1;
   r[1]=ts;
   double hl=hsl_t(ts);
   double hg=hsg_t(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]=ssl_t(v1)*(1-r[6])+ ssg_t(v1)*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 double P(double t, double v)
{
// Değiştirilmiş Benedict Webb Rubin Denklemi
// v m^3/kg özgül hacim
// V  m^3/kmol veya litre/mol özgül hacim
// t  Sıcaklık derece C
// T  Sıcaklık derece K
// P Basınç kPa
double V=v*rd.M;
double Vc=rd.vc*rd.M;
double T=t+273.15;
double T2=T*T;
double T3=T2*T;
double T4=T2*T2;
double pi=0;
double a[];
a=new double[15];
a[0]=rd.R*T;
a[1]=rd.B[0]*T + rd.B[1]*Math.sqrt(T) + rd.B[2] + rd.B[3]/T + rd.B[4]/T2;
a[2]=rd.B[5]*T + rd.B[6] + rd.B[7]/T + rd.B[8]/T2;
a[3]=rd.B[9]*T + rd.B[10] + rd.B[11]/T;
a[4]=rd.B[12];
a[5]=rd.B[13]/T + rd.B[14]/T2;
a[6]=rd.B[15]/T;
a[7]=rd.B[16]/T + rd.B[17]/T2;
a[8]=rd.B[18]/T2;
a[9]=rd.B[19]/T2 + rd.B[20]/T3;
a[10]=rd.B[21]/T2 + rd.B[22]/T4;
a[11]=rd.B[23]/T2 + rd.B[24]/T3;
a[12]=rd.B[25]/T2 + rd.B[26]/T4;
a[13]=rd.B[27]/T2 + rd.B[28]/T3;
a[14]=rd.B[29]/T2 + rd.B[30]/T3 + rd.B[31]/T4;
double vn[];
vn=new double[15];
double v1=V;
for(int i=1;i<14;i++){vn[i]=v1;v1*=V;}
double C=Math.exp(-rd.vc*rd.vc/v/v);
pi=a[0]/V+a[1]/vn[2]+a[2]/vn[3]+a[3]/vn[4]+a[4]/vn[5]+a[5]/vn[6]+a[6]/vn[7]+a[7]/vn[8]+a[8]/vn[9];
pi+=C*(a[9]/vn[3]+a[10]/vn[5]+a[11]/vn[7]+a[12]/vn[9]+a[13]/vn[11]+a[14]/vn[13]);
return pi*rd.mf;
}

double Pi(double t, double v)
{
double vsl=vsl_t(t);
double vsg=vsg_t(t);
if(v < vsl || v > vsg) {return P(t,v);}
else return Ps(t);
}


  public double dPdT(double t, double v)
  {
  double pi=0;
  double V=v*rd.M;
  double T=t+273.15;
  double T2=T*T;double T3=T2*T;double T4=T2*T2;double T5=T4*T;
  double vc2=rd.vc*rd.vc;
  double da[];
  da=new double[15];
  da[0]=rd.R;
  da[1]=rd.B[0] +0.5* rd.B[1]/Math.sqrt(T) - rd.B[3]/T2 -2*rd.B[4]/T3;
  da[2]=rd.B[5] - rd.B[7]/T2 -2* rd.B[8]/T3;
  da[3]=rd.B[9] - rd.B[11]/T2;
  da[4]=0;
  da[5]=-rd.B[13]/T2 - 2.0*rd.B[14]/T3;
  da[6]=-rd.B[15]/T2;
  da[7]=-rd.B[16]/T2 -2.0*rd.B[17]/T3;
  da[8]=-2.0*rd.B[18]/T3;
  da[9]=-2.0*rd.B[19]/T3 - 3.0* rd.B[20]/T4;
  da[10]=-2.0*rd.B[21]/T3 - 4.0*rd.B[22]/T5;
  da[11]=-2.0*rd.B[23]/T3 -3.0* rd.B[24]/T4;
  da[12]=-2.0*rd.B[25]/T3 -4.0* rd.B[26]/T5;
  da[13]=-2.0*rd.B[27]/T3 -3.0* rd.B[28]/T4;
  da[14]=-2.0*rd.B[29]/T3 -3.0* rd.B[30]/T4 -4.0* rd.B[31]/T5;
  double vn[];
  vn=new double[15];
  double v1=V;
  for(int i=1;i<14;i++){vn[i]=v1;v1*=V;}
  double C=Math.exp(-rd.vc*rd.vc/v/v);
  pi=da[0]/V+da[1]/vn[2]+da[2]/vn[3]+da[3]/vn[4]+da[4]/vn[5]+da[5]/vn[6]+da[6]/vn[7]+
  da[7]/vn[8]+da[8]/vn[9];
  pi+=C*(da[9]/vn[3]+da[10]/vn[5]+da[11]/vn[7]+da[12]/vn[9]+
  da[13]/vn[11]+da[14]/vn[13]);
  return pi*rd.mf;
  }

  public double t1(double Pin,double v)
  {
  // Newton - Raphson kök bulma metodu
  // ikinci derece türev hesabını da yapar
  // fonksiyon ve birinci türevi gereklidir
  double x=Ts(Pin);
  int nmax=500;
  double dx=1e-5;
  double x1m;
  double tolerance=1.0e-15;
  double fx,fxm,dfx,dfxm,d2fx;
  for(int i=0;i<nmax;i++)
  {
  fx=P(x,v)-Pin;
  fxm=P((x-dx),v)-Pin;
  dfx=dPdT(x,v);
  dfxm=dPdT((x-dx),v);
  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;
  }

// ======= v relations==========

public double v_tp(double t, double Pin)
{
//avaluate whole range to start root finding
// v m^3/kg
// V  m^3/kmol or liters/mol
// t temperature degree C
// T temperature degree K
// P pressure kPa
double maxit=50;
double xl=0.0006;
//double xu=1.0;
double xu=3.0;
double dx=(xu-xl)/maxit;
double min=99e50;
 double fx;
 double x;
 double xmin=xl;
 fx=P(t,xl)-Pin;
 for(x=xl;x<xu;x+=dx)
 {
 if(Math.abs(fx)<min) {min=Math.abs(fx);xmin=x;}
 fx=P(t,x)-Pin;
 } //end of while
// volume as a function of temperature and pressure
 double f0,df0;
 int nmax=150;
//you can adjust tolerance if you require better accuracy
 int i=0;
 //use result of search to start newton iteration
 x=xmin;
 dx=1.0e-5;
 double dx1=2.1144e-5;
 double x1m;
 double tolerance=1.0e-15;
 double fxm,fxm1,dfx,dfxm,d2fx;
 for(i=0;i<nmax;i++)
 {
    fx=P(t,x)-Pin;
    fxm=P(t,(x-dx))-Pin;
    fxm1=P(t,(x-dx1))-Pin;
    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_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=(T*8.3145/Pi)/rd.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_braket(double t,double Pi,double x1,double x2)
{
	//returns a bracketed root 
	//if a root is not bracketed, enlarge region outwards each time
	int NTRY=50;
	double a[]=new double[2];
	double FACTOR=1.343253;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Bad initial range in ZBRAC");
	f1=P(t,x1)-Pi;
	f2=P(t,x2)-Pi;
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;}
		if (Math.abs(f1) < Math.abs(f2))
			f1=P(t,(x1 += FACTOR*(x1-x2)))-Pi;
		else
			f2=P(t,(x2 += FACTOR*(x2-x1)))-Pi;
	}
	{a[0]=x1;a[1]=x2;return a;}
}

public double v_tp2(double t,double Pi)
{
//specific volume as a function of temperature and pressure	
	double T=t+273.15;
	double x=(T*rd.R/Pi);
	double dx=1.5436543624e-3;
	double x1=x-dx;
	double x2=x+dx;
	double aa[]=v_tp_braket(t,Pi,x1,x2);
	x1=aa[0];
	x2=aa[1];
	double ITMAX=100;
    double EPS=3.0e-8;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=P(t,a)-Pi,fb=P(t,b)-Pi,fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) System.out.println("xxx Root must be bracketed in ZBRENT");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=P(t,b)-Pi;
	}
	System.out.println("Maximum number of iterations exceeded in BRENT");
return b;
}

public double[] v_th_braket(double t,double hi,double x1,double x2)
{
	//returns a bracketed root 
	//if a root is not bracketed, enlarge region outwards each time
	int NTRY=50;
	double a[]=new double[2];
	double FACTOR=1.343253;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Bad initial range in ZBRAC");
	f1=h(t,x1)-hi;
	f2=h(t,x2)-hi;
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;}
		if (Math.abs(f1) < Math.abs(f2))
			f1=h(t,(x1 += FACTOR*(x1-x2)))-hi;
		else
			f2=h(t,(x2 += FACTOR*(x2-x1)))-hi;
	}
	{a[0]=x1;a[1]=x2;return a;}
}

public double v_th(double t,double hi)
{
//specific volume as a function of temperature and enthalpy	
	double T=t+273.15;
	double x=vg_app(t);
	double dx=1.5436543624e-3;
	double x1=x-dx;
	double x2=x+dx;
	//find a bracket with the root
	double aa[]=v_th_braket(t,hi,x1,x2);
	x1=aa[0];
	x2=aa[1];
	double ITMAX=100;
    double EPS=3.0e-8;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=h(t,a)-hi,fb=h(t,b)-hi,fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) System.out.println("Root must be bracketed in ZBRENT");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=h(t,b)-hi;
	}
	System.out.println("Maximum number of iterations exceeded in BRENT");
return b;
}

public double[] v_tu_braket(double t,double ui,double x1,double x2)
{
	//returns a bracketed root 
	//if a root is not bracketed, enlarge region outwards each time
	int NTRY=50;
	double a[]=new double[2];
	double FACTOR=1.343253;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Bad initial range in ZBRAC");
	f1=u(t,x1)-ui;
	f2=u(t,x2)-ui;
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;}
		if (Math.abs(f1) < Math.abs(f2))
			f1=u(t,(x1 += FACTOR*(x1-x2)))-ui;
		else
			f2=u(t,(x2 += FACTOR*(x2-x1)))-ui;
	}
	{a[0]=x1;a[1]=x2;return a;}
}


public double v_tu(double t,double ui)
{
//specific volume as a function of temperature and internal energy	
	double T=t+273.15;
	double x=vg_app(t);
	double dx=1.5436543624e-3;
	double x1=x-dx;
	double x2=x+dx;
	//find a bracket with the root
	double aa[]=v_tu_braket(t,ui,x1,x2);
	x1=aa[0];
	x2=aa[1];
	double ITMAX=100;
    double EPS=3.0e-8;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=u(t,a)-ui,fb=u(t,b)-ui,fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) System.out.println("Root must be bracketed in ZBRENT");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=u(t,b)-ui;
	}
	System.out.println("Maximum number of iterations exceeded in BRENT");
return b;
}
public double[] v_ts_braket(double t,double si,double x1,double x2)
{
	//returns a bracketed root 
	//if a root is not bracketed, enlarge region outwards each time
	int NTRY=50;
	double a[]=new double[2];
	double FACTOR=1.343253;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Bad initial range in ZBRAC");
	f1=s(t,x1)-si;
	f2=s(t,x2)-si;
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;}
		if (Math.abs(f1) < Math.abs(f2))
			f1=s(t,(x1 += FACTOR*(x1-x2)))-si;
		else
			f2=s(t,(x2 += FACTOR*(x2-x1)))-si;
	}
	{a[0]=x1;a[1]=x2;return a;}
}


public double v_ts(double t,double si)
{
//specific volume as a function of temperature and entropy	
	double T=t+273.15;
	double x=vg_app(t)+0.1;
	double dx=1.5436543624e-3;
	double x1=x-dx;
	double x2=x+dx;
	//find a bracket with the root
	double aa[]=v_ts_braket(t,si,x1,x2);
	x1=aa[0];
	x2=aa[1];
	double ITMAX=100;
    double EPS=3.0e-8;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=s(t,a)-si,fb=s(t,b)-si,fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) System.out.println("Root must be bracketed");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=s(t,b)-si;
	}
	System.out.println("Maximum number of iterations exceeded in BRENT");
return b;
}

public double v_tx(double t,double xi)
{
//specific volume as a function of temperature and quality	
return (1-xi)*vsl_t(t)+xi*vsg_t(t);
}

public double v_ph(double pi,double hi)
{
//specific volume as a function of pressure and enthalpy	
return v_tp(t_ph(pi,hi),pi);
}

public double v_pu(double pi,double ui)
{
//specific volume as a function of pressure and enthalpy	
return v_tp(t_pu(pi,ui),pi);
}
public double v_ps(double pi,double si)
{
//specific volume as a function of pressure and enthalpy	
return v_tp(t_ps(pi,si),pi);
}

public double v1(double t, double Pin)
{
double maxit=50;
double xl=0.0004;
//double xu=1.0;
double xu=2.0;
double dx=(xu-xl)/maxit;
double min=99e50;
double test;
double xr=0;
double es,ea;
double fxl,fxr,fxu;
int iter=0;
int sign1,sign2;
es=0.000001;
ea=1.1*es;
fxu=P(t,xu)-Pin;
fxl=P(t,xl)-Pin;
fxr=1.0;
//false position method iteration
while((fxr>es)&&(iter<maxit))
{
iter++;
if((xl+xu)!=0)
xr=xu-fxu*(xl-xu)/(fxl-fxu);
fxr= P(t,xr)-Pin;
sign1=(int)(fxl/Math.abs(fxl));
sign2=(int)(fxr/Math.abs(fxr));
if(sign1==sign2) {xl=xr;fxl=fxr;}
else             {xu=xr;fxu=fxr;}
test= fxl*fxr;
if(test==0.0)       ea=0;
else if(test<0.0)  xu=xr;
else               xl=xr;
} //end of while
return xr;
}

public double v2(double t, double Pin)
{	
//starting to search from big value through the small one
//avaluate whole range to start root finding
// v m^3/kg
// V  m^3/kmol or liters/mol
// t temperature degree C
// T temperature degree K
// P pressure kPa

double maxit=50;
double xl=0.0004;
//double xu=1.0;
double xu=2.0;
double dx=(xu-xl)/maxit;
double min=99e50;
double fx;
double x;
double xmin=xl;
   fx=P(t,xl)-Pin;
   for(x=xu;x>xl;x-=dx)
   {
   if(Math.abs(fx)<min) {min=Math.abs(fx);xmin=x;}
   fx=P(t,x)-Pin;
   //System.out.println("fx = "+fx+"x="+x);
   } //end of while
// volume as a function of temperature and pressure
   double f0,df0;
   int nmax=150;
//you can adjust tolerance if you require better accuracy
   int i=0;
   //use result of search to start newton iteration
   x=xmin;
   dx=1.0e-5;
   double dx1=2.1144e-5;
   double x1m;
   double tolerance=1.0e-15;
   double fxm,fxm1,dfx,dfxm,d2fx;
   for(i=0;i<nmax;i++)
   {
      fx=P(t,x)-Pin;
       fxm=P(t,(x-dx))-Pin;
      fxm1=P(t,(x-dx1))-Pin;
      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_pv(double Pin, double v)
 {
 // v m^3/kg
 // V  m^3/kmol or liters/mol
 // t temperature degree C
 // T temperature degree K
 // P pressure kPa
 // Temperature as a function of pressure and volume
  double f0,df0;
  int nmax=200;
//you can adjust tolerance if you require better accuracy
  int i=0;
  double x=Ts(Pin);
  //System.out.println("x="+x);
  double dx=1.0e-5;
  double dx1=2.4653124653e-5;
  double x1m;
  double tolerance=1.0e-10;
  double fx,fxm,fxm1,dfx,dfxm,d2fx;
  for(i=0;i<nmax;i++)
  {
     fx=P(x,v)-Pin;
     if(fx==0) return x; 
     fxm=P((x-dx),v)-Pin;
     fxm1=P((x-dx1),v)-Pin;
     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_braket(double pi,double ui,double x1,double x2)
{
	//returns a bracketed root 
	//if a root is not bracketed, enlarge region outwards each time
	int NTRY=50;
	double a[]=new double[2];
	double FACTOR=1.343253;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Bad initial range in ZBRAC");
	f1=u(x1,v_tp(x1,pi))-ui;
	f2=u(x2,v_tp(x2,pi))-ui;
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;}
		if (Math.abs(f1) < Math.abs(f2))
			f1=u((x1 += FACTOR*(x1-x2)),v_tp(x1,pi))-ui;
		else
			f2=u((x2 += FACTOR*(x2-x1)),v_tp(x2,pi))-ui;
	}
	{a[0]=x1;a[1]=x2;return a;}
}

public double t_pu(double pi,double ui)
{
//specific volume as a function of pressure and enthalpy	
	double x=Ts(pi);
	double dx=2.3e-1;
	double x1=x-dx;
	double x2=x+dx;
	//find a bracket with the root
	double aa[]=t_pu_braket(pi,ui,x1,x2);
	x1=aa[0];
	x2=aa[1];
	double ITMAX=100;
    double EPS=3.0e-8;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=u(a,v_tp(a,pi))-ui,fb=u(b,v_tp(b,pi))-ui,fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) System.out.println("Root must be bracketed");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=u(b,v_tp(b,pi))-ui;
	}
	System.out.println("Maximum number of iterations exceeded in BRENT");
return b;
}

public double[] t_ps_braket(double pi,double si,double x1,double x2)
{
	//returns a bracketed root 
	//if a root is not bracketed, enlarge region outwards each time
	int NTRY=50;
	double a[]=new double[2];
	double FACTOR=1.343253;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Bad initial range in ZBRAC");
	f1=s(x1,v_tp(x1,pi))-si;
	f2=s(x2,v_tp(x2,pi))-si;
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;}
		if (Math.abs(f1) < Math.abs(f2))
			f1=s((x1 += FACTOR*(x1-x2)),v_tp(x1,pi))-si;
		else
			f2=s((x2 += FACTOR*(x2-x1)),v_tp(x2,pi))-si;
	}
	{a[0]=x1;a[1]=x2;return a;}
}
public double[] t_ph_braket(double pi,double hi,double x1,double x2)
{
	//returns a bracketed root 
	//if a root is not bracketed, enlarge region outwards each time
	int NTRY=50;
	double a[]=new double[2];
	double FACTOR=1.343253;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Bad initial range in ZBRAC");
	f1=h(x1,v_tp(x1,pi))-hi;
	f2=h(x2,v_tp(x2,pi))-hi;
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;}
		if (Math.abs(f1) < Math.abs(f2))
			f1=h((x1 += FACTOR*(x1-x2)),v_tp(x1,pi))-hi;
		else
			f2=h((x2 += FACTOR*(x2-x1)),v_tp(x2,pi))-hi;
	}
	{a[0]=x1;a[1]=x2;return a;}
}
public double t_ph(double pi,double hi)
{
//specific volume as a function of pressure and enthalpy	
	double x=Ts(pi);
	double dx=2.3e-1;
	double x1=x-dx;
	double x2=x+dx;
	//find a bracket with the root
	double aa[]=t_ph_braket(pi,hi,x1,x2);
	x1=aa[0];
	x2=aa[1];
	double ITMAX=100;
    double EPS=3.0e-8;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=h(a,v_tp(a,pi))-hi,fb=h(b,v_tp(b,pi))-hi,fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) System.out.println("Root must be bracketed");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=h(b,v_tp(b,pi))-hi;
	}
	System.out.println("Maximum number of iterations exceeded in BRENT");
return b;
}
public double t_ps(double pi,double si)
{
//specific volume as a function of pressure and enthalpy	
	double x=Ts(pi);
	double dx=2.3e-1;
	double x1=x-dx;
	double x2=x+dx;
	//find a bracket with the root
	double aa[]=t_ps_braket(pi,si,x1,x2);
	x1=aa[0];
	x2=aa[1];
	double ITMAX=100;
    double EPS=3.0e-8;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=s(a,v_tp(a,pi))-si,fb=s(b,v_tp(b,pi))-si,fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) System.out.println("Root must be bracketed");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=s(b,v_tp(b,pi))-si;
	}
	System.out.println("Maximum number of iterations exceeded in BRENT");
return b;
}

public double vg_app(double t)
{
return Math.exp(rd.vg[0]+rd.vg[1]*t+rd.vg[2]*t*t+rd.vg[3]*t*t*t+
   rd.vg[4]*t*t*t*t+rd.vg[5]*t*t*t*t*t);
}

public double vsg_t(double t)
{
// saturation vapor volume m^3/kg
  if(t>=rd.Tc) return rd.vc;
//else   
// volume as a function of temperature and pressure
   double fx,f0,df0;
   double Pin=Ps(t);
   int nmax=150;
   double x1;
//you can adjust tolerance if you require better accuracy
   int i=0;
   double x=vg_app(t);
   x1=x;
   double dx=1.0e-6;
   double dx1=2.013e-6;
   double x1m;
   double tolerance=1.0e-8;
   double fxm,fxm1,dfx,dfxm,d2fx;
   for(i=0;i<nmax;i++)
   {
      fx=P(t,x)-Pin;
      fxm=P(t,(x-dx))-Pin;
      fxm1=P(t,(x-dx1))-Pin;
      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 x1;
}


public double hfg(double t)
{
double x;
if(t>=rd.Tc) x=0;
else x=(t+273.15)*(vsg_t(t)-vsl_t(t))*dPsdT(t);
return x;
}

public double sfg(double t)
{
double x;
if(t>=rd.Tc) x=0;
else x=(vsg_t(t)-vsl_t(t))*dPsdT(t);
return x;	
}

public double rogs(double t)
{ 
//density of saturated vapor	
return 1.0/vsg_t(t);
}
	
public double rols(double t)
{
// saturation liquid density as a function of saturation temperature
double T=t+273.15;
double log10=Math.log(10);
double Tr=T/(rd.Tc+273.15);
if(Tr > 1) return rd.dc;
double df;
double t1=Math.pow((1.0-Tr),(1.0/3.0));
double t2=Math.pow((1.0-Tr),(2.0/3.0));
double t4=Math.pow((1.0-Tr),(4.0/3.0));
df=rd.adf[0]+rd.adf[1]*t1+rd.adf[2]*t2+rd.adf[3]*(1.0-Tr)+rd.adf[4]*t4;
return df;
}

public double vsl_t(double t)
{
//saturation liquid volume
	if(t >= rd.Tc) return rd.vc;
    return 1.0/rols(t);
}

  public double dPdT1(double t,double v)
  {
// This method calculates derivatives of a  function
// accuracy of method can be adjusted by changing
// variables h0 and n
  double x=t;
  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]=( P((x + h[i]),v) - P(( x - h[i]),v))/(2.0*h[i]);
  }
  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;
  }


public double dPdV(double t, double v)
{
 double pi=0;
 double V=v*rd.M;
 double T=t+273.15;
 double T2=T*T;double T3=T2*T;double T4=T2*T2;
 double Vc=rd.vc*rd.M;
 double Vc2=Vc*Vc;
 double a[];
 a=new double[15];
 a[0]=rd.R*T;
 a[1]=rd.B[0]*T + rd.B[1]*Math.pow(T,0.5) + rd.B[2] + rd.B[3]/T + rd.B[4]/T2;
 a[2]=rd.B[5]*T + rd.B[6] + rd.B[7]/T + rd.B[8]/T2;
 a[3]=rd.B[9]*T + rd.B[10] + rd.B[11]/T;
 a[4]=rd.B[12];
 a[5]=rd.B[13]/T + rd.B[14]/T2;
 a[6]=rd.B[15]/T;
 a[7]=rd.B[16]/T + rd.B[17]/T2;
 a[8]=rd.B[18]/T2;
 a[9]=rd.B[19]/T2 + rd.B[20]/T3;
 a[10]=rd.B[21]/T2 + rd.B[22]/T4;
 a[11]=rd.B[23]/T2 + rd.B[24]/T3;
 a[12]=rd.B[25]/T2 + rd.B[26]/T4;
 a[13]=rd.B[27]/T2 + rd.B[28]/T3;
 a[14]=rd.B[29]/T2 + rd.B[30]/T3 + rd.B[31]/T4;
 double vn[];
 vn=new double[15];
 double v1=V;
 for(int i=1;i<15;i++){vn[i]=v1;v1*=V;}
 double C=Math.exp(-rd.vc*rd.vc/v/v);
 pi=-(a[0]/vn[2]+2.0*a[1]/vn[3]+3.0*a[2]/vn[4]+4.0*a[3]/vn[5]+5.0*a[4]/vn[6]+
     6.0*a[5]/vn[7]+7.0*a[6]/vn[8]+8.0*a[7]/vn[9]+9.0*a[8]/vn[10]);
 pi+=2.0*Vc2/vn[3]*C*(a[9]/vn[3]+a[10]/vn[5]+a[11]/vn[7]+a[12]/vn[9]+
     a[13]/vn[11]+a[14]/vn[13]);
 pi-=C*(3.0*a[9]/vn[4]+5.0*a[10]/vn[6]+7.0*a[11]/vn[8]+9.0*a[12]/vn[10]+
        11.0*a[13]/vn[12]+13.0*a[14]/vn[14]);
 return pi*rd.mf;
}

public double dPdV1(double t,double v)
      {
// This method calculates derivatives of a  function
// accuracy of method can be adjusted by changing
// variables h0 and n
      double x=v;
      double h0=0.001*v;
      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]=( P(t,(x + h[i])) - P(t,( x - h[i])))/(2.0*h[i]);
      }
      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;
}

    public  double dA1(double t,double v)
   {
   //numerical integration of A just to check value
   double T=t+273.15;
    double a=v;
    double b=50.0;
    double s[],w[];
    int i;
    s=new double[30];
    w=new double[30];
    s[ 0] = .15532579626752470000E-02;
    s[ 1] = .81659383601264120000E-02;
    s[ 2] = .19989067515846230000E-01;
    s[ 3] = .36899976285362850000E-01;
    s[ 4] = .58719732103973630000E-01;
    s[ 5] = .85217118808615820000E-01;
    s[ 6] = .11611128394758690000E+00;
    s[ 7] = .15107475260334210000E+00;
    s[ 8] = .18973690850537860000E+00;
    s[ 9] = .23168792592899010000E+00;
    s[10] = .27648311523095540000E+00;
    s[11] = .32364763723456090000E+00;
    s[12] = .37268153691605510000E+00;
    s[13] = .42306504319570830000E+00;
    s[14] = .47426407872234120000E+00;
    s[15] = .52573592127765890000E+00;
    s[16] = .57693495680429170000E+00;
    s[17] = .62731846308394490000E+00;
    s[18] = .67635236276543910000E+00;
    s[19] = .72351688476904450000E+00;
    s[20] = .76831207407100990000E+00;
    s[21] = .81026309149462140000E+00;
    s[22] = .84892524739665800000E+00;
    s[23] = .88388871605241310000E+00;
    s[24] = .91478288119138420000E+00;
    s[25] = .94128026789602640000E+00;
    s[26] = .96310002371463720000E+00;
    s[27] = .98001093248415370000E+00;
    s[28] = .99183406163987350000E+00;
    s[29] = .99844674203732480000E+00;
   w[ 0] = .39840962480827790000E-02;
   w[ 1] = .92332341555455000000E-02;
   w[ 2] = .14392353941661670000E-01;
   w[ 3] = .19399596284813530000E-01;
   w[ 4] = .24201336415292590000E-01;
   w[ 5] = .28746578108808720000E-01;
   w[ 6] = .32987114941090080000E-01;
   w[ 7] = .36877987368852570000E-01;
   w[ 8] = .40377947614710090000E-01;
   w[ 9] = .43449893600541500000E-01;
   w[10] = .46061261118893050000E-01;
   w[11] = .48184368587322120000E-01;
   w[12] = .49796710293397640000E-01;
   w[13] = .50881194874202750000E-01;
   w[14] = .51426326446779420000E-01;
   w[15] = .51426326446779420000E-01;
   w[16] = .50881194874202750000E-01;
   w[17] = .49796710293397640000E-01;
   w[18] = .48184368587322120000E-01;
   w[19] = .46061261118893050000E-01;
   w[20] = .43449893600541500000E-01;
   w[21] = .40377947614710090000E-01;
   w[22] = .36877987368852570000E-01;
   w[23] = .32987114941090080000E-01;
   w[24] = .28746578108808720000E-01;
   w[25] = .24201336415292590000E-01;
   w[26] = .19399596284813530000E-01;
   w[27] = .14392353941661670000E-01;
   w[28] = .92332341555455000000E-02;
   w[29] = .39840962480827790000E-02;
   int n=30;
   double z=0;
   double x,y;
   for(i=0;i<n;i++)
   {
   x=(b+a)/2.0+(b-a)/2.0*s[i];
   y=P(t,x)-rd.R*rd.mf*T/x;
   z+=(b-a)/2*w[i]*y;
   }
   for(i=0;i<n;i++)
   {
   x=(b+a)/2.0+(b-a)/2.0*(-s[i]);
   y=P(t,x)-rd.R*rd.mf*T/x;
   z+=(b-a)/2.0*w[i]*y;
   }
   return z;
   }

   public double dA(double t, double v)
   {
      double V0=vsg_t(0);
      double V=v*rd.M;
      double T=t+273.15;
      double a[];
      a=new double[15];
      a[0]=rd.R*T;
      double T2=T*T;double T3=T2*T;double T4=T3*T;
      a[1]=rd.B[0]*T + rd.B[1]*Math.sqrt(T) + rd.B[2] + rd.B[3]/T + rd.B[4]/T2;
      a[2]=rd.B[5]*T + rd.B[6] + rd.B[7]/T + rd.B[8]/T2;
      a[3]=rd.B[9]*T + rd.B[10] + rd.B[11]/T;
      a[4]=rd.B[12];
      a[5]=rd.B[13]/T + rd.B[14]/T2;
      a[6]=rd.B[15]/T;
      a[7]=rd.B[16]/T + rd.B[17]/T2;
      a[8]=rd.B[18]/T2;
      a[9]=rd.B[19]/T2 + rd.B[20]/T3;
      a[10]=rd.B[21]/T2 + rd.B[22]/T4;
      a[11]=rd.B[23]/T2 + rd.B[24]/T3;
      a[12]=rd.B[25]/T2 + rd.B[26]/T4;
      a[13]=rd.B[27]/T2 + rd.B[28]/T3;
      a[14]=rd.B[29]/T2 + rd.B[30]/T3 + rd.B[31]/T4;
      double V2=V*V;double V3=V2*V;double V4=V2*V2;double V5=V2*V3;double V6=V3*V3;double V7=V4*V3;
      double V8=V6*V2;double V9=V*V8;double V10=V2*V8;
      double Vc=rd.vc*rd.M;
      double Vc2=Vc*Vc;double Vc3=Vc2*Vc;double Vc4=Vc2*Vc2;double Vc5=Vc2*Vc3;double Vc6=Vc3*Vc3;
      double Vc7=Vc4*Vc3;double Vc8=Vc6*Vc2;double Vc9=Vc*Vc8;double Vc10=Vc2*Vc8;double Vc11=Vc3*Vc8;double Vc12=Vc2*Vc10;
      double C=Math.exp(-Vc2/V2);
      double dAi;
      dAi=0.5*C/V10/Vc12*(
      a[9]*V10*Vc10+
      a[10]*(V8*Vc10+V10*Vc8)+
      a[11]*(V6*Vc10+2.0*V8*Vc8+2.0*V10*Vc6)+
      a[12]*(V4*Vc10+3.0*V6*Vc8+6.0*V8*Vc6+6.0*V10*Vc4)+
      a[13]*(V2*Vc10+4.0*V4*Vc8+12.0*V6*Vc6+24.0*V8*Vc4+24.0*V10*Vc2)+
      a[14]*(Vc10+5*V2*Vc8+20.0*V4*Vc6+60.0*V6*Vc4+120.0*V8*Vc2+120.0*V10))-
      (a[1]/V+0.5*a[2]/V2+1.0/3.0*a[3]/V3+1.0/4.0*a[4]/V4+1.0/5.0*a[5]/V5+1.0/6.0*a[6]/V6+1.0/7.0*a[7]/V7+1.0/8.0*a[8]/V8);
      double dAinf=0.5/Vc12*(a[9]*Vc10+a[10]*Vc8+2.0*a[11]*Vc6+6.0*a[12]*Vc4+24.0*a[13]*Vc2+120*a[14]);
      return rd.mf*((dAinf-dAi)/rd.M)-rd.R*T*Math.log(V/V0);  //KJ/kg
      }

double dG(double t,double v)
{
  //gibs free energy
double T=t+273.15;
double V=v*rd.M;
return dA(t,v)-(P(t,v)*V-rd.R*T);
}

public double dAdT(double t, double v)
{
        double V0=vsg_t(0);
        double V=v*rd.M;
        double T=t+273.15;
        double T2=T*T;double T3=T2*T;double T4=T2*T2;double T5=T4*T;
        double da[];
        da=new double[15];
        da[0]=rd.R;
        da[1]=rd.B[0] +0.5* rd.B[1]/Math.sqrt(T) - rd.B[3]/T2 -2*rd.B[4]/T3;
        da[2]=rd.B[5] - rd.B[7]/T2 -2* rd.B[8]/T3;
        da[3]=rd.B[9] - rd.B[11]/T2;
        da[4]=0;
        da[5]=-rd.B[13]/T2 - 2.0*rd.B[14]/T3;
        da[6]=-rd.B[15]/T2;
        da[7]=-rd.B[16]/T2 -2.0*rd.B[17]/T3;
        da[8]=-2.0*rd.B[18]/T3;
        da[9]=-2.0*rd.B[19]/T3 - 3.0* rd.B[20]/T4;
        da[10]=-2.0*rd.B[21]/T3 - 4.0*rd.B[22]/T5;
        da[11]=-2.0*rd.B[23]/T3 -3.0* rd.B[24]/T4;
        da[12]=-2.0*rd.B[25]/T3 -4.0* rd.B[26]/T5;
        da[13]=-2.0*rd.B[27]/T3 -3.0* rd.B[28]/T4;
        da[14]=-2.0*rd.B[29]/T3 -3.0* rd.B[30]/T4 -4.0* rd.B[31]/T5;
        double V2=V*V;double V3=V2*V;double V4=V2*V2;double V5=V2*V3;double V6=V3*V3;double V7=V4*V3;
        double V8=V6*V2;double V9=V*V8;double V10=V2*V8;
        double Vc=rd.vc*rd.M;
        double Vc2=Vc*Vc;double Vc3=Vc2*Vc;double Vc4=Vc2*Vc2;double Vc5=Vc2*Vc3;double Vc6=Vc3*Vc3;
        double Vc7=Vc4*Vc3;double Vc8=Vc6*Vc2;double Vc9=Vc*Vc8;double Vc10=Vc2*Vc8;double Vc11=Vc3*Vc8;double Vc12=Vc2*Vc10;
        double C=Math.exp(-rd.vc*rd.vc/v/v);
        double dAi;
        dAi=0.5*C/V10/Vc12*(
        da[9]*V10*Vc10+
        da[10]*(V8*Vc10+V10*Vc8)+
        da[11]*(V6*Vc10+2.0*V8*Vc8+2.0*V10*Vc6)+
        da[12]*(V4*Vc10+3.0*V6*Vc8+6.0*V8*Vc6+6.0*V10*Vc4)+
        da[13]*(V2*Vc10+4.0*V4*Vc8+12.0*V6*Vc6+24.0*V8*Vc4+24.0*V10*Vc2)+
        da[14]*(Vc10+5*V2*Vc8+20.0*V4*Vc6+60.0*V6*Vc4+120.0*V8*Vc2+120.0*V10))-
        (da[1]/V+0.5*da[2]/V2+1.0/3.0*da[3]/V3+1.0/4.0*da[4]/V4+1.0/5.0*da[5]/V5+1.0/6.0*da[6]/V6+1.0/7.0*da[7]/V7+1.0/8.0*da[8]/V8);
        double dAinf=0.5/Vc12*(da[9]*Vc10+da[10]*Vc8+2.0*da[11]*Vc6+6.0*da[12]*Vc4+24.0*da[13]*Vc2+120*da[14]);
        return rd.mf*((dAinf-dAi)/rd.M)-rd.R*Math.log(V/V0); //converted to KJ/kgK
}


public double dH(double t, double v)
{
//enthalpy departure function
double V=v*rd.M;
double T=t+273.15;
double da=dA(t,v);
double tds=-T*dAdT(t,v);
double pv=P(t,v)*v;
double rt=rd.R*rd.mf/rd.M*T;
double pvrt=pv-rt;
double tot=da+tds+pvrt;
return tot;
}


public double dS(double t, double v)
{
//entropy departure function
return -dAdT(t,v);
}


public double H0(double t)
{
//ideal gas enthalpy
        double T=t+273.15;
        double Tr=rd.Tr+273.15;
        //System.out.println("T="+T+"Tr="+Tr);
        double T2=T*T;double T3=T2*T;double T4=T2*T2;
        double Tr2=Tr*Tr;double Tr3=Tr2*Tr;double Tr4=Tr2*Tr2;
        double hh =(rd.acp[0]*(T-Tr)+rd.acp[1]/2.0*(T2-Tr2)+rd.acp[2]/3.0*(T3-Tr3))/rd.M;
        return hh; //KJ/kg K
}


public double h1(double t, double v)
{
double dh1=dH(t,v);
double tr=rd.Tr;
double vr=vsg_t(tr);
double dh2=dH(tr,vr);
double h0=H0(t);
double h0r=H0(tr);
double tot=rd.hr+hfg(tr)+h0+dh1-dh2;
//System.out.println("rd.hr="+rd.hr+"hfg="+hfg(Tr)+"dh1="+dh1+"dh2="+dh2+"h0="+h0+"h0r="+h0r+"tot"+tot);
return tot;
}

public double h(double t, double v)
{
double h1=h1(t,v);
double x=x_tv(t,v);
if(x>1)
{
 h1=h1(t,v);
}
else if(x<0)
{
 double P=P(t,v);
 h1=hsl_t(t)+vsl_t(t)*(P(t,v)-Ps(t));
}
else if(x==0)
{
 h1=hsl_t(t);
}
else if(x==1)
{
 h1=hsg_t(t);
}
else if(x>0 || x<1)
{	
 h1=hsl_t(t)*(1-x)+hsg_t(t)*x;
}
return h1;
}

public double u(double t, double v)
{
double h1=h1(t,v);
double v1=0;
double x=x_tv(t,v);
if(x>1)
{
 h1=h1(t,v)-P(t,v)*v;
}
else if(x<0)
{
 double P=P(t,v);
 h1=hsl_t(t)-vsl_t(t)*Ps(t);
}
else if(x==0)
{
 h1=hsl_t(t)-vsl_t(t)*Ps(t);
}
else if(x==1)
{
 h1=hsg_t(t)-vsg_t(t)*Ps(t);
}
else if(x>0 || x<1)
{	
 h1=hsl_t(t)*(1-x)+hsg_t(t)*x;
 v1=vsl_t(t)*(1-x)+vsg_t(t)*x;
 h1=h1-v1*Ps(t);
}
return h1;
}

public double s1(double t,double v)
{
double ds1=dS(t,v);
double tr=rd.Tr;
double vr=vsg_t(0);
double ds2=dS(0,vr);
double s0=S0_tv(t,v);
double s0r=S0_tv(tr,vr);
//double tot=rd.sr+sfg(0)+s0+ds1-ds2;
double tot=sfg(0)+s0+ds1-ds2;
//System.out.println("rd.sr="+rd.sr+"sfg="+sfg(tr)+"ds1="+ds1+"ds2="+ds2+"s0="+s0+"s0r="+s0r+"tot"+tot);
return tot;
}

public double s(double t, double v)
{
double s1=s1(t,v);
double x=x_tv(t,v);
if(x>1)
{
 s1=s1(t,v);
}
else if(x<0)
{
 double P=P(t,v);
 s1=ssl_t(t)+vsl_t(t)*(P(t,v)-Ps(t));
}
else if(x==0)
{
 s1=ssl_t(t);
}
else if(x==1)
{
 s1=ssg_t(t);
}
else if(x>0 || x<1)
{	
 s1=ssl_t(t)*(1-x)+ssg_t(t)*x;
}
return s1;
}


public double ssg_t(double t)
{
  return s1(t,vsg_t(t));
}

public double ssl_t(double t)
{
  return s1(t,vsg_t(t))-sfg(t);
}

public double hsg_t(double t)
{
  return h1(t,vsg_t(t));
}

public double hsl_t(double t)
{
  return h1(t,vsg_t(t))-hfg(t);
}

public double S0(double t)
{
//ideal gas entropy
double T=t+273.15;
double Tr=273.15;
double T2=T*T;double T3=T2*T;double T4=T2*T2;
double Tr2=Tr*Tr;double Tr3=Tr2*Tr;double Tr4=Tr2*Tr2;
double ss = rd.sr+(rd.acp[0]*Math.log(T/Tr)+rd.acp[1]*(T-Tr)+rd.acp[2]/2.0*(T2-Tr2))/rd.M;
return ss; //KJ/kmol K
}

public double S0_tv(double t,double v)
{
//ideal gas entropy
//double v0=vsg_t(t);
double T=t+273.15;
double Tr=273.15;
double T2=T*T;double T3=T2*T;double T4=T2*T2;
double Tr2=Tr*Tr;double Tr3=Tr2*Tr;double Tr4=Tr2*Tr2;
double ss = rd.sr+(rd.acp[0]*Math.log(T/Tr)+rd.acp[1]*(T-Tr)+rd.acp[2]/2.0*(T2-Tr2))/rd.M;
//System.out.println("ss="+ss);
return ss; //KJ/kmol K
}

public double Ps(double t)
      {
      // saturation pressure as a function of saturation temperature
      double MM=rd.aVP[0][9];
      double T=t+273.16;           
      //if( T >= rd.Tc) return rd.Pc; 
      double pi=0.0;
      if(t<=rd.th[0])
      {
      pi=rd.aVP[0][0]+rd.aVP[0][1]/T+rd.aVP[0][2]/(T*T)+rd.aVP[0][3]*Math.log10(T)+rd.aVP[0][4]*T
      +rd.aVP[0][5]*T*T+rd.aVP[0][6]*T*T*T;
      if(rd.aVP[0][8]!=0)
         pi+=(rd.aVP[0][7]*(rd.aVP[0][8]-T)/T)*Math.log10((rd.aVP[0][8]-T)*MM);
      }
      else if(t>rd.th[0] && t<= rd.th[1])
      {
      pi=rd.aVP[1][0]+rd.aVP[1][1]/T+rd.aVP[1][2]/(T*T)+rd.aVP[1][3]*Math.log10(T)+rd.aVP[1][4]*T+
      rd.aVP[1][5]*T*T+rd.aVP[1][6]*T*T*T;
         if(rd.aVP[1][8]!=0)
         pi+=(rd.aVP[1][7]*(rd.aVP[1][8]-T)/T)*Math.log10((rd.aVP[1][8]-T)*MM);
      }   
      else if(t>rd.th[1] && t<= rd.th[2])
     {
      pi=rd.aVP[2][0]+rd.aVP[2][1]/T+rd.aVP[2][2]/(T*T)+rd.aVP[2][3]*Math.log10(T)+rd.aVP[2][4]*T+
      rd.aVP[2][5]*T*T+rd.aVP[2][6]*T*T*T;
         if(rd.aVP[2][8]!=0)
         pi+=(rd.aVP[2][7]*(rd.aVP[2][8]-T)/T)*Math.log10((rd.aVP[2][8]-T)*MM);
      }   
      else if(t>rd.th[2] && t<= rd.th[3])
     {
      pi=rd.aVP[3][0]+rd.aVP[3][1]/T+rd.aVP[3][2]/(T*T)+rd.aVP[3][3]*Math.log10(T)+rd.aVP[3][4]*T+
      rd.aVP[3][5]*T*T+rd.aVP[3][6]*T*T*T;
         if(rd.aVP[3][8]!=0)
         pi+=(rd.aVP[3][7]*(rd.aVP[3][8]-T)/T)*Math.log10((rd.aVP[3][8]-T)*MM);
      }   
      else if(t>rd.th[3] && t<= rd.th[4])
     {
      pi=rd.aVP[4][0]+rd.aVP[4][1]/T+rd.aVP[4][2]/(T*T)+rd.aVP[4][3]*Math.log10(T)+rd.aVP[4][4]*T+
      rd.aVP[4][5]*T*T+rd.aVP[4][6]*T*T*T;
         if(rd.aVP[4][8]!=0)
         pi+=(rd.aVP[4][7]*(rd.aVP[4][8]-T)/T)*Math.log10((rd.aVP[4][8]-T)*MM);
      }   
      else if(t> rd.th[4])
     {
      pi=rd.aVP[5][0]+rd.aVP[5][1]/T+rd.aVP[5][2]/(T*T)+rd.aVP[5][3]*Math.log10(T)+rd.aVP[5][4]*T+
      rd.aVP[5][5]*T*T+rd.aVP[5][6]*T*T*T;
         if(rd.aVP[5][8]!=0)
         pi+=(rd.aVP[5][7]*(rd.aVP[5][8]-T)/T)*Math.log10((rd.aVP[5][8]-T)*MM);
      }   
      //System.out.println("pi="+pi);
      pi=Math.pow(10.0,pi);
      //System.out.println("pi="+pi);
      //if(rd.refrigerantName=="R123" || rd.refrigerantName=="R134a") pi*=1e-3; 
      return pi; //kPa
      }

public double dPsdT(double t)
{
//derivative of saturation pressure as a function of saturation temperature
double MM=rd.aVP[0][9];
double T=t+273.15;
double log10=Math.log(10);
double pi=0.0;
if(t<=rd.th[0])
{
pi=(-rd.aVP[0][1]/(T*T)-2.0*rd.aVP[0][2]/(T*T*T)+rd.aVP[0][4]+2.0*rd.aVP[0][5]*T+3.0*rd.aVP[0][6]*T*T)*Math.log(10.0)
+rd.aVP[0][3]/T;
if(rd.aVP[0][7]!=0.0)
{
pi+=rd.aVP[0][7]*(-1.0/T-(rd.aVP[0][8]-T)/(T*T))*Math.log((rd.aVP[0][8]-T)*MM)+
(rd.aVP[0][7]*(rd.aVP[0][8]-T)/T)*(-1.0/(rd.aVP[0][8]-T));
}
}
else if(t>rd.th[0] && t<= rd.th[1])
{
pi=(-rd.aVP[1][1]/(T*T)-2.0*rd.aVP[1][2]/(T*T*T)+rd.aVP[1][4]+2.0*rd.aVP[1][5]*T+3.0*rd.aVP[1][6]*T*T)*Math.log(10.0)
+rd.aVP[1][3]/T;
if(rd.aVP[1][7]!=0.0)
{
pi+=rd.aVP[1][7]*(-1.0/T-(rd.aVP[1][8]-T)/(T*T))*Math.log((rd.aVP[1][8]-T)*MM)+
(rd.aVP[1][7]*(rd.aVP[1][8]-T)/T)*(-1.0/(rd.aVP[1][8]-T));
}
}
else if(t>rd.th[1] && t<= rd.th[2])
{
pi=(-rd.aVP[2][1]/(T*T)-2.0*rd.aVP[2][2]/(T*T*T)+rd.aVP[2][4]+2.0*rd.aVP[2][5]*T+3.0*rd.aVP[2][6]*T*T)*Math.log(10.0)
+rd.aVP[2][3]/T;
if(rd.aVP[2][7]!=0.0)
{
pi+=rd.aVP[2][7]*(-1.0/T-(rd.aVP[2][8]-T)/(T*T))*Math.log((rd.aVP[2][8]-T)*MM)+
(rd.aVP[2][7]*(rd.aVP[2][8]-T)/T)*(-1.0/(rd.aVP[2][8]-T));
}
}
else if(t>rd.th[2] && t<= rd.th[3])
{
pi=(-rd.aVP[3][1]/(T*T)-2.0*rd.aVP[3][2]/(T*T*T)+rd.aVP[3][4]+2.0*rd.aVP[3][5]*T+3.0*rd.aVP[3][6]*T*T)*Math.log(10.0)
+rd.aVP[3][3]/T;
if(rd.aVP[3][7]!=0.0)
{
pi+=rd.aVP[3][7]*(-1.0/T-(rd.aVP[3][8]-T)/(T*T))*Math.log((rd.aVP[3][8]-T)*MM)+
(rd.aVP[3][7]*(rd.aVP[3][8]-T)/T)*(-1.0/(rd.aVP[3][8]-T));
}
}
else if(t>rd.th[3] && t<= rd.th[4])
{
pi=(-rd.aVP[4][1]/(T*T)-2.0*rd.aVP[4][2]/(T*T*T)+rd.aVP[4][4]+2.0*rd.aVP[4][5]*T+3.0*rd.aVP[4][6]*T*T)*Math.log(10.0)
+rd.aVP[4][3]/T;
if(rd.aVP[4][7]!=0.0)
{
pi+=rd.aVP[4][7]*(-1.0/T-(rd.aVP[4][8]-T)/(T*T))*Math.log((rd.aVP[4][8]-T)*MM)+
(rd.aVP[4][7]*(rd.aVP[4][8]-T)/T)*(-1.0/(rd.aVP[4][8]-T));
}
}
else if(t> rd.th[4])      
{
pi=(-rd.aVP[5][1]/(T*T)-2.0*rd.aVP[5][2]/(T*T*T)+rd.aVP[5][4]+2.0*rd.aVP[5][5]*T+3.0*rd.aVP[5][6]*T*T)*Math.log(10.0)
+rd.aVP[5][3]/T;
if(rd.aVP[5][7]!=0.0)
{
pi+=rd.aVP[5][7]*(-1.0/T-(rd.aVP[5][8]-T)/(T*T))*Math.log((rd.aVP[5][8]-T)*MM)+
(rd.aVP[5][7]*(rd.aVP[5][8]-T)/T)*(-1.0/(rd.aVP[5][8]-T));
}
}         
pi*=Ps(t);
return pi; //kPa
}

public double Ts(double ps)
{
double x;
if(ps>rd.Pc) x=rd.Tc;
else
{
// Saturation Temperature as a function of saturation pressure
double f0,df0;
int nmax=150;
//you can adjust tolerance if you require better accuracy
int i=0;
double p1=Math.log(ps);
double p2=p1*p1;
double p3=p2*p1;
double p4=p2*p2;
double p5=p4*p1;
x=rd.aa[0]+rd.aa[1]*p1+rd.aa[2]*p2+rd.aa[3]*p3+rd.aa[4]*p4+rd.aa[5]*p5+
rd.aa[6]*p5*p1+rd.aa[7]*p5*p2+rd.aa[8]*p4*p4;
//System.out.println("x="+x);
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;
   if(fx==0) return x;
   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 Cv0(double t)
{
//ideal gas constant volume heat capacity

return Cp0(t)-8.314471/rd.M; //KJ/kg K
}

public double Cp0(double t)
{
//ideal gas constant pressure heat capacity
double T=t+273.15;
double Cp0=(rd.acp[0]+rd.acp[1]*T+rd.acp[2]*T*T)/rd.M;
return Cp0;  //KJ/kg K
}

//===============================
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 Ps=Ps(t);
double vsl=vsl_t(t);
double vsg=vsg_t(t);
if(Ps >= rd.Pc*1e3)
  xx=3;
if(p < Ps)  xx=2;
else if (p > Ps) xx=-1;
else //p==Ps
{
  if(v==vsl)
    xx=0;
  else if(v==vsg)
    xx=1;
  else if(v>vsl && v<vsg)
    xx=(v-vsl)/(vsg-vsl);
}    
  return xx;  
}

public double x_pv(double p, double v)
{
//phase and/or quality of the two pahase system
double xx=0;
double t=t_pv(p,v);
double Ps=Ps(t);
double vsl=vsl_t(t);
double vsg=vsg_t(t);
if(Ps >= rd.Pc*1e3)
  xx=3;
if(p < Ps)  xx=2;
else if (p > Ps) xx=-1;
else //p==Ps
{
  if(v==vsl)
    xx=0;
  else if(v==vsg)
    xx=1;
  else if(v>vsl && v<vsg)
    xx=(v-vsl)/(vsg-vsl);
}    
return xx;  
}
//===============================
/*
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 Ts=Ts(p); 
if(t> rd.Pc)
  xx=3;
else if(t > Ts)
  xx=2;        
else if( t < Ts)
  xx=-1; 
else if( t==Ts)
 {  double vsl=vsl_t(t);
    double vsg=vsg_t(t);
    if(v==vsl)
       xx=0;
    else if(v==vsg)
       xx=1;
    else if(v>vsl && v<vsg)
       xx=(v-vsl)/(vsg-vsl);
 }
  return xx;  
}

public double x_pv(double p, double v)
{
//phase and/or quality of the steam-water system
double ts=Ts(p);
if(p > rd.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);
 }
}
*/
// ======= u relations =============



//==========toString output methods
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, iç enerji       ";
   s1[5][0]="s, entropi         ";
   s1[6][0]="x, kuruluk derecesi";
   s1[7][0]="   yoğunluk        ";
   s1[8][0]="   faz             ";
   if(unit.equals("SI"))
   {
   s1[0][2]=" kPa           ";
   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 buh/kg kar.";
   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 buh/lbm kar";
   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]=" derece C      ";
   s2[2]=" m^3/kg        ";
   s2[3]=" KJ/kg         ";
   s2[4]=" KJ/kg         ";
   s2[5]=" KJ/kg         ";
   s2[6]=" kg buh/kg kar ";
   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 buh/lbm kar";
   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]=" derece C       ";
   s2[2]=" m^3/kg         ";
   s2[3]=" KJ/kg          ";
   s2[4]=" KJ/kg          ";
   s2[5]=" KJ/kg K        ";
   s2[6]=" kg buh/kg kar. ";
   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 buh/lbm kar";
   s2[7]=" lbm/ft^3       ";
   }
for(int i=0;i<8;i++)
{
s3+=s1[i]+s2[i]+"\n";
}
return s3;
}
public String phase(double x)
{
//phase of the steam-water system
if( x==3)
  return "kritik üstü sıvı";
else
 {
 if(x<0)
   return "sıvı";
 else if( x==0.0)
   return "doymuş sıvı";
 else if(x==1.0)
   return "doymuş buhar";
 else if( x>0.0 && x<1.0)
   return "doymuş sıvı-buhar karışımı";
 else // (x>1)
   return "kızgın buhar";
 }
}
public String phase(double t, double v)
{
//sıvı-buhar sisteminin fazı
if( Ps(t) > rd.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";
 }
}
void setUnit(String x)
{
//sets SI or English units
if(x.equals("EN")) unit=u[1];
else unit=u[0];
}

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

public String readref(String rName)
{
String s="";
Vector gi=new Vector();
    rd=new sogutMBWRData("R134a","1,1,1,2-tetrafluoroethane",100.0,100.0,
    102.03,0.08314471,-26.06,
    101.108,4060.3,515.3,0.00194,
    -6.5455235227E-02,
     5.8893751817E+00,
    -1.3761788409E+02,
    2.2693168845E+04,
    -2.9262613296E+06,
    -1.1923776190E-04,
    -2.7214194543E+00,
    1.6295253680E+03,
    7.2942203182E+05,
    -1.1724519115E-04,
    8.6864510013E-01,
    -3.0660168246E+02,
    -2.5664047742E-02,
    -2.4381835971E+00,
    -3.1603163961E+02,
    3.4321651521E-01,
    -1.0154368796E-02,
    1.1734233787E+00,
    -2.7301766113E-02,
    -6.6338502898E+05,
    -6.4754799101E+07,
    -3.7295219382E+04,
    1.2614735899E+09,
    -6.4742200070E+02,
    1.2362450399E+05,
    -1.5699196293E+00,
    -5.1848932204E+05,
    -8.1396321392E-02,
    3.0325168842E+01,
    1.3399042297E-04,
    -1.5856192849E-01,
    9.0679583743E+00,
    19.4006,0.258531,-1.29665e-4,
    -2.679733546158988,-0.03529415434499572,1.341475686349758E-4,-5.143558797023507E-7,2.9575357272122354E-10,-4.88776240453033E-11,
    //5.281464E+02,7.551834E+02,1.028676E+03,-9.491172E+02,5.935660E+02,
//    536.221006032219,616.352543182115,1676.510938011078,-2096.876408754982,1282.310820150414,
518.405669486274,787.522450665451,1117.355047146022,-1342.056141309297,921.106359488362,
-94.37312278041645,
10.030208098702133,
0.7191056597512215,
0.0893279209267348,
-0.052688395138829947,
0.02434737340081483,
-0.004819940673769452,
4.7109431397984736E-4,
-1.7773917875274492E-5,
    0.0,200.0,1.0,198.36,0.7262);
    rd.vapor_pressure_in(
40.522800283269,-2357.438894824160, -0.089065965808,-13.065860972714,  0.008392230792, -0.000000055394, -0.000000002942,  0.249822319940,371.009384430072,1.0,-102.0,-50.0,
40.699854935181,-2354.790437386798,  0.004524205291,-13.059792264105,  0.007528007120, -0.000000010072,  0.000000000068,  0.232577537221,363.153073423871,1.0,-50.0 ,  0.0,
40.727590901535,-2385.172015920361, -0.020630472765,-13.077978823494,  0.007796546367,  0.000000014654, -0.000000000395,  0.261507125628,381.806492958325,1.0,  0.0 , 50.0, 
40.701201519379,-2362.752882389259,  0.000729358084,-13.067723238875,  0.007616849378,  0.000000001856, -0.000000000073,  0.234456803819,374.516259905957,1.0, 50.0 , 80.0,
40.699848569178,-2363.997386188621,  0.098879745992,-13.065540674524,  0.007607121638,  0.000000000547, -0.000000000040,  0.233150514766,375.850034550997,1.0, 80.0 , 95.0,
40.658257980565,-2469.694729012021,  0.086021390747,-12.923013853986,  0.008028696845, -0.000000041928, -0.000000003745,  0.208545513549,374.258243825424,1.0, 95.0 ,101.0);
    try{
    gi.addElement(rd);
       } catch(ArrayIndexOutOfBoundsException e) {System.out.println("Array Index Out Of Bounds Exception");}
    rd=new sogutMBWRData("R123","2,2 dichloro-1,1,1-trifluoroethane",0.9980793432775332,1.0,
    152.93,8.314471,27.85,
    183.68,3668.0,550.0,0.0018181818,
   -1.9737179010E+00,
    1.4351227567E+02,
   -2.9709876039E+03,
    4.5057867093E+05,
   -4.0472167920E+07,
    1.0191935793E-01,
   -1.9633752269E+02,
    9.2512206444E+04,
   -4.0076301858E+07,
   -2.0404730428E-02,
    3.6819648916E+01,
   -1.1445252703E+04,
   -5.5096827243E-01,
   -7.7622951034E+00,
    1.2172130502E+05,
   -8.2869171955E+00,
    7.2619253573E-01,
   -1.4363099944E+03,
    6.4252739335E+01,
    3.3524849016E+07,
   -2.0683821999E+09,
    1.8256503429E+06,
    4.2391447763E+10,
    2.6036695332E+04,
    1.3075003968E+07,
    1.1568895787E+03,
   -7.4177654079E+06,
   -2.4495004199E+00,
    4.9612079844E+03,
    4.0451631009E-03,
    1.2149154232E+02,
    2.2624786818E+04,
    28.9811,0.304711,-1.95477,
    -0.7990443910133962,-0.040602051348168036,1.4487013236958585e-4,-1.1319636530766604e-6,1.0117822724317965e-8,-3.430833666112072e-11,
    523.751847403986,1117.707515854108,220.747895457507,-37.569995382442,234.080949368258,   
   -58.94032112615241,
    13.011386898149794,
    0.8057115166593576,
    0.0690016385974559,
    0.007798621601978084,
   -0.001061467427007301,
   -3.69960220041321E-5,
    5.8996032777258733E-5,
   -4.314058230647788E-6,
    0.0,200.0,1.0,179.7,1.6580);
rd.vapor_pressure_in(1658.695405450931,-2483856.580538052600,-2.101805152254,18.033121525459,-0.088609090559,-0.000007971647,-0.000000079758,461.119202148665,1670.294704674823,1.0,-100.0,-50.0,
1656.338125951897,-2480583.209380854000,-0.001172796865,17.925084650340,-0.088692876569,-0.000000000058,-0.000000000035,461.785925991842,1666.666926993161,1.0,-50.0   ,0.0,
1656.344035650972,-2480575.541480165000, 0.000117880595,17.927174613479,-0.088688920001,-0.000000001516,-0.000000000208,461.782662717192,1666.669549145737,1.0,  0.0   ,50.0,
1656.321518615394,-2480614.220769612600,-0.002050325658,17.924621732037,-0.088686944053,-0.000000000068, 0.000000000125,461.787616118262,1666.683728384168,1.0, 50.0,  100.0,
1656.367348347346,-2480577.740957461300,-0.001676458222,17.925265842914,-0.088686418576,-0.000000000308,-0.000000000111,461.787092935390,1666.654666564091,1.0,100.0,  170.0,
1656.333544001960,-2480584.630558614600, 0.000369694724,17.925305616845,-0.088685075172, 0.000000000386, 0.000000000157,461.785793654575,1666.664478260489,1.0,170.0,  184.0);
    try{
    gi.addElement(rd);
       } catch(ArrayIndexOutOfBoundsException e) {System.out.println("Array Index Out Of Bounds Exception");}
    int i;
    Enumeration n1=gi.elements();
    while(n1.hasMoreElements())
    {   rd=(sogutMBWRData)n1.nextElement();
    s+=rd.refrigerantName+" ";
    if (rd.refrigerantName.equals(rName)) {break;}
    }
    return s;
    }
}