import java.io.*; 
public class Darcy_friction
{
 
 public static double f_Haaland(double Re,double eod)
 {
 //Haaland equation
 double f1=-1.8*Math.log10(Math.pow((eod/3.7),1.11)+6.9/ Re);
 f1=1.0/(f1*f1);
 return f1;
 }
 
 public static double f_Moody(double Re,double eod)
 {
 // Moody equation
 // 4000<Re<107 and e/D <0.01
 double f1=5.5e-3*(1+Math.pow((2e4*eod+1e6/Re),(1.0/3.0)));
 return f1;
 }
 
 public static double f_Wood(double Re,double eod)
 {
 // Wood equation
 // Re>10000 and any e/D
 double a=0.52*eod+0.094*Math.pow(eod,0.225);
 double b=88.0*Math.pow(eod,0.44);
 double C=1.62*Math.pow(eod,0.134);
 double f1=a+b*Math.pow(Re,-C);
 return f1;
 }
 
 public static double f_Churchill(double Re,double eod)
 {
 // Churchill equation
 // for all values of Re and e/D
 double A=Math.pow((-2.457*Math.log(Math.pow((7.0/Re),0.9)+0.27*eod)),16);
 double B=Math.pow((37530.0/Re),16);
 double f1=8.0*Math.pow((Math.pow((8.0/Re),12)+1.0/Math.pow((A+B),1.5)),(1.0/12.0));
 return f1;
 }
 
 public static double f_Chen(double Re,double eod)
 {
 // Chen equation
 // for all values of Re and e/D
 double A = Math.log10(Math.pow(eod,1.1098)/2.8257 + (5.8506/Math.pow(Re,0.8981)));
 double f1=-2.0*Math.log10(eod/3.7065-5.0452*A/Re);
 f1=1.0/(f1*f1);
 return f1;
 }
 
 public static double f_Swamee(double Re,double eod)
 {
 //Swamee-Jain equation
 double A=Math.log10(eod/3.7+5.74/Math.pow(Re,0.9));
 double f1=0.25/(A*A);
 return f1;
 }
 
 public static double f_Zigrand(double Re,double eod)
 {
 //Zigrang and Sylvester Equation
 // for 4000<Re<108  and 0.00004<e/D<0.05
 double A=Math.log10(eod/3.7-5.02/Re*Math.log10(eod/3.7+13/Re));
 double f1=-2.0*Math.log10((eod)/3.7-5.02*A/Re);
 f1=1.0/(f1*f1);
 return f1;
 }
 
 public static double f_Romeo(double Re,double eod)
 {
 // Romeo - Royo - Monzon Equation (2002)

 double A=Math.log10(eod/3.827-4.567/Re*Math.log10(Math.pow((eod/7.7918),0.9924)+Math.pow(5.3326/(208.815+Re),0.9345)));
 double f1=-2.0*Math.log10((eod)/3.7065-5.0272*A/Re);
 f1=1.0/(f1*f1);
 return f1;
 }
 
 public static double f_Serghides(double Re,double eod)
 {
 //Serghides equation
 // for Re>2100 and any e/D
 double A1=-2.0*Math.log10(eod/3.7+12.0/Re);
 double B1=-2.0*Math.log10(eod/3.7+2.51*A1/Re);
 double C1=-2.0*Math.log10(eod/3.7+2.51*B1/Re);
 double f1=A1-((B1-A1)*(B1-A1))/(C1-2.0*B1+A1);
 f1=1.0/(f1*f1);
 return f1;
 }
 
 public static double f_Goudar(double Re,double eod)
 {
 //Goudar Sonnad equation 
 // for Re>2100 and any e/D
 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);
 return f1;
 }
 
 public static double f_turhan(double Re,double eod)
 {
 //Turhan equation 
      double ai[]={-0.522742801770317,	        
-0.862597376548066,	         
0.623145107999888,	        
-0.090643514190191,	        
-0.312237235469110,	         
0.627834486425472,	         
3.304936649648337};

      double f1=ai[0]*Math.log10(Math.pow(eod,ai[1]))+
             ai[2]*Math.log10(Math.pow(Re,ai[3] ))+
             ai[4]*Math.log10(Math.pow(eod,ai[5])*Math.pow(Re,ai[6] ));
             f1=1.0/(f1*f1);
      return f1;
 }
 
public static double fx(double X,double Re,double eod)
{
// colebrook equation to solve
// friction factor for turbulent flow 2000 < Re
double xx=2.0*Math.log10(eod/3.7+2.51/Re*X)+X;
return xx;
}

public static double dfx(double X,double Re,double eod)
{
//derivative of colebrook equation
double xx;
xx = 1+2.0/(eod/3.7+2.51/Re*X)/Math.log(10.0)*2.51/Re;
return  xx;
}

public static double f_colebrook(double Re,double eod)
{
// solution of the colebrook equation
// by using newton method
double fi=f_Goudar(Re,eod);
double xx=1.0/ Math.pow(fi,0.5);
int nmax=50;
double tolerance=1.0e-20;
for(int i=0;i<nmax;i++)
{
double fx1=fx(xx,Re,eod);
double dfx1=dfx(xx,Re,eod);
xx-=fx1/dfx1;
if(Math.abs(fx1)<tolerance) 
{double ff=1.0/(xx*xx);
return ff;}
}
double ff=1.0/(xx*xx);
return ff;
}

public static double[][][] calculate()
{ double eod[]={1e-8,1e-6,1e-4,1e-2,1e-1};
  double Re[]=new double[51];
  double dx1=(10000.0-2100)/10.0;
  double dx2=(100000.0-10000.0)/10.0;
  double dx3=(1e6-1e5)/10.0;
  double dx4=(1e7-1e6)/10.0;
  double dx5=(1e7-1e6)/10.0;
  double dx[]={dx1,dx2,dx3,dx4,dx5,dx5};
  double xx=0;
  int i=0,j=0,k=0,l=0;
  int n=eod.length;
  double a[][][]=new double[12][n][51];
  double x[]=new double[51];
  xx=2100.0;
  for(i=0;i<51;i++)
  {  l=i/10;
     x[i]=xx;  
     xx+=dx[l]; 
     //System.out.println("i="+i+"dx="+dx[l]+"x="+x[i]); 
  } 
  //System.out.println("x=\n"+Matrix.toStringT(x)); 
  for(k=0;k<12;k++)
  {
  for(j=0;j<n;j++)
  {
  for(i=0;i<51;i++)
  {
  if(k==0) a[k][j][i]=f_colebrook(x[i],eod[j]);
  else if(k==1) a[k][j][i]=f_Haaland(x[i],eod[j]);
  else if(k==2) a[k][j][i]=f_Moody(x[i],eod[j]);
  else if(k==3) a[k][j][i]=f_Wood(x[i],eod[j]);
  else if(k==4) a[k][j][i]=f_Churchill(x[i],eod[j]);
  else if(k==5) a[k][j][i]=f_Chen(x[i],eod[j]);
  else if(k==6) a[k][j][i]=f_Swamee(x[i],eod[j]);
  else if(k==7) a[k][j][i]=f_Zigrand(x[i],eod[j]);
  else if(k==8) a[k][j][i]=f_Serghides(x[i],eod[j]);
  else if(k==9) a[k][j][i]=f_Goudar(x[i],eod[j]);
  else if(k==10) a[k][j][i]= f_Romeo(x[i],eod[j]);
  else if(k==11) a[k][j][i]= f_turhan(x[i],eod[j]);
  }//end of i Re
  }//end of j eod
  }//end of k(equation)
  return a;
}
public static double[][][] error()
{ 
  double eod[]={1e-8,1e-6,1e-4,1e-2,1e-1};
  int n=eod.length;
  double a[][][]=calculate();
  double error[][][]=new double[12][n+1][52];
  int i=0,j=0,k=0,l=0;
  double xx=0;
  double x[]=new double[51];
  double dx1=(10000.0-2100)/10.0;
  double dx2=(100000.0-10000.0)/10.0;
  double dx3=(1e6-1e5)/10.0;
  double dx4=(1e7-1e6)/10.0;
  double dx5=(1e7-1e6)/10.0;
  double dx[]={dx1,dx2,dx3,dx4,dx5,dx5};
  xx=2100.0;
  for(i=0;i<51;i++)
  {  l=i/10;
     x[i]=xx;  
     xx+=dx[l]; 
     //System.out.println("i="+i+"dx="+dx[l]+"x="+x[i]); 
  } 
  for(k=1;k<12;k++)
  {
  for(j=0;j<n;j++)
  {
  error[k-1][j+1][0]=eod[j];    
  for(i=0;i<51;i++)
  {
  error[k-1][0][i+1]=x[i];  
  error[k-1][j+1][i+1]=(a[k][j][i]-a[0][j][i])/a[0][j][i]*100;
  }//end of i Re
  }//end of j eod
  }//end of k(equation)
  return error;
}
public static void create_data() throws IOException 
{
 PrintWriter cfout=new PrintWriter(new BufferedWriter(new FileWriter("a.txt"))); 
	double eod[]={1e-8,5e-8,1e-7,5e-7,1e-6,5e-6,1e-5,5e-5,1e-4,5e-4,1e-3,5e-3,1e-2,5e-2};
int n2=eod.length;
  double x[]=new double[2001];
  double dx1=(10000.0-2100)/500.0;
  double dx2=(100000.0-10000.0)/500.0;
  double dx3=(1e6-1e5)/500.0;
  double dx4=(1e7-1e6)/500.0;
  double dx5=(1e7-1e6)/500.0;
  double dx[]={dx1,dx2,dx3,dx4,dx5,dx5};
  double xx=2100.0;
  int l=0;
  for(int i=0;i<2001;i++)
  {  l=i/500;
     x[i]=xx;  
     xx+=dx[l]; 
     //System.out.println("i="+i+"dx="+dx[l]+"x="+x[i]); 
  } 
int n3=x.length;
double a[]=new double[n3];
int i,j;   
  for(j=1;j<(n2-1);j++)
  {   
  for(i=0;i<2001;i++)
  {  
  cfout.println(x[i]+" "+eod[j]+" "+f_colebrook(x[i],eod[j]));
  }//end of i Re
  }//end of j eod
  cfout.close(); 
}

public static void plot()
{
double e[][][]=error();
int n1=e.length;
int n2=e[0].length;
int n3=e[0][0].length;
Plot pp[]=new Plot[n1];
double a[][]=new double[n2][n3];
double eod[]=new double[n2];
double Re[]=new double[n3];
String s[]={"Haaland","Moody","Wood","Churchill","Chen","Swamee","Zigrand","Serghides","Goudor","Romeo"}; 

  double xx=0;
  double x[]=new double[51];
  double dx1=(10000.0-2100)/10.0;
  double dx2=(100000.0-10000.0)/10.0;
  double dx3=(1e6-1e5)/10.0;
  double dx4=(1e7-1e6)/10.0;
  double dx5=(1e7-1e6)/10.0;
  double dx[]={dx1,dx2,dx3,dx4,dx5,dx5};
  xx=2100.0;
  int l=0;
  for(int i=0;i<51;i++)
  {  l=i/10;
     x[i]=xx;  
     xx+=dx[l]; 
     //System.out.println("i="+i+"dx="+dx[l]+"x="+x[i]); 
  } 
int i,j,k;
  for(k=1;k<11;k++)
  {
  for(j=0;j<(n2-1);j++)
  {
  eod[j]=e[k-1][j+1][0];    
  for(i=0;i<51;i++)
  {  
  a[j][i]=e[k][j+1][i+1];
  }//end of i Re
  if(j==0) 
  {pp[k]=new Plot(x,a[j]);
   pp[k].setPlabel(s[k]+" denklemi");
   pp[k].setXlabel("Re");
   pp[k].setYlabel("f_hata %");
  }
  else pp[k].addData(x,a[j]);
  System.out.println(Matrix.toString(x));
  }//end of j eod
  pp[k].plot();
  }//end of k(equation)
}

public static void plot_f()
{
double eod[]={1e-8,5e-8,1e-7,5e-7,1e-6,5e-6,1e-5,5e-5,1e-4,5e-4,1e-3,5e-3,1e-2,5e-2};
int n2=eod.length;
  double x[]=new double[2001];
  double dx1=(10000.0-2100)/10.0;
  double dx2=(100000.0-10000.0)/10.0;
  double dx3=(1e6-1e5)/10.0;
  double dx4=(1e7-1e6)/10.0;
  double dx5=(1e7-1e6)/10.0;
  double dx[]={dx1,dx2,dx3,dx4,dx5,dx5};
  double xx=2100.0;
  int l=0;
  for(int i=0;i<201;i++)
  {  l=i/200;
     x[i]=xx;  
     xx+=dx[l]; 
     //System.out.println("i="+i+"dx="+dx[l]+"x="+x[i]); 
  } 
int n3=x.length;
double a[]=new double[n3];
int i,j;   
  for(i=0;i<2001;i++)
  {  
  a[i]=f_colebrook(x[i],eod[0]);
  }//end of i Re
   Plot pp=new Plot(x,a);
   pp.setPlabel("Moody diagramę");
   pp.setXlabel("Re");
   pp.setYlabel("f");
   pp.setXlogScaleOn();
   pp.setYlogScaleOn();
  for(j=1;j<(n2-1);j++)
  {   
  for(i=0;i<2001;i++)
  {  
  a[i]=f_colebrook(x[i],eod[j]);
  }//end of i Re
  pp.addData(x,a);
  }//end of j eod
  pp.plot();
}

public static void main(String arg[]) throws IOException 
{
	double Re=107163.52439399;
	double eod=0.00184000;
    double f=f_colebrook(Re,eod);
    System.out.println("f : "+f);	
    /*
String s[]={"Haaland","Moody","Wood","Churchill","Chen","Swamee","Zigrand","Serghides","Goudor","Romeo","Turhan"}; 
double e[][][]=error();
for(int i=0;i<12;i++)
{
System.out.println("Equation : "+s[i]);
System.out.println(Matrix.toString(e[i]));
System.out.println("=========================");

}

//create_data();
*/
} 
}