//======================================================
// Numerical Analysis package in java
// Adaptive Gauss-Lobatto integral
// Referans : Adapted Quadrature - Revisited
// Walter Gander - Walter Gautchi
// ftp.inf.ethz.ch
// Dr. Turhan Coban
// =====================================================

import java.io.*;

class f1 extends f_x
{ double func(double x)
  {return 4.0/Math.PI*Math.sqrt(1.0-x*x);}
}

class integ15
{
public static double adaptive_Lobatto_integral(f_x ff,double a,double b)
{double tol=Double.MIN_VALUE;
 return adaptive_Lobatto_integral(ff,a,b,tol);
}

public static double adaptive_Lobatto_integral(f_x ff,double a,double b,double tol)
{   int i;
	double m=(a+b)/2.0; double h=(b-a)/2.0;
    double alpha=Math.sqrt(2.0/3.0); 
    double beta=1.0/Math.sqrt(5.0);
    double x1,x2,x3;
    double eps =Double.MIN_VALUE;
    //System.out.println("eps="+eps);  
    x1=0.94288241569547971905635175843185720232; 
    x2=0.64185334234578130578123554132903188354;
    x3=0.23638319966214988028222377349205292599;
    double A,B,C,D,E,F,G;
    A=0.015827191973480183087169986733305510591;
    B=0.094273840218850045531282505077108171960;
    C=0.15507198733658539625363597980210298680;
    D=0.18882157396018245442000533937297167125;
    E=0.19977340522685852679206802206648840246;
    F=0.22492646533333952701601768799639508076;
    G=0.24261107190140773379964095790325635233;
	double x[]={a,(m-x1*h),(m-alpha*h),(m-x2*h),(m-beta*h),(m-x3*h),m,(m+x3*h),(m+beta*h),(m+x2*h),(m+alpha*h),(m+x1*h),b};
	int n=x.length;
	double y[]=new double[n];
	for(i=0;i<n;i++) y[i]=ff.func(x[i]);
	double fa=y[0]; 
	double fb=y[12];
    double i2=(h/6.0)*(y[0]+y[12]+5.0*(y[4]+y[8]));
    double i1=(h/1470.0)*(77.0*(y[0]+y[12])+432.0*(y[2]+y[10])+
     625.0*(y[4]+y[8])+672.0*y[6]);
    double is=h*(A*(y[0]+y[12])+B*(y[1]+y[11])+C*(y[2]+y[10])+D*(y[3]+y[9])+E*(y[4]+y[8])+F*(y[5]+y[7])+G*y[6]);    
    int s=(int)(is/Math.abs(is)); 
    if(s==0) s=1;
    double erri1=Math.abs(i1-is);
    double erri2=Math.abs(i2-is);
    double R; 
    R=erri1/erri2;
    if(R>0 && R<1) tol=tol/R;
    //System.out.println("is="+is+"s="+s+"tol="+tol+"eps="+eps);
    if(is==0) is=b-a;
  double Q=0;  
  Q=Lobatto(ff,a,b,fa,fb,is,tol);
  return Q;	
}
	
public static double Lobatto(f_x ff,double a,double b,double fa,double fb,double is,double tol)
{
  double h=(b-a)/2; 
  double m=(a+b)/2;
  double alpha=Math.sqrt(2.0/3.0); 
  double beta=1.0/Math.sqrt(5.0);
  double mll=m-alpha*h; 
  double ml=m-beta*h; 
  double mr=m+beta*h; 
  double mrr=m+alpha*h;
  double x[]={mll,ml,m,mr,mrr};
  int n=x.length;
  double y[]=new double[n];
  for(int i=0;i<n;i++) y[i]=ff.func(x[i]);
  double fmll=y[0]; 
  double fml=y[1]; 
  double fm=y[2]; 
  double fmr=y[3]; 
  double fmrr=y[4];
  double Q;
  double i2=(h/6)*(fa+fb+5*(fml+fmr));
  double i1=(h/1470)*(77*(fa+fb)+432*(fmll+fmrr)+625*(fml+fmr)+672*fm);
  double di=Math.abs(i1-i2);
  if (((is + di) == is) || (mll <= a) || (b<=mrr) || di<tol)
  {
   if (((m <= a) || (b<=m)))
   { System.out.println("istenilen toleransa ulaşılamadı delta i = "+di);}
   Q=i1;
  }  
  else
  {
    Q=Lobatto(ff,a,mll,fa,fmll,is,tol) + Lobatto(ff,mll,ml,fmll,fml,is,tol)+
      Lobatto(ff,ml,m,fml,fm,is,tol)   + Lobatto(ff,m,mr,fm,fmr,is,tol)+
      Lobatto(ff,mr,mrr,fmr,fmrr,is,tol)+Lobatto(ff,mrr,b,fmrr,fb,is,tol);
  }
  return Q;
}

  public static void main(String args[]) throws IOException
  {
  f1 b1=new f1();
  System.out.println("PI="+Math.PI+" Ayarlanabilir Gauss-Lobatto integrali f1 : "+adaptive_Lobatto_integral(b1,0.0,1.0));
  }
}
