// ======================================================
// Numerical Analysis package in java
// example to show utilisation of integration (integral)
// and differentiation (derivative) functions
// Adaptive Gauss-legendre integration 
// Dr. Turhan Coban
// ======================================================
import java.io.*;


class f4 extends f_x
{ double a[][];
  public double T;
  static double sigma=5.670503208914327e-8;
  static double C1=3.74177489e8;
  static double C2=1.438769e4;
    
  public f4(double x[][],double Ti)
  {a=x;T=Ti;}
    
  public f4(String name,double Ti)
  {double x[][]=Text.readDoubleT(name);a=x;T=Ti;
   Text.print(x);
  }
  
  public f4(double xlimiti,double Ti)
  { 
    double x[][]=Text.readDoubleT("a.txt");
    a=x;
    T=Ti;
    Text.printT(x);
  }
  
  public double func(double x)
  {
   double y=C2/(x*T);  
   double eps=0;
   eps=movingLSTSQR((x*T),a[0],a[1],2,2);
   //System.out.println("x="+x+"xT="+x*T+"eps="+eps);
   double integ=eps/sigma/(T*T*T*T)*C1/(x*x*x*x*x)/(Math.exp(y)-1.0);
   return integ;
  }
  
//polynomial moving least square 
  public static double[] pivotlugauss(double a[][],double b[])
  { //gauss elimination with pivoting
  int n=b.length;  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gauss eleme
  int i,j,k,p,ii,jj;
  for(k=0;k<(n-1);k++)
  { //pivotlama
	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;
	}
	//gauss elemeyi çözme  
    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];
	}
  }
  //geriye doğru yerine koyma
  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 movingLSTSQR(double x,double xi[],double yi[],int n,int kk)
{
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];
double w[]=new double[l];
for(i=0;i<l;i++) {if(x==xi[i]) return yi[i];}
for(i=0;i<l;i++) {w[i]=1.0/Math.pow(Math.abs(x-xi[i]),kk);}
for(i=0;i<np1;i++)
  {  for(j=0;j<np1;j++)
     {if(i==0 && j==0) A[i][j]=l;
     else  for(k=0;k<l;k++) A[i][j] += w[k]*Math.pow(xi[k],(i+j));
     }
     for(k=0;k<l;k++) { if(i==0)  B[i]+= w[k]*yi[k];
                        else      B[i] += w[k]*Math.pow(xi[k],i)*yi[k];}
  } 
X=pivotlugauss(A,B);
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;
return funcEKK(X,x);
}

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

}


class emissivity4
{
public static double[][] gauss_legendre_coefficients(double x1,double x2,int n)
{
	// calculates legendre gauss-coefficients as coefficients of the integral
	// for n terms
	double EPS=3.0e-5;
	int m,j,i;
	double z1,z,xm,xl,pp,p3,p2,p1;
                 double a[][]=new double[2][n];
	m=(n+1)/2;
	xm=0.5*(x2+x1);
	xl=0.5*(x2-x1);
	for (i=1;i<=m;i++)  {
		z=Math.cos(Math.PI*((i-0.25)/(n+0.5)));
		do {
			p1=1.0;
			p2=0.0;
			for (j=1;j<=n;j++) {
				p3=p2;
				p2=p1;
				p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
			}
			pp=n*(z*p1-p2)/(z*z-1.0);
			z1=z;
			z=z1-p1/pp;
		} while (Math.abs(z-z1) > EPS);
		a[0][i-1]=xm-xl*z;
		a[0][n-i]=xm+xl*z;
		a[1][i-1]=2.0*xl/((1.0-z*z)*pp*pp);
		a[1][n-i]=a[1][i-1];
			}
return a;
}

public static double gauss_legendre_integral(f_x f,double x1,double x2,int n)
{
//n : number of integral coefficients
// this routine first generates gauss legendre coefficients
// for [x1,x2] band 
// then calculates gauss legendre integral
double a[][]=new double[2][n];
a=gauss_legendre_coefficients(x1,x2,n);
//Text.print(Text.T(a));
//System.out.println(Matrix.toString(a));
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*f.func(a[0][i]);}
return z;
}

public static double adaptive_gauss_legendre_integral(f_x f,double x1,double x2)
{return adaptive_gauss_legendre_integral(f,x1,x2,1.0e-12,100);}

public static double adaptive_gauss_legendre_integral(f_x f,double x1,double x2,double eps)
{return adaptive_gauss_legendre_integral(f,x1,x2,eps,500);}

public static double adaptive_gauss_legendre_integral(f_x f,double x1,double x2,double eps,double MAXITER)
{
int n=1;
double ans1=gauss_legendre_integral(f,x1,x2,n);
double ans2;
do
{ans2=ans1;n++;ans1=gauss_legendre_integral(f,x1,x2,n);}
while(Math.abs(ans2-ans1)>eps && n<MAXITER);
if(n==MAXITER) System.out.println("istenilen hata seviyesi tutturulamadı n="+n );
return ans2;
}

  public static double F(double x)
  {
  double C2=1.438769e4;
  double y=C2/x;
  double total=0;
  for(int n=1;n<40;n++)
  {total+=Math.exp(-n*y)/n*(y*y*y+3.0*y*y/n+6.0*y/(n*n)+6.0/(n*n*n));}  
  double pi=Math.PI;
  total=15.0/(pi*pi*pi*pi)*total;
  return total;
  }
  
 
  	  
  public static void main(String args[]) throws IOException
  {
  //Gauss-Legendre integrali
  double T=800.0; //K
  f4 b1=new f4("a.txt",T);
  double e=adaptive_gauss_legendre_integral(b1,0.28,40.0);
  System.out.println("emmissivity  : "+e);
}
}
