public class pipe_pressure_drop_TP
{ refrigerant rf;
  double eod;
  double d;
  double m;
  double L[];
  double P[];
  int n;
  double A;
  public pipe_pressure_drop_TP(String name, double di, double mi, double eodi)
  {rf=new refrigerant(name);
  d=di;
  m=mi;
  eod=eodi;
  n=100;
  L=new double[n+1];
  P=new double[n+1];
  A=Math.PI*di*di/4.0;
  }

//---- Çı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

  public void dP(double P1,double t1, double P2)
  {
   //start as one phase fluid and then continue as two phase
    double dP=(P1-P2)/n;
    double dPmomentum,dPmomentum1,dPmomentum2; 
    double dPL;
    double dPT;
    double dl;
    double x=0;
    double rol=0,rog=0;//yoğunluk
    double visl,visg;//vizkozite
    double Rel,Reg;//Reynolds sayısı
    double t=t1; //sıcaklık
    double h=0; //entalpi
    double sigma;//yüzey gerilim kuvveti
    double Vl,Vg;
    double p=P1;
    double epsilon;
    double E,F,H,WeL,fi,FrH;
    double fl,fg;
    P[0]=P1;
    L[0]=0.0;
    double l=0;
    double xo,to,ho,vislo,visgo,rolo,rogo,sigmao,epsilono;    
    double a1[]=rf.property("tp",t1,P1);
    double roH;
	xo=a1[6];
    to=a1[1];
    ho=a1[3];
    h=ho;
    vislo=rf.r.viscosity_l(t);   // sıvı vizkozite
	visgo=rf.r.viscosity_v(t);   // buhar vizkozite
	rolo=rf.r.rol(t);            // sıvı yoğunluk
	rogo=rf.r.rov(t);            // buhar yoğunluk
	sigmao=rf.r.yuzeygerilim(t); // yüzey gerilim kuvveti
    double ps=p_l(ho,P1,P2);
    double ts=rf.r.Tsb(ps);
    double hs=rf.r.h_l(ts);
    System.out.println("ps="+ps+"ts="+ts+"hs="+hs);
    double ts1,hs1;
    for(int i=0;i<n;i++)
    { 
	xo=a1[6];
    to=a1[1];
    ho=a1[3];
    vislo=rf.r.viscosity_l(t);   // sıvı vizkozite
	visgo=rf.r.viscosity_v(t);   // buhar vizkozite
	rolo=rf.r.rol(t);            // sıvı yoğunluk
	rogo=rf.r.rov(t);            // buhar yoğunluk
	sigmao=rf.r.yuzeygerilim(t); // yüzey gerilim kuvveti
    p-=dP;
    System.out.println("i="+i+"t="+t+"x="+x+"P="+p+"l="+l);
	a1=rf.property("ph",p,h);
    x=a1[6];
    t=a1[1];
    //h=a1[3];
    visl=rf.r.viscosity_l(t);   // sıvı vizkozite
	visg=rf.r.viscosity_v(t);   // buhar vizkozite
	rol=rf.r.rol(t);            // sıvı yoğunluk
	rog=rf.r.rov(t);            // buhar yoğunluk
	sigma=rf.r.yuzeygerilim(t); // yüzey gerilim kuvveti
	 if(p>ps)  //sıvı fazı
	 {   Vl=m/(rol*A);                 //  hız
	     rol=a1[7];
	     Rel=rol*Vl*d/visl;                //  sıvı Reynolds sayıs
	     double f1=f(Rel);
	     dl=dP*1e3*d*2*9.81/(rol*Vl*Vl*f1);
	     System.out.println("M="+m+"rol="+rol+"A="+A+"dP="+dP+"d="+d+"rol="+rol+"V1="+Vl+"Rel="+Rel+"f1="+f1);
	     l+=dl;    
	     L[i+1]=l;
	     P[i+1]=p;
	 }
	 else   //sıvı buhar karışımı(p<=ps)
	 {
	  epsilon=x/rog/((1+0.2*(1-x)*Math.pow((9.806*d*rol*rol/(m*m)),0.25))*(x/rog+(1-x)/rol)+1.18*(1-x)*Math.pow(9.806-sigma*(rol-rog),0.25)/(m*Math.sqrt(rol)));
	  epsilono=xo/rogo/((1+0.2*(1-xo)*Math.pow((9.806*d*rolo*rolo/(m*m)),0.25))*(xo/rogo+(1-x)/rolo)+1.18*(1-x)*Math.pow(9.806-sigma*(rolo-rogo),0.25)/(m*Math.sqrt(rolo)));
	  dPmomentum1=m*m*(((1-x)*(1-x)/(rol*(1-epsilon))+x*x/(rog*epsilon)));
	  dPmomentum2=m*m*(((1-xo)*(1-xo)/(rolo*(1-epsilono))+x*x/(rogo*epsilono)));
	  dPmomentum=dPmomentum2-dPmomentum1;
	  System.out.println("epsilon="+epsilon+"epsilono="+epsilono+"dPmomentum="+dPmomentum);	  
	  Vl=m/(rol*A);                 //  hız
	  Rel=rol*Vl*d/visl;
	  fl=f(Rel);
	  Vg=m/(rog*A);                 //  hız
	  Reg=rog*Vg*d/visg;
	  fg=f(Reg);  
	  E=(1-x)*(1-x)+x*x*rol*fg/(rog*fl);
	  F=Math.pow(x,0.78)*Math.pow((1-x),0.224);
	  H=Math.pow((rol/rog),0.91)*Math.pow((visg/visl),0.19)*Math.pow((1.0-visg/visl),0.7);
	  roH=1.0/(x/rog+(1-x)/rol);
	  WeL=m*m*d/(sigma*roH);
	  FrH=m*m/(9.806*d*roH*roH);
	  fi=E+3.2*F*H/(Math.pow(FrH,0.045)*Math.pow(WeL,0.035));
      //dPL=(dP-dPmomentum*1e-3)/fi;
      dPL=dP/fi;
      dl=dPL*1e3*d*2*9.81/(rol*Vl*Vl*fl);
      l+=dl;
	  L[i+1]=l;
	  P[i+1]=p;
	  } 
	} 

	}
  //boru içi 

	
  // eps     : pipe equivalent roughness m
  //          -------------------------------------------
  //          value of equivalent roughness for some pipes:
  //          Pipe             equivalent rougness m
  //          -------------    ---------------------
  //          Riveted steel    0.9e-3 - 9e-3
  //          Concrete         0.3e-3 - 3e-3
  //          Wood stave       0.18e-3 - 0.9 e-3
  //          Cast iron        0.26e-3
  //          Galvanised iron  0.15e-3
  //          Commercial steel 0.045e-3
  //          wrought iron     0.045e-3
  //          drawn tubing     0.0015e-3
  //          Plastic,glass    0.0 (smooth)
  //          -------------------------------------------
// boru içi 
 public double f(double Re)
 {
 //Goudar equation 
 double a=2/Math.log(10);
 double b=eod/3.7;
 double d=Math.log(10.0)/5.02*Re;
 double s=b*d+Math.log(d);
 double q=Math.pow(s,(s/(s+1)));
 double g=b*d+Math.log(d/q);
 double z=Math.log(q/g);
 double dLA=(g/(g+1))*z;
 double dCFA=dLA*(1+z/2.0/((g+1)*(g+1)+(z/3)*(2.0*g-1)));
 double f1=a*(Math.log(d/q)+dCFA);
 f1=1.0/(f1*f1);
 System.out.println("Goudar equation f1="+f1);
 return f1;
 }
 
double fps(double h, double p)
{return rf.r.h_l(rf.r.Tsb(p))-h;}

double h_tp(double p,double h,double t)
{return rf.r.h_l(t)+(rf.r.Psb(t)-p)/rf.r.rol(t)-h;}

double t_l(double P,double h,double t1,double t2)
{ 
double x1=t1,x2=t2;
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.0001;
f1=h_tp(P,h,x1);
f2=h_tp(P,h,x2);
double x3=(x1+x2)/2.0;f3=fps(h,x3);
if(f1==0) return x1;
else if(f2==0) return x2;
else if(f3==0) return x3; 
p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
fp=h_tp(P,h,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=h_tp(P,h,p);
    ea=Math.abs(x3-p);
    }
    iter++;
} 
return p;
}
 
double p_l(double h,double p1,double p2)
{ 
double x1=p1,x2=p2;
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.0001;
f1=fps(h,x1);
f2=fps(h,x2);
double x3=(x1+x2)/2.0;f3=fps(h,x3);
if(f1==0) return x1;
else if(f2==0) return x2;
else if(f3==0) return x3; 
p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
fp=fps(h,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=fps(h,p);
    ea=Math.abs(x3-p);
    }
    iter++;
} 
return p;
}
public void plot()
{ Plot pp=new Plot(L,P);
pp.plot();
}

public static void main(String arg[])
{pipe_pressure_drop_TP pd=new pipe_pressure_drop_TP("R22",1.0e-3,0.1,1.0e-5);
 pd.dP(1000.0,20.0,100.0);
 pd.plot();
}

}