import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
import java.io.*;

class f2 extends f_xj
{//adaptation function
fi_xi ff1;
public f2(fi_xi ffi) {ff1=ffi;}	

public double func(double x[]) 
{	
double ff=0.0;
double fa[]=ff1.func(x);
for(int i=0;i<fa.length;i++)  ff+=fa[i]*fa[i];
return ff; 
}
}

class fa extends fi_xi
{ // Nonlinear system of equation
double y[];
double w;
fm2 b3;
double a[][];
public fa(double Pr,double yi[],double wi)
{ b3=new fm2(Pr);
  y=new double[yi.length];
  for(int i=0;i<y.length;i++) y[i]=yi[i];w=wi;
}

public double[] func(double x[])
{		
// shooting method
y[2]=x[0];
y[4]=x[1];
//System.out.println(y[0]+" "+y[1]+" "+y[2]+" "+y[3]+" "+y[4]+" ");
a=RK6(b3,0.0,w,y,10000);
//System.out.println(a[0][1000]+" "+a[1][1000]+" "+a[2][1000]+" "+a[3][1000]+" "+a[4][1000]+" "+a[5][1000]+" ");
double ff[]=new double[2];
ff[0]=a[2][10000]-0.0;
ff[1]=a[4][10000]-0.0;
return ff;
} 
 //differenetial equation solution
public static double[][] RK6(f_xi fp,double x0,double xn,double f0[],int N)
{
//6th order Runge Kutta Method
//fp  : given set of derivative functions dfj/dxi(fj,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
int M=f0.length;
double fi[][];
fi=new double[M][N+1];
double xi[]=new double[M+1];
double k[]=new double[6];
int i,j;
double x;
for(j=0;j<M;j++)
  {
  fi[j][0]=f0[j];
  xi[j+1]=f0[j];
  }
for(x=x0,i=0;i<N;x+=h,i++)
  {
  for(j=1;j<=M;j++)
    {
    xi[0]=x;
    xi[j]=fi[j-1][i];
    k[0]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/2;
    k[1]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/4.0+k[1]/4.0;
    k[2]=h*fp.func(xi,j-1);
    xi[0]=x+h;
    xi[j]=fi[j-1][i]-k[1]+2.0*k[2];
    k[3]=h*fp.func(xi,j-1);
    xi[0]=x+2.0/3.0*h;
    xi[j]=fi[j-1][i]+7.0/27.0*k[0]+10.0/27.0*k[1]+1.0/27.0*k[3];
    k[4]=h*fp.func(xi,j-1);
    xi[0]=x+1.0/5.0*h;
    xi[j]=fi[j-1][i]+28.0/625.0*k[0]-1.0/5.0*k[1]+546.0/625.0*k[2]+54.0/625.0*k[3]-378/625.0*k[4];
    k[5]=h*fp.func(xi,j-1);
    fi[j-1][i+1]=fi[j-1][i]+k[0]/24.0+5.0*k[3]/48.0+27.0*k[4]/56.0+125.0*k[5]/336.0;
    xi[j]=fi[j-1][i];
    }
  }
double a[][]=new double[M+1][N+1];
for(x=x0,i=0;i<=N;x+=h,i++)
{
  a[0][i]=x;
  for(j=1;j<=M;j++)
  {
  a[j][i]=fi[j-1][i];
  }
}
return a;
}
}

//Blassius flat plate differential equation
class fm2 extends f_xi
{ double Pr; 
  fm2(double Pri) {Pr=Pri;}
// multivariable function
 double func(double x[],int x_ref)
 {
   //Natural convection similarity solution
   double a=0;
   if(x_ref==4) a= -3.0*Pr*x[1]*x[5];
   if(x_ref==3) a= x[5];
   if(x_ref==2) a= -3.0*x[1]*x[3]+2.0*x[2]*x[2]-x[4];
   if(x_ref==1) a= x[3];
   if(x_ref==0) a= x[2];
   return a;
 }
}
class fb extends f_xr
{
	double func(double Pr,int i)
	{	
     double xx=0.0;
      if(i==0) xx=Math.pow(Pr,0.25);
      else if(i==1) xx=Math.sqrt(Pr);
      else if(i==2) xx=Math.pow(Pr,0.75);
      else if(i==3) xx=Pr;
      else if(i==4) xx=Math.pow(Pr,1.25);
      else if(i==5) xx=Math.pow(Pr,1.5);
      else if(i==6) xx=Math.pow(Pr,1.75);
   return xx;
	}
}    

