
public class droplet_condensation
{
Gas1 g1=new Gas1("air");
ref_CS3 rfcs=new ref_CS3();
steam1 st=new steam1();	
yashava yh=new yashava();
moist_air ma=new moist_air();
double Tdot=1.0;//K
double Pdot=1.0;//Mpa


public double sat_temp_K(double psat)
	{// International steam tables , Springer(2008)
	 // saturation temperature equation
	 // page 26
	 // psat in Mpa
	double n1[]={0.11670521452767e4,-0.72421316703206e6,-0.17073846940092e2,0.12020824702470e5,-0.32325550322333e7,0.14915108613530e2,-0.48232657361591e4,0.40511340542057e6,-0.23855557567849,0.65017534844798e3}	;
	double beta=Math.pow((psat/Pdot),0.25);
	double E=beta*beta+n1[2]*beta+n1[5];
	double F=n1[0]*beta*beta+n1[3]*beta+n1[6];
	double G=n1[1]*beta*beta+n1[4]*beta+n1[7];
	double D=(2*G)/(-F-Math.sqrt(F*F-4*E*G));
	double dum=Math.sqrt((n1[9]+D)*(n1[9]+D)-(4*(n1[8]+n1[9]*D)));
	double Ts=(n1[9]+D-dum)/2;
	return Ts;	
	}

public double sat_pres_kPa(double Tsat)	
   {// International steam tables , Springer(2008)
	 // saturation temperature equation
	 // page 25
	 // Tsat in K
	 double n[]={0.11670521452767e4,-0.72421316703206e6,-0.17073846940092e2,0.12020824702470e5,-0.32325550322333e7,0.14915108613530e2,-0.48232657361591e4,0.40511340542057e6,-0.23855557567849,0.65017534844798e3}	;
	 double v=Tsat+(n[8]/(Tsat-n[9]));
	 double A=v*v+n[0]*v+n[1];
	 double B=n[2]*v*v+n[3]*v+n[4];
	 double C=n[5]*v*v+v*n[6]+n[7];
	 double nomin=2*C;
	 double denomin=-B+Math.sqrt(B*B-4*A*C);
	 double psat=Math.pow((nomin/denomin),4);
	 return psat*1000;//kPa  
	} 

   public double[] water_sat_prop(double psat)
	{
	double top[]=new double[17]; 
	double ts=sat_temp_K(psat);
	top[0]=rfcs.Psb(ts-273.15);
	top[1]=rfcs.rol(ts-273.15);
	top[2]=rfcs.rov(ts-273.15);
	top[3]=rfcs.h_l(ts-273.15);
	top[4]=rfcs.h_v(ts-273.15);
	top[5]=rfcs.h_lv(ts-273.15);
	top[6]=rfcs.s_l(ts-273.15);
	top[7]=rfcs.s_v(ts-273.15);
	top[8]=rfcs.viscosity_l(ts-273.15);
	top[9]=rfcs.viscosity_v(ts-273.15);
	top[10]=rfcs.k_l(ts-273.15);
    top[11]=rfcs.k_v(ts-273.15);
    top[12]=rfcs.Cpl(ts-273.15);
    top[13]=rfcs.Cpv(ts-273.15);
    top[14]=rfcs.soundv_l(ts-273.15);
    top[15]=rfcs.soundv_v(ts-273.15);
    top[16]=rfcs.yuzeygerilimi(ts-273.15);
    return top;	
	}	
	

public double h_vap(double tamb,double t,double P,double rh,double V,double D)
{
// Tbulk  Kelvin
// P bar	
// D m
double Nu=0.0;
double tbulk=(t+tamb)/2;	
double ro_vap=ma.ro_moist(tbulk,P,rh);
double[] vis=ma.vis_moist(t,P,rh);
double vis_m=vis[2];
double kv=ma.k_moist(t,P,rh);
double Re=ro_vap*V*D/vis_m;
if((Re>0)&&(Re<=3e5))
{Nu=2+Math.sqrt(0.25*Re+3e-4*Math.pow(Re,1.6));}
else if((Re>3e5)&&(Re<5e6))
{Nu=430+Re*5e-3+0.025e-9*Re*Re-3.1e-17*Re*Re*Re;}
double h=Nu*kv/D;
return h;	
}	
	



public static double [] thomas(double f[],double e[],double g[],double r[])
  {
  // 3 band matrix system of equation solving algorithm
  int n=f.length;
  double x[]=new double[n];
  for(int k=1;k<n;k++)
     {e[k]=e[k]/f[k-1];
      f[k]=f[k]-e[k]*g[k-1];
     }
  for(int k=1;k<n;k++)
     {r[k]=r[k]-e[k]*r[k-1];
     }
  x[n-1]=r[n-1]/f[n-1];   
  for(int k=(n-2);k>=0;k--)
     {x[k]=(r[k]-g[k]*x[k+1])/f[k];}
  return x;      
  }
  
  
  public static double [][] cubic_spline(double xi[],double yi[],double c0,double cn)
  {
  int n=xi.length;
  double h[]=new double[n];
  double w[]=new double[n];
  double f[]=new double[n];
  double e[]=new double[n];
  double g[]=new double[n];
  double d[]=new double[n];
  double x[]=new double[n];
  double S[][]=new double[4][n];
  int k;	  
  for(k=0;k<(n-1);k++)
     {h[k]=xi[k+1]-xi[k];
      w[k]=(yi[k+1]-yi[k])/h[k];      
	 }
  d[0]=c0;
  d[n-1]=cn;
  for(k=1;k<(n-1);k++)
     {d[k]=6.0*(w[k]-w[k-1]);}
  f[0]=1.0;
  f[n-1]=1.0;
  g[0]=0.0;
  g[n-1]=0.0;
  e[0]=0.0;
  e[n-1]=0.0;
  for(k=1;k<(n-1);k++)
  {f[k]=2.0*(h[k]+h[k-1]);e[k]=h[k-1];g[k]=h[k];}
  S[2]=thomas(f,e,g,d);
  S[3]=xi;
  for(k=0;k<(n-1);k++)
    {S[0][k]=(6.*yi[k+1]-h[k]*h[k]*S[2][k+1])/(6.0*h[k]);
     S[1][k]=(6.*yi[k]-h[k]*h[k]*S[2][k])/(6.0*h[k]);
    } 
  return S;    	 	     	 
  }
  
