//======================================================
// Numerical Analysis package in java
// Adaptive Simpson integral 
// Referans : Adapted Quadrature - Revisited
// Walter Gander - Walter Gautchi
// ftp.inf.ethz.ch
// error control :
// (((is + (I1-I2)) == is) || (m <= a) || (b<=m))
// Dr. Turhan Coban
// ======================================================
import java.io.*;

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


class integ16
{
public static double adaptive_simpson_integral2(f_x ff,double a,double b)
{   int i;
	Plot pp=new Plot(ff,a,b,100);
    double h=(b-a);
    double h2=h/4.0;
    double h3=h/8.0;
    double m=(a+b)/2.0;
    double x[]=new double[9];
    double y[]=new double[9];
    for(i=0;i<9;i++) {x[i]=a+i*h3;y[i]=ff.func(x[i]);}
    double is=h/28350.0*(989.0*y[0]+5888.0*y[1]-928.0*y[2]+10496.0*y[3]-4540.0*y[4]+
           10496.0*y[5]-928.0*y[6]+5888.0*y[7]+989.0*y[8]); 
    double fa = y[0]; 
    double fm = y[4]; 
    double fb = y[8];
    if(is==0) is = b-a;
    double Q = adaptsimstp(ff,a,b,fa,fm,fb,is);
    return Q;	
}
	
public static double adaptsimstp(f_x ff,double a,double b,double fa,double fm,double fb,double is)
{ double m = (a + b)/2; 
  double h = (b - a)/4;
  double x[] = {a + h, b - h};
  int n=x.length;
  double y[] = new double[n];
  for(int i=0;i<n;i++) {y[i]=ff.func(x[i]);}
  double fml = y[0];
  double fmr = y[1];
  double  i1 = h/1.5*(fa+4.0*fm+fb);
  double  i2 = h/3.0*(fa+4.0*(fml+fmr)+2.0*fm+fb);
  double Q=0;
  i1 = (16.0*i2 - i1)/15.0; //romberg integrasyonu ile iyileştirme
  if (((is + (i1-i2)) == is) || (m <= a) || (b<=m))
  { if ((m <= a) || (b<=m))
    { System.out.println("istenilen toleransa ulaşılamadı");}
    Q = i1;
  }  
  else
  {Q = adaptsimstp(ff,a,m,fa,fml,fm,is)+adaptsimstp(ff,m,b,fm,fmr,fb,is);}
  return Q;
}

  public static void main(String args[]) throws IOException
  {
  f1 b1=new f1();
  System.out.println("integral of class         f1 : "+adaptive_simpson_integral2(b1,0.0,1.0));
  }
}