public class natural_convection_similarity
{
	 public static double[] inputdata(String s)
     {
     String s1=JOptionPane.showInputDialog(s);        
     StringTokenizer token=new StringTokenizer(s1);
     int n=token.countTokens()-1;
     int m=n+1;
     double a[]=new double[m];
     int j=0;           
     while(token.hasMoreTokens())
     {
     Double ax=new Double(token.nextToken());
     a[j++]=ax.doubleValue();
     }
     return a; 
     }	
public static String toString(double x[])
{ String s="";
int n=x.length;
  for(int k=0;k<n;k++)
  {s+=x[k]+" ";}
  return s;
}
public static double[] nelder(f_xj fnelder,double a[],double da[],int maxiteration,double tolerance,int printlist)
{
	int i,j;
double x[][]=new double[a.length+1][a.length];
double p[][]=new double[a.length+1][a.length+1];
for(i=0;i<x.length;i++)
   {for(j=0;j<x[0].length;j++)
     {if(i==j){x[i][j]=a[j]+da[j];p[i][j]=x[i][j];}
      else    {x[i][j]=a[j];p[i][j]=x[i][j]; }
     }
     p[i][j] = fnelder.func(p[i]);
   }


   // Inlet variable definitions
   // fnelder : abstract multivariable function f(x)
   // x : independent variable set of n+1 simplex elements
   // maxiteration : maximum iteration number
   // tolerance : 
         int NDIMS = x.length-1;
         int NPTS = x.length;
         int FUNC = NDIMS;
         int ncalls = 0;
            ////// construct the starting simplex //////////////////
            //double p[][]=new double[NPTS][NPTS]; // [row][col] = [whichvx][coord,FUNC]
            double z[]=new double[NDIMS];
            double best = 1E99;
            //////////////// calculate the first function values for the simplex ////////////////
 
            int iter=0;
            
            for (iter=1; iter<maxiteration; iter++)
            {
	            ///////////  define lo, nhi, hi (low high next_to_high //////////////
                int  ilo=0, ihi=0, inhi = -1; // -1 means missing
                double flo = p[0][FUNC];
                double fhi = flo;
                double pavg,sterr;
                for (i=1; i<NPTS; i++)
                {
                   if (p[i][FUNC] < flo)
                     {flo=p[i][FUNC]; ilo=i;}
                   if (p[i][FUNC] > fhi)
                     {fhi=p[i][FUNC]; ihi=i;}
                }
                double fnhi = flo;
                inhi = ilo;
                for (i=0; i<NPTS; i++)
                  if ((i != ihi) && (p[i][FUNC] > fnhi))
                    {fnhi=p[i][FUNC]; inhi=i;}
                ////////// exit criteria //////////////
                if ((iter % 4*NDIMS) == 0)
                {
                    // calculate the avarage (including maximum value)
                    pavg=0;
                    for(i=0;i<NPTS;i++)
                      pavg+=p[i][FUNC];
                    pavg/=NPTS;
                    double tot=0;
                    if(printlist!=0)
                    {  System.out.print(iter);
                       for (j=0; j<=NDIMS; j++)
                          { System.out.print(p[ilo][j]+" ");}
                       System.out.println("");
                    }
                    for(i=0;i<NPTS;i++)
                    { tot=(p[i][FUNC]-pavg)*(p[i][FUNC]-pavg);}
                    sterr=Math.sqrt(tot/NPTS);
                    //if(sterr < tolerance)
                      { for (j=0; j<NDIMS; j++)
                           { z[j]=p[ilo][j];}
                        //break;
                      }
                    best = p[ilo][FUNC];
                }

                ///// calculate avarage without maximum value //////

                double ave[] = new double[NDIMS];
                for (j=0; j<NDIMS; j++)
                  ave[j] = 0;
                for (i=0; i<NPTS; i++)
                  if (i != ihi)
                    for (j=0; j<NDIMS; j++)
                       ave[j] += p[i][j];
                for (j=0; j<NDIMS; j++)
                   ave[j] /= (NPTS-1);


                ///////// reflect ////////////////

                double r[] = new double[NDIMS];
                for (j=0; j<NDIMS; j++)
                  r[j] = 2*ave[j] - p[ihi][j];
                double fr = fnelder.func(r);

                if ((flo <= fr) && (fr < fnhi))  // in zone: accept
                {
                    for (j=0; j<NDIMS; j++)
                      p[ihi][j] = r[j];
                    p[ihi][FUNC] = fr;
                    continue;
                }

                if (fr < flo)  //// expand
                {
                    double e[] = new double[NDIMS];
                    for (j=0; j<NDIMS; j++)
                      e[j] = 3*ave[j] - 2*p[ihi][j];
                    double fe = fnelder.func(e);
                    if (fe < fr)
                    {
                       for (j=0; j<NDIMS; j++)
                         p[ihi][j] = e[j];
                       p[ihi][FUNC] = fe;
                       continue;
                    }
                    else
                    {
                       for (j=0; j<NDIMS; j++)
                         p[ihi][j] = r[j];
                       p[ihi][FUNC] = fr;
                       continue;
                    }
                }

                ///////////// shrink:

                if (fr < fhi)  
                {
                    double c[] = new double[NDIMS];
                    for (j=0; j<NDIMS; j++)
                      c[j] = 1.5*ave[j] - 0.5*p[ihi][j];
                    double fc = fnelder.func(c);
                    if (fc <= fr)
                    {
                        for (j=0; j<NDIMS; j++)
                          p[ihi][j] = c[j];
                        p[ihi][FUNC] = fc;
                        continue;
                    }
                    else   /////// daralt
                    {
                        for (i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (j=0; j<NDIMS; j++)
                                p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                              p[i][FUNC] = fnelder.func(p[i]);
                          }
                        continue;
                    }
                }

                if (fr >= fhi)   /// 
                {
                    double cc[] = new double[NDIMS];
                    for (j=0; j<NDIMS; j++)
                      cc[j] = 0.5*ave[j] + 0.5*p[ihi][j];
                    double fcc = fnelder.func(cc);
                    if (fcc < fhi)
                    {
                        for (j=0; j<NDIMS; j++)
                          p[ihi][j] = cc[j];
                        p[ihi][FUNC] = fcc;
                        continue;
                    }
                    else    ///////// 
                    {
                        for (i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (j=0; j<NDIMS; j++)
                                p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                              p[i][FUNC] = fnelder.func(p[i]);
                          }
                    }
                }
            }
       return z;
}
      public static double[] nelder(f_xj fnelder,double a[],double da[],double tolerance)
      {return nelder(fnelder,a,da,500,tolerance,0);}

