///////////////////////////////////////////////////////////////////////
// Aqueous Solutions of Lithium and Calcium Chlorides
// Property formulatýons : For use in air conditioning equipment design
// M.Conde Engineering,2009
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
// Nomenclature////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
// eps=mass fraction of solute/////////////////////////////////////////
// teta=reduced temperature of water///////////////////////////////////
///////////////////////////////////////////////////////////////////////
public class LiClprop
{
double rocw=322.0;//critical density of water[kg /m3]	
double Tcw=647.226;// Critical temperature of H2O [K]	
double Pcritw=220.9e5;//critical pressure of H2O [Pa]	
double dvisc=55.071e-6;// Critical dynamic visc of H20[Pa.s]
double[] A_ice={0.422088,-0.090410,-2.936350};	
double[] A_5H2O={-0.005340,2.015890,-3.114590};
double[] A_3H2O={-0.560360,4.723080,-5.811050};
double[] A_2H2O={-0.315220,2.882480,-2.624330};
double[] A_H2O={-1.312310,6.177670,-5.034790};
double[] A_LiCl={-1.356800,3.448540,0.0};	
double[] NNL={0.0,0.0,1.0,0.0,-1.0,-2.0,-3.0,-4.0,1.0,0.0,-1.0,-2.0,-3.0};
double[] MML={0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0,1.0};
double[] AAL={0.0,5.88137e1,-9.18938,-8.33307e1,1.84228e2,-1.70534e2,-1.05457e2,-2.94856e1,1.96926e-1,-8.68008e-1,2.25497,-2.60177,1.17349};
double[] NNG={0.0,0.0,1.0,0.0,-3.0,6.0,7.0,8.0,15.0,0.0,0.0,0.0,0.0};
double[] MMG={0.0,0.0,0.0,0.0,0.0,1.0,1.0,1.0,3.0,0.0,0.0,0.0,0.0};
double[] AAG={1.0,3.98432,6.68909,-8.78439,-5.18911e-2,-6.25248e-1,5.07144e-1,-1.24364e-1,-4.49013e-1,0.0,0.0,0.0,0.0};


// Converts mass fraction to mole fraction
public double mass_to_mole(double w)
{
double M_LiCl=42.394;//kg/kmol
double M_H2O=18.015;//kg/kmol
return (w/M_LiCl)/((w/M_LiCl)+((1-w)/M_H2O));	
}






//Calculation of Solubility boundary 
//consist of Ice Line - LiCl_5H20 - LiCl_3H20 - LiCl_2H20 - LiCl_H20 - LiCl
public  double sol_boundary(double eps)
{
  double sumlicl=0.0;
  double sumH20=0.0;
  double sum2H20=0.0;
  double sum3H20=0.0;
  double sum5H20=0.0;
  double sumice=0.0;

        for(int i=0;i<=2;i++)
        {sumlicl+=A_LiCl[i]*Math.pow(eps,i);
         sumH20+=A_H2O[i]*Math.pow(eps,i);
         sum2H20+=A_2H2O[i]*Math.pow(eps,i);
         sum3H20+=A_3H2O[i]*Math.pow(eps,i);
         sum5H20+=A_5H2O[i]*Math.pow(eps,i);
        }
        sumice=A_ice[0]+A_ice[1]*eps+A_ice[2]*Math.pow(eps,2.5);	
	
double summice=(sumice*Tcw)-273.15;
double summlicl=(sumlicl*Tcw)-273.15;
double summH20=(sumH20*Tcw)-273.15;
double summ2H20=(sum2H20*Tcw)-273.15;
double summ3H20=(sum3H20*Tcw)-273.15;
double summ5H20=(sum5H20*Tcw)-273.15;
	
double out=0.00;
if((eps>0)&&(0.253>eps))
{out=summice;}
else if((eps>=0.253)&&(0.287>eps))
{out =summ5H20;}
else if((eps>=0.287)&&(0.369>eps))
{out =summ3H20;}
else if((eps>=0.369)&&(0.452>eps))
{out=summ2H20;}
else if((eps>=0.452)&&(0.558>eps))
{out=summH20;}
else if((eps>=0.558)&&(0.6>eps))
{out=summlicl;}
return out;
}

public void sol_plot() // plots solution boundary vs mass fraction
{
double[] xx=new double[600];
double[] yy=new double[600];
double j=0.000;
for(int i=0;i<600;i++)
{xx[i]=j;j=j+0.001;}	
for(int i=0;i<600;i++)
{yy[i]=sol_boundary(xx[i]);}
Plot pp=new Plot(xx,yy);
pp.plot();	
pp.setXlabel("Mass fraction of LiCl");
pp.setYlabel("Temperature [C]");
pp.setPlabel("Solubility boundary");
} 

public double vap_press_water(double T)
{
double A[]={-7.858230,1.839910,-11.781100,22.670500,-15.939300,1.775160};	
double to=1-(T/Tcw);
return Pcritw*Math.exp((A[0]*to+(A[1]*Math.pow(to,1.5))+(A[2]*to*to*to)+(A[3]*Math.pow(to,3.5))+(A[4]*to*to*to*to)+(A[5]*Math.pow(to,7.5)))/(1-to));	
	//in pascal
}

public double vapor_press(double eps,double T)
{
// T in kelvin
	
	
double sumlicl=0.0;
double sumH20=0.0;
double sum2H20=0.0;
double sum3H20=0.0;
double sum5H20=0.0;
double sumice=0.0;

for(int i=0;i<=2;i++)
{sumlicl+=A_LiCl[i]*Math.pow(eps,i);
 sumH20+=A_H2O[i]*Math.pow(eps,i);
 sum2H20+=A_2H2O[i]*Math.pow(eps,i);
 sum3H20+=A_3H2O[i]*Math.pow(eps,i);
 sum5H20+=A_5H2O[i]*Math.pow(eps,i);
}
sumice=A_ice[0]+A_ice[1]*eps+A_ice[2]*Math.pow(eps,2.5);	

double teta=0.00;
if((eps>0)&&(0.253>eps))
{teta=sumice;}
else if((eps>=0.253)&&(0.287>eps))
{teta =sum5H20;}
else if((eps>=0.287)&&(0.369>eps))
{teta =sum3H20;}
else if((eps>=0.369)&&(0.452>eps))
{teta=sum2H20;}
else if((eps>=0.452)&&(0.558>eps))
{teta=sumH20;}
else if((eps>=0.558)&&(0.6>eps))
{teta=sumlicl;}	
double tetaa=T/Tcw;
double[] pi ={0.28,4.30,0.60,0.21,5.10,0.49,0.362,-4.75,-0.40,0.03};	
double A=2.0-(Math.pow((1.0+(Math.pow((eps/pi[0]),pi[1]))),pi[2]));	
double B=Math.pow((1.0+(Math.pow((eps/pi[3]),pi[4]))),pi[5])-1.0;	
double x1=Math.pow((1.0+Math.pow((eps/pi[6]),pi[7])),pi[8]);
double x2=pi[9]*Math.exp(-(eps-0.1)*(eps-0.1)/0.005);
double pi25=1-x1-x2;
double f=A+B*tetaa;
double ff=pi25*f;
double vpwat=vap_press_water(T);//Pascal
return ff*vpwat/1000; // return in kPa!!!!
}

public void vap_press_plot(double T)
{
   double yy[]=new double[600];	
   double xx[]=new double[600];
   double eps=0.000;
   for(int i=0;i<600;i++)
   {
   xx[i]=eps; 	   	
   yy[i]=vapor_press(eps,T);
   eps=eps+0.001;
   }
   Plot pp=new Plot(xx,yy);
   pp.plot();
  pp.setXlabel("Mass fraction of LiCl");
  pp.setYlabel("Vapor Pressure [kPa]");
  pp.setPlabel("Vapor Pressure vs Mass Fraction at "+T+" Kelvin");
}

public double sol_density(double eps,double T)
{
double t=T-273.15;
double teta=T/Tcw;
double to=1-teta;	
double ro[]={1.0,0.540966,-0.303792,0.100791};
double B[]={1.9937718430,1.0985211604,-0.5094492996,-1.7619124270,-45.9005480267,-723692.2618632};	
double ro_sat=rocw*(1+(B[0]*Math.pow(to,0.3333))+(B[1]*Math.pow(to,0.666))+(B[2]*Math.pow(to,1.666))+(B[3]*Math.pow(to,16/3))+(B[4]*Math.pow(to,43/3))+(B[5]*Math.pow(to,110/3)));	
double s1=0.0;
ref_CS3 recs33=new ref_CS3();
double ro_satt=recs33.rol(t);
for(int i=0;i<=3;i++)
{s1+=ro[i]*Math.pow((eps/(1.0-eps)),i);}	
return  ro_sat*s1;	
}

public void sol_density_plot(double T)
{
double xx[]=new double[600];	
double yy[]=new double[600];	
double eps=0.000;
for(int i=0;i<600;i++)
{
xx[i]=eps;	
yy[i]=sol_density(eps,T);
eps=eps+0.001;	
}
Plot pp=new Plot(xx,yy);
pp.plot();
pp.setXlabel("Mass fraction of LiCl");
pp.setYlabel("Solution density [kg/m3]");
pp.setPlabel("Solution density vs Mass fraction at "+T+" Kelvin");
}

public double sol_surf_tens(double eps,double T)
{
double teta=T/Tcw;	
double s[]={2.757115,-12.011299,14.751818,2.443204,-3.147739};	
double s0=235.8e-3;//N/m
double b=-0.625;
double mu=1.256;	
double sH20=s0*(1.0-b*(1.0-teta))*Math.pow((1-teta),mu);
return sH20*(1+(s[0]*eps)+(s[1]*eps*teta)+(s[2]*eps*teta*teta)+(s[3]*eps*eps)+(s[4]*eps*eps*eps));	
}

public void sol_surf_tens_plot(double T)
{
double xx[]=new double[600];	
double yy[]=new double[600];	
double eps=0.000;
for(int i=0;i<600;i++)
{
xx[i]=eps;	
yy[i]=sol_surf_tens(eps,T);
eps=eps+0.001;	
}
Plot pp=new Plot(xx,yy);
pp.plot();
pp.setXlabel("Mass fraction of LiCl");
pp.setYlabel("Surface tension [N/m]");
pp.setPlabel("Surface tension vs Mass fraction at "+T+" Kelvin");
}

public double d_vis_sol(double eps,double T)
{
double Tsolt=T-273.15;	
double out=0.0;
ref_CS3 recs33=new ref_CS3();
double ro50=recs33.rol(0.01);
//System.out.println(ro50);	
if(Tsolt>0)
{
	
double[] H={1.000,0.978197,0.579829,-0.202354};
double[][] G=new double[6][7];
G[0][0]=0.5132047;
G[0][1]=0.2151778;
G[0][2]=-0.2818107;
G[0][3]=0.1778064;
G[0][4]=-0.0417661;
G[0][5]=0.0;
G[0][6]=0.0;
G[1][0]=0.3205656;
G[1][1]=0.7317883;
G[1][2]=-1.070786;
G[1][3]=0.4605040;
G[1][4]=0.0;
G[1][5]=-0.01578386;
G[1][6]=0.0;
G[2][0]=0.0;
G[2][1]=1.241044;
G[2][2]=-1.263184;
G[2][3]=0.2340379;
G[2][4]=0.0;
G[2][5]=0.0;
G[2][6]=0.0;
G[3][0]=0.0;
G[3][1]=1.476783;
G[3][2]=0.0;
G[3][3]=-0.4924179;
G[3][4]=0.1600435;
G[3][5]=0.0;
G[3][6]=-0.003629481;
G[4][0]=-0.7782567;
G[4][1]=0.0;
G[4][2]=0.0;
G[4][3]=0.0;
G[4][4]=0.0;
G[4][5]=0.0;
G[4][6]=0.0;
G[5][0]=0.1885447;
G[5][1]=0.0;
G[5][2]=0.0;
G[5][3]=0.0;
G[5][4]=0.0;
G[5][5]=0.0;
G[5][6]=0.0;
double royt=recs33.rol(Tsolt);

double Tr=T/Tcw;
double ror=royt/rocw;

double no=0.0;
for(int i=0;i<=3;i++)	
{no+=H[i]*Math.pow(Tr,-i);}
no=Math.sqrt(Tr)*(1/no);	

double allsum=0.0;
for(int i=0;i<=5;i++)
{
    for(int j=0;j<=6;j++)
    {allsum+=G[i][j]*Math.pow(((1.0/Tr)-1.0),i)*Math.pow((ror-1.0),j);}	
}

double n11=Math.exp(ror*allsum);	

double n2=1.0;
double nred=no*n11*n2;
double ndvis=nred*dvisc;///////////

double teta=T/Tcw;
double[] n1={0.090481,1.390262,0.675875,-0.583517};
double etta=eps/(Math.pow((1-eps),(1/0.6)));
double nsol=ndvis*Math.exp((n1[0]*Math.pow(etta,3.6))+(n1[1]*etta)+(n1[2]*etta/teta)+(n1[3]*etta*etta));
out=nsol;
} 
if(Tsolt<=0.0){

double A=1.0261862;
double B=12481.702;
double C=-19510.923;
double D=7065.286;
double E=-395.561;
double F=143922.996;
//double k0=recs33.viscosity_l(0.01);

double tetta=(T/228.0)-1.0;
double nsol0=(1750e-6)*(A+(B*Math.pow(tetta,0.02))+(C*Math.pow(tetta,0.04))+(D*Math.pow(tetta,0.08))+(E*Math.pow(tetta,2.85))+(F*Math.pow(tetta,8.0)));
out=nsol0;
}

return out;
}

public void d_vis_sol_plot(double T)
{
double xx[]=new double[500];	
double yy[]=new double[500];
double eps=0.000;

for(int i=0;i<500;i++)
{
xx[i]=eps;
yy[i]=d_vis_sol(eps,T);
eps=eps+0.001;	
}	
Plot pp=new Plot(xx,yy);
pp.plot();	
pp.setXlabel("Mass fraction of LiCl");
pp.setYlabel("Dynamic viscosity [Pa.s]");
pp.setPlabel("Dynamic viscosity vs Mass fraction at "+T+" Kelvin");
}


public double k_sol(double eps,double T)
{
double t=T-273.15;	
double alfa[]={10.8958e-3,-11.7882e-3};	
double alfar=alfa[0]+alfa[1]*eps;
double Is=1.0;
double M=42.394;
double soldens=sol_density(eps,T);	
double eps_eq=eps*soldens*Is/M;	
	
ref_CS3 recs33=new ref_CS3();
double kl=recs33.k_l(t);	
double val=(kl-eps_eq*alfar);
return val;	
} 

public void k_sol_plot(double T)
{
double xx[]=new double[500];	
double yy[]=new double[500];
double eps=0.000;

for(int i=0;i<500;i++)
{
xx[i]=eps;
yy[i]=k_sol(eps,T);
eps=eps+0.001;	
}	
Plot pp=new Plot(xx,yy);
pp.plot();	
pp.setXlabel("Mass fraction of LiCl");
pp.setYlabel("Thermal conductivity [W/mK]");
pp.setPlabel("Thermal conductivity vs Mass fraction at "+T+" Kelvin");
}
public double Cp(double eps,double T)
{
double t=T-273.15;
ref_CS3 recs33=new ref_CS3();
double cp=recs33.Cpl(t)	;
double A=1.43980;
double B=-1.24317;
double C=-0.12070;
double D=0.12825;
double E=0.62934;
double F=58.5225;
double G=-105.6343;
double H=47.7948;
double f1=0.0;
double f2=0.0;
if(eps<=0.31)
{f1=A*eps+B*eps*eps+C*eps*eps*eps;}
else
{f1=D+E*eps;}
double teta=(T/228.0)-1.0;
f2=F*Math.pow(teta,0.02)+G*Math.pow(teta,0.04)+H*Math.pow(teta,0.06);
cp=cp*(1-f1*f2);
return cp*1000;	
}

public void Cp_sol_plot(double T)
{
double xx[]=new double[600]; 	
double yy[]=new double[600]; 
double eps=0.000;
for(int i=0;i<600;i++)
{
xx[i]=eps;
yy[i]=Cp(eps,T);
eps=eps+0.001;	
}	
Plot pp=new Plot(xx,yy);
pp.plot();	
pp.setXlabel("Mass fraction of LiCl");
pp.setYlabel("Specific Heat [J/kgK]");
pp.setPlabel("Specific Heat vs Mass fraction at "+T+" Kelvin");	
}

public double differantial_enthalpy(double eps,double T)
{
double teta=T/Tcw;	
double H[]={0.845,-1.965,-2.265,0.6,169.105,457.850};
double eta=eps/(H[3]-eps);
double deltahd=H[4]+H[5]*teta;
return deltahd*Math.pow((1.0+Math.pow((eta/H[0]),H[1])),H[2]);
}

public void diff_ent_plot(double T)
{
double xx[]=new double[600]; 	
double yy[]=new double[600]; 
double eps=0.000;
for(int i=0;i<600;i++)
{
xx[i]=eps;
yy[i]=differantial_enthalpy(eps,T);
eps=eps+0.001;	
}		
Plot pp=new Plot(xx,yy);
pp.plot();
pp.setXlabel("Mass fraction of LiCl");
pp.setYlabel("Differantial Enthalpy of Dilution [kJ/kg]");
pp.setPlabel("Differantial Enthalpy of Dilution vs Mass fraction at "+T+" Kelvin");	
}

// Diffusion Coefficient of water into solutions of LiCl 
public double diff_coeff(double eps,double T)
{
double A=0.11353e-16;	
double t=T-273.15;
ref_CS3 recs33=new ref_CS3();
double vis_l=recs33.viscosity_l(t);
double roll=recs33.rol(t);
double Ru=8314.4621;//J/kmolK
double M=18.01528;//kg/kmol
double Vmolar=M/roll;//Molar volume at t [m3/kmol];
double Vcritic=M/rocw ;// Critical molar volume at ctitical temp [m3/kmol]
double Do=A*Math.pow(Vcritic,0.6666)*Ru*T/(vis_l*Vmolar);	
double eps1=0.52;
double eps2=-4.92;
double eps3=-0.56;
double DD=(1.0-Math.pow((1.0+Math.pow((Math.sqrt(eps)/eps1),eps2)),eps3));
return Do*DD;
}


public void diff_coeff_plot(double T)
{
double xx[]=new double[600]; 	
double yy[]=new double[600]; 
double eps=0.000;
for(int i=0;i<600;i++)
{
xx[i]=eps;
yy[i]=diff_coeff(eps,T);
eps=eps+0.001;	
}		
Plot pp=new Plot(xx,yy);
pp.plot();	
pp.setXlabel("Mass fraction of LiCl");
pp.setYlabel("Diffusion coefficient of water into solution of LiCl ");
pp.setPlabel("Diffusion coefficient of water into solution of LiCl vs Mass fraction at "+T+" Kelvin");		
}




public static void main(String args[])
{

double T=353.00;//K 
LiClprop LC=new LiClprop();
//LC.vap_press_plot(T);
//LC.sol_plot();
//LC.sol_density_plot(T);
//LC.d_vis_sol_plot(T);
//LC.Cp_sol_plot(T);
//LC.diff_ent_plot(T);
//LC.diff_coeff_plot(T);
//LC.k_sol_plot(T);



}















} 