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

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

class f1 extends f_x
{
   double func(double x)
  {return 100.0/(x*x)*Math.sin(10.0/x);
   //return 1.0/(1+x*x);
  }
}

class integ10B
{
public static double adaptive_newton_cotes_5_integral(f_x ff,double a,double b)
{return adaptive_newton_cotes_5_integral(ff,a,b,1.0e-20,30);}

public static double adaptive_newton_cotes_5_integral(f_x ff,double a,double b,double eps)
{return adaptive_newton_cotes_5_integral(ff,a,b,eps,30);}

public static double adaptive_newton_cotes_5_integral(f_x ff,double a,double b,int MAX_ITERATION)
{return adaptive_newton_cotes_5_integral(ff,a,b,1.0e-12,MAX_ITERATION);}

public static double adaptive_newton_cotes_5_integral(f_x ff,double a,double b,double eps,int MAX_ITERATION)
{   int i;
    double h=(b-a);
    double h1=h/4.0;
    double h2=h/8.0;
	double x[]=new double[9];
    double f[]=new double[9];
    double w1[]=new double[5];
    double z1[]=new double[5]; 
    int j;
    for(i=0;i<9;i+=2) {j=i/2;x[i]=a+i*h2;f[i]=ff.func(x[i]);w1[j]=x[i];z1[j]=f[i];}
	double sum =2.0*h1/45.0*(7.0*f[0]+32.0*f[2]+12.0*f[4]+32.0*f[6]+7.0*f[8]); 
	System.out.println("int1 = "+sum);
	double Q=0;
	Plot pp=new Plot(ff,a,b,100);
	pp.setPlotType(0,0);// 	 
	pp.addData(w1,z1,28);
	int k=0;
    Q=newton_cotes_5(ff,x,f,sum,eps,k,MAX_ITERATION,Q,pp);
    pp.plot();
    return Q;	
}
	
public static double newton_cotes_5(f_x ff,double x[],double f[],double sum,double eps,int k,int MAX_ITERATION,double Q,Plot pp)
{
        
	    double h=(x[8]-x[0]);
        double h1=h/4.0;
        double h2=h/8.0;        
        int p,i;
        double x1[]=new double[9];
        double f1[]=new double[9];
        double x2[]=new double[9];
        double f2[]=new double[9];
        double w1[]=new double[4];
        double z1[]=new double[4]; 
        int j;
        for(i=0;i<9;i++) {x[i]=x[0]+i*h2;f[i]=ff.func(x[i]);}
        for(j=0;i<4;i++) {w1[j]=x[2*j];z1[j]=f[2*j];}           
        for(i=0;i<9;i++) {System.out.println(i+" "+x[i]+" "+f[i]);}
        {System.out.println("==================");}
        double sum1=(2.0*h2/45.0)*(7.0*f[0]+32.0*f[1]+12.0*f[2]+32.0*f[3]+7.0*f[4]);
        double sum2=(2.0*h2/45.0)*(7.0*f[4]+32.0*f[5]+12.0*f[6]+32.0*f[7]+7.0*f[8]);     
        double sum3=sum1+sum2;
        x1[0]=x[0];f1[0]=f[0];
        x1[2]=x[1];f1[2]=f[1];
        x1[4]=x[2];f1[4]=f[2];
        x1[6]=x[3];f1[6]=f[3];
        x1[8]=x[4];f1[8]=f[4];        
        x2[0]=x[4];f2[0]=f[4];
        x2[2]=x[5];f2[2]=f[5];
        x2[4]=x[6];f2[4]=f[6];
        x2[6]=x[7];f2[2]=f[7];
        x2[8]=x[8];f2[4]=f[8]; 
        //System.out.println("sum="+sum);
        //if(1.0/63.0*Math.abs(sum3-sum)<eps || k>=MAX_ITERATION) {Q+=sum3;}
        if( k>=MAX_ITERATION) {Q+=sum3;}
        else
        {         
        k++;  
        System.out.println("giriž k = "+k+"Q="+Q);  
        Q+=newton_cotes_5(ff,x1,f1,sum1,eps,k,MAX_ITERATION,Q,pp)+
           newton_cotes_5(ff,x2,f2,sum2,eps,k,MAX_ITERATION,Q,pp);        
        pp.addData(w1,z1,28);
        }
        //pp.setPlotType(k,28);//              
        //System.out.println("k = "+k+"Q="+Q+"int2 = "+(sum1+sum2)+"sum1="+sum1+"sum2="+sum2);
        return Q;
}

  public static void main(String args[]) throws IOException
  {
  f1 b1=new f1();
  System.out.println("integral of class         f1 : "+adaptive_newton_cotes_5_integral(b1,1.0,3.0,1.0e-10,5));
  }
}