      public static double[] nelder(f_xj fnelder,double a[],double da[])
      {return nelder(fnelder,a,da,500,1.0e-10,0);}
      
      public static double[] nelder(f_xj fnelder,double a[])
      {
	          double [] da=new double[a.length];
              for(int i=0;i<a.length;i++) da[i]=0.1*a[i];
              return nelder(fnelder,a,da); 
      }	
//general least square curve fitting	
  public static double[] gausswithpartialpivot(double a[][],double b[])
  { //Gauss elimination with partial pivoting
  int n=b.length;	  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gauss elimination
  int i,j,k,p,ii,jj;
  for(k=0;k<(n-1);k++)
  { //partial pivoting
	p=k;
	buyuk=Math.abs(a[k][k]);
	for(ii=k+1;ii<n;ii++)
	{  dummy=Math.abs(a[ii][k]);
	   if(dummy > buyuk) {buyuk=dummy;p=ii;}
	}
	if(p!=k)
	{ for(jj=k;jj<n;jj++)
	  { dummy=a[p][jj];
	    a[p][jj]=a[k][jj];
	    a[k][jj]=dummy;
	  }
	  dummy=b[p];
	  b[p]=b[k];
	  b[k]=dummy;
	}
	//  
    for(i=k+1;i<n;i++)
    {  carpan=a[i][k]/a[k][k];
       a[i][k]=0;
       for(j=k+1;j<n;j++)
       {   a[i][j]-=carpan*a[k][j]; }
           b[i]   =b[i]   -carpan*b[k];
	}
  }
  //backward substitution
  x[n-1]=b[n-1]/a[n-1][n-1];
  for(i=n-2;i>=0;i--)
  {
    toplam=0;
    for(j=i+1;j<n;j++)
    {  toplam+=a[i][j]*x[j];}
    x[i]=(b[i]-toplam)/a[i][i];
  }
  return x;
}	  

  
public static double[] GeneralLeastSquare(double xi[],double yi[],int n)
{
fb f=new fb();
int l=xi.length;
int i,j,k;
int np1=n+1;
double A[][];
A=new double[np1][np1];
double B[];
B=new double[np1];
double X[];
X=new double[np1];
for(i=0;i<n+1;i++)
  {  B[i]=0;
	 for(j=0;j<np1;j++)
     {
	 A[i][j]=0.0;    
     for(k=0;k<l;k++) A[i][j]+=f.func(xi[k],i)*f.func(xi[k],j);
     }
     for(k=0;k<l;k++) B[i]+= f.func(xi[k],i)*yi[k];
  }
//System.out.println("A = \n"+Matrix.toString(A));
//System.out.println("B = \n"+Matrix.toStringT(B));  
X=gausswithpartialpivot(A,B);
//X=B/A;
double max=0;
for(i=0;i<n+1;i++)
 if(Math.abs(X[i]) > max) max = Math.abs(X[i]);
for(i=0;i<n+1;i++)
 if((Math.abs(X[i]/max) > 0) && (Math.abs(X[i]/max) < 1.0e-100)) X[i]=0;
Text.printT(X); 
return X;
}

public static double funcGeneralLeastSquare(double e[],double x)
{
// this function calculates the value of
// least square curve fitting function
fb f=new fb();
int n=e.length;
double ff=0;
if(n!=0.0)
  { 
  for(int i=n-1;i>=0;i--)
   {ff+=e[i]*f.func(x,i);}
  } 
return ff;
}

public static double hata(double x[],double y[],double e[])
{
//calculates absolute square root error of a least square approach
double n=x.length;
 int k;
 double total=0;
 for(k=0;k<n;k++)
 {
 total+=(y[k]-funcGeneralLeastSquare(e,x[k]))*(y[k]-funcGeneralLeastSquare(e,x[k]));
 }
total=Math.sqrt(total);
return total;
}

 public static double[][]  funcGeneralLeastSquare(double xi[],double yi[],int polinomkatsayisi,int aradegersayisi)
 {
	//aradegersayisi:  x--o--o--x--o--o--x zincirinde x deneysel noktalar ise
	// ara dešer sayisi 2 dir
	 int n=xi.length;
	int nn=(n-1)*(aradegersayisi+1)+1;
	double z[][]=new double[2][nn];
    double E[]=GeneralLeastSquare(xi,yi,polinomkatsayisi);	
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {z[0][k]=xi[i];z[1][k]=funcGeneralLeastSquare(E,z[0][k]);k++;
       for(int j=0;j<aradegersayisi;j++) 
       {dx=(xi[i+1]-xi[i])/((double)aradegersayisi+1.0);
       z[0][k]=z[0][k-1]+dx;z[1][k]=funcGeneralLeastSquare(E,z[0][k]);k++;}
       }
       z[0][k]=xi[i];z[1][k]=funcGeneralLeastSquare(E,z[0][k]);
    return z;
 } 


public static void main(String arg[]) throws IOException
{ //initial guess for differential equation set
  double y[]={0,0,0.7,1.0,-0.5};
 //String s="leave a space between initial root guess values: xmin0 xmin1.. ";
 //double x[]=inputdata(s);
 double x[]={0.7,-0.5};//inputdata(s);
 fa ff;
 f2 f; 
 double r1[];
 //selected Prandtl number values
 double Pr[]={0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10,11,12,13,14,15};
 //selected values for eta corresponding to the Pr numbers
 double w[]={16.0,15.0,15.0,15.0,15.0,15.0,15.0,15.0,15.0,10.0,10.0,9.0,8.0,8.0,9.0,7.0,7.0,7.0,6.0,5.0,5.0,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5,4.5};
 //output data for Nusselt's number curve fitting
 double xi[]=new double[Pr.length];
 double yi[]=new double[Pr.length];
 //Initial calculations
 //non linear eqation function
 ff=new fa(Pr[0],y,w[0]);
 //optimisation function
 f=new f2(ff);
 r1=nelder(f,x);
 xi[0]=Pr[0]; 
 yi[0]=-r1[1];
 Plot pp=new Plot();
 Plot pp1=new Plot();
 for(int i=0;i<xi.length;i++)
 {ff=new fa(Pr[i],y,w[i]);
  f=new f2(ff);
  r1=nelder(f,x);
   xi[i]=Pr[i]; 
   yi[i]=-r1[1];
  System.out.println("wall surface Pr="+Pr[i]+" eta="+ff.a[0][0]+" f= "+ff.a[1][0]+" f'= "+ff.a[2][0]+" f''= "+ff.a[3][0]+" T = "+ff.a[4][0]+" T' = "+ff.a[5][0]+" "); 
  System.out.println("boundary layer limit Pr= "+Pr[i]+" eta = "+ff.a[0][10000]+" f = "+ff.a[1][10000]+" f' = "+ff.a[2][10000]+" f'' =  "+ff.a[3][10000]+" T = "+ff.a[4][10000]+" T' = "+ff.a[5][10000]+" ");
   String s1[]={""+'\u03B7',"f","f'","f''",""+'\u03B8',""+'\u03B8'+"'"};
   if(i==0) {Plot pp2=new Plot(ff.a[0],ff.a[2]);pp=pp2;
             Plot pp3=new Plot(ff.a[0],ff.a[4]);pp1=pp3;
            }
   else
   { pp.addData(ff.a[0],ff.a[2]);
     pp1.addData(ff.a[0],ff.a[4]);
   }
 }
  double gPr1[]=new double[Pr.length];
  double gPr2[]=new double[Pr.length];
  double egPr1[]=new double[Pr.length];
  double egPr2[]=new double[Pr.length];
  double b[]=GeneralLeastSquare(xi,yi,6);  
  for(int i=0;i<Pr.length;i++)
  {gPr1[i]=b[0]*Math.pow(Pr[i],0.25)+b[1]*Math.pow(Pr[i],0.5)+b[2]*Math.pow(Pr[i],0.75)+b[3]*Pr[i]+b[4]*Math.pow(Pr[i],1.25)+b[5]*Math.pow(Pr[i],1.5)+b[6]*Math.pow(Pr[i],1.75);
   gPr2[i]=0.75*Math.sqrt(Pr[i])/Math.pow((0.609+1.221*Math.sqrt(Pr[i])+1.238*Pr[i]),0.25);
   egPr1[i]=(gPr1[i]-yi[i])/yi[i]*100;
   egPr2[i]=(gPr2[i]-yi[i])/yi[i]*100;
   }
  
  Text.printT(b);
  double d[][]={xi,yi};
  Text.printT(d,"d");
  System.out.println("f(Pr) ="+b[0]+b[1]+"*Pr^0.25+"+b[2]+"*Pr^0.5+"+b[3]+"*Pr^0.75+"+b[4]+"*Pr +"+b[5]+"*Pr^1.25+"+b[6]+"*Pr^1.5");
 pp.plot();
 pp1.plot();
 Plot pp5=new Plot(xi,yi);
 pp5.addData(xi,gPr1);
 pp5.setColor(1,0,0,255);
 pp5.setPlotType(1,23);
 pp5.addData(xi,gPr2);
  pp5.setColor(2,255,0,0);
 pp5.setPlotType(2,28);
 pp5.plot();
 Plot pp6=new Plot(xi,egPr1);
 pp6.setColor(1,0,0,255);

 pp6.addData(xi,egPr2);
  pp6.setColor(2,255,0,0);
 pp6.plot();
}
}