  public static double funcSpline(double S[][],double x)
  {
  int n=S[0].length;
  double xx1=0;
  double xx2=0;
  double y=0;
  double hk=0;
  for(int k=0;k<(n-1);k++)
     {if(S[3][k]<=x && x<=S[3][k+1]) 
      {hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       break;
       }
     }
       if(y==0 && S[3][n-2]<=x )
       {    
	   int k=n-2;
	   hk=(S[3][k+1]-S[3][k]);
       xx1=(x-S[3][k]);
       xx2=(S[3][k+1]-x);
       y=S[0][k]*xx1+S[1][k]*xx2+(xx1*xx1*xx1*S[2][k+1]+xx2*xx2*xx2*S[2][k])/(6.0*hk);
       }
  return y;
  }     

public double kl_val(double l_h)
{
double[] lh={20,40,60,80,100,200};	
double[] val={1.56,1.18,1.08,1.04,1.00,1.00};
double[][] roint=cubic_spline(lh,val,0,0);
double newkl=funcSpline(roint,l_h); 	
return newkl;	
} 
 

public double h_cond(double tdeww,double tamb,double P,double rh,double l_h,double V,double D)
{

double k_tp=ma.k_moist(tdeww,P,1.0);	
double[] vis=ma.vis_moist(tdeww,P,1.0);	
double vis_tp=vis[2];
double ro_tp=ma.ro_moist(tdeww,P,1.0)	;
double Re=ro_tp*V*D/vis_tp;
double Nu=0.065*Math.pow(Re,0.66)*kl_val(l_h);	
return Nu*k_tp/D;
}
	
public double[] tekfaz_t(double tamb,double t,double P,double rh,double V,double D,double Vroom,double m_su,double total_length)
{
double[] prop=yh.property("tdb_rh",tamb,rh,P);
double w=prop[5];
double twp=prop[10];//
double Pair=prop[4]*100;//kpa
double mair=(Pair*Vroom)/(0.287*(tamb+273));// havanın kütlesi
double m_yh=mair*w;//yaş havanın kütlesi
double tt=0;
double dt=0.001;//sec
double h;
double A=Math.PI*D*D;//damlacık yüzey alanı
double dq=0.0;
double h_ent=ma.enthalpy_moist(tamb,P,rh);//kj/kg moist air
double h_droplet=rfcs.h_l(t);//kj/kg // su damlacığı entalpi
double prp[];
double Vol=(4/3)*Math.PI*(D/2)*(D/2)*(D/2);//su damlacığı hacmi
double mdroplet=0.0;
double psat;//=sat_pres_kPa(t+273.15);
double pp[];
double n_droplet;
double xc=0.0;
int i=0;

double[] dist=new double[50000];
double[] hdroplet=new double[50000];
double[] tdroplet=new double[50000];
double[] timecounter=new double[50000]; 
double[] tambient=new double[50000];
double[] hentalphy=new double[50000];
while(tamb>twp)
{
tambient[i]=tamb;	
timecounter[i]=tt;
hdroplet[i]=h_droplet;
tdroplet[i]=t;
dist[i]=(V*timecounter[i])+(0.5*9.81*timecounter[i]*timecounter[i]);
xc=(V*timecounter[i])+(0.5*9.81*timecounter[i]*timecounter[i]);
if(xc>total_length)
{
System.out.println("Damlacık toplam yüksekliğe ulaştı, katedilen yol "+xc+" m");
break;	
}
h=h_vap(tamb,t,P,rh,V,D);
dq=h*A*(tamb-t)*(dt-tt);	
hentalphy[i]=h_ent;
mdroplet=rfcs.rol(t)*Vol;
n_droplet=m_su/mdroplet;
h_droplet=((dq/mdroplet)+h_droplet*1000)/1000;
psat=sat_pres_kPa(t+273.15);
pp=st.property_SI("ph",psat,h_droplet);
t=pp[1];
System.out.print("tdew= "+twp+" tamb= "+tambient[i]+" zaman= "+timecounter[i]+" yaş hava entalpi= "+hentalphy[i]);
System.out.println("  katedilen yol= "+dist[i]+" t_droplet= "+tdroplet[i]);
h_ent=((-dq*n_droplet/mair)+h_ent*1000)/1000;
prp=yh.property("w_h",w,h_ent,P);
tamb=prp[1];
dt+=0.001;
tt+=0.001;
i++;
}
double xx=(V*timecounter[i-1])+(0.5*9.81*timecounter[i-1]*timecounter[i-1]);
System.out.println("Damlacığın aldığı yol "+xx+" m");
double[] depo=new double[7]; 
depo[0]=timecounter[i-1];
depo[1]=twp;
depo[2]=hdroplet[i-1];
depo[3]=tdroplet[i-1];
depo[4]=xx;
depo[5]=tambient[i-1];
depo[6]=mair;
return depo;	
}

public void two_phase_iter(double tamb,double t,double P,double rh,double V,double D,double Vroom,double m_su,double total_length,double l_h)
{ 
double[] values=tekfaz_t(tamb,t,P,rh,V,D,Vroom,m_su,total_length);
double time=values[0];//sn
double tdew=values[1];// C
double droplet_ent=values[2];//kj/kg
double temp_droplet=values[3];// C
double distance=values[4];// m	
double tambient=values[5];// C
double mair=values[6];//kg
double tt=0.0;
double dt=0.001;
double pp[];
double dq;
double h_cnd;
double dm;
double psat;
double mdroplet;
double Volume=(4/3)*(Math.PI)*(D/2)*(D/2)*(D/2);
double A=Math.PI*D*D;
double h_moistair=ma.enthalpy_moist(tdew,P,1.0);//kj/kg rh=1.0!!!
double hfg_moistair=rfcs.h_lv(tdew)*1000;
double ndroplet;

double[] dent_counter=new double[10000];
double[] dtemp_counter=new double[10000];
double[] time_counter=new double[10000];
double[] dist_counter=new double[10000];
double[] hma_counter=new double[10000];

int i=0;
while(total_length>distance)
{
	
hma_counter[i]=h_moistair;
dist_counter[i]=distance;
time_counter[i]=tt;
dtemp_counter[i]=temp_droplet;
dent_counter[i]=droplet_ent;
System.out.println(" Katedilen yol= "+dist_counter[i]+" zaman="+time_counter[i]+" Damlacık sıcaklığı= "+dtemp_counter[i]+" Yaş hava entalpi= "+hma_counter[i]);
h_cnd=h_cond(tdew,tamb,P,rh,l_h,V,D);
dq=h_cnd*A*(tdew-temp_droplet)*(dt-tt);
dm=dq/hfg_moistair;
mdroplet=rfcs.rol(temp_droplet)*Volume;
droplet_ent=(mdroplet*droplet_ent+dm*h_moistair)/(mdroplet+dm);
psat=sat_pres_kPa(temp_droplet+273.15);
pp=st.property_SI("ph",psat,droplet_ent);
temp_droplet=pp[1];
ndroplet=m_su/mdroplet;
h_moistair=((-dq*ndroplet/mair)+h_moistair*1000)/1000;
distance+=((V*tt)+(0.5*9.81*tt*tt));
tt+=0.001;
dt+=0.001;
i++;
}
System.out.println("Toplam zaman= "+(time+time_counter[i-1])+" sn");
}

public static void main(String args[])
{
droplet_condensation dc=new droplet_condensation();
double t=11.0;// C ilk damlacık sıcaklığı
double tamb=80.0;// C ilk nemli hava sıcaklığı
double P=1.0;//bar nemli havanın basıncı
double rh=0.55;//bağıl nem
double V=0.5;//m/s damlacık hızı
double D=1e-4;// m damlacık çapı
double l_h=40;
double m_su=0.2;//kg su(tepeden püskürtülen) 
double Vroom=1.5;//m3 ortam hacmi
double total_length=5.5;// m püskürtme yüksekliği

dc.two_phase_iter(tamb,t,P,rh,V,D,Vroom,m_su,total_length,l_h);
}	
	
}