// ======================================================
// 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 f1 extends f_x
{ public double T;
  static double sigma=5.670503208914327e-8;
  static double C1=3.74177489e8;
  static double C2=1.438769e4;
  f1(double Ti)
  {T=Ti;}  
  public double func(double x)
  { return (1.27917-0.10833*x)/sigma/(T*T*T*T)*C1/(x*x*x*x*x)/(Math.exp(C2/(T*x))-1.0);}
}

class emissivity2
{
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-15;
	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,100);}

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
  f1 b1=new f1(650.0);
  double e=0.9*F(2275)+adaptive_gauss_legendre_integral(b1,3.5,9.5)+0.25*(1-F(6175));
  System.out.println("emmissivity  : "+e);
}
}
