// test_aquad.java 
public class test_aquad
{
  public test_aquad()
  {
    double xmin = pass.a;
    double xmax = pass.b;
    double eps  = 0.001;
    double area, err, exact;

    System.out.println("test_aquad.java testing aquad.java x^2, later 1/x");
    area = aquad(xmin, xmax, eps);
    // exact = Math.log(xmax)-Math.log(xmin); // for 1/x
	exact = 8.0/3.0-xmin*xmin*xmin/3.0; // for x^2
    err = area-exact;
    System.out.println("xmin="+xmin+", xmax="+xmax+", area="+area+
                       ", exact="+exact+", err="+err+"\n");

    xmin = 0.01;
    area = aquad(xmin, xmax, eps); // f(x) passed below
    // exact = Math.log(xmax)-Math.log(xmin); // for 1/x
	exact = 8.0/3.0-xmin*xmin*xmin/3.0; // for x^2
    err = area-exact;
    System.out.println("xmin="+xmin+", xmax="+xmax+", area="+area+
                       ", exact="+exact+", err="+err+"\n");

    xmin = 0.001;
    area = aquad(xmin, xmax, eps);
    // exact = Math.log(xmax)-Math.log(xmin); // for 1/x
	exact = 8.0/3.0-xmin*xmin*xmin/3.0; // for x^2
    err = area-exact;
    System.out.println("xmin="+xmin+", xmax="+xmax+", area="+area+
                       ", exact="+exact+", err="+err+"\n");

    xmin = 0.0001;
    area = aquad(xmin, xmax, eps);
    // exact = Math.log(xmax)-Math.log(xmin); // for 1/x
	exact = 8.0/3.0-xmin*xmin*xmin/3.0; // for x^2
    err = area-exact;
    System.out.println("xmin="+xmin+", xmax="+xmax+", area="+area+
                       ", exact="+exact+", err="+err+"\n");

    xmin = 0.00001;
    area = aquad(xmin, xmax, eps);
    // exact = Math.log(xmax)-Math.log(xmin); // for 1/x
	exact = 8.0/3.0-xmin*xmin*xmin/3.0; // for x^2
    err = area-exact;
    System.out.println("xmin="+xmin+", xmax="+xmax+", area="+area+
                       ", exact="+exact+", err="+err+"\n");

    //xmin = 0.0;
    //xmax = 1.0;
    //eps = 0.001;
    //area = aquad(f1, xmin, xmax, eps);
    //exact = 29.85832540;
    //err = area - exact;  
    //System.out.println("xmin="+xmin+", xmax="+xmax+", area="+area+
    //                   ", exact="+exact+", err="+err+"\n");

    System.out.println("test_aquad.java finished");
  }

  // aquad   adaptive quadrature  numerical integration test 
  //         for a desired accuracy on irregular functions.   
  // define the function to be integrated as:                 
  //  double f(double x)   in file pass_funct.java                                   
  //  {                                                       
  //    // compute value                                      
  //    return value;                                         
  //  }                                                       
  //                                                          
  // integrate from xmin to xmax                              
  // approximate desired absolute error in all intervals, eps 
  // accuracy absolutely not guaranteed                       

  double aquad(double xmin, double xmax, double eps)
               // function f(x) comes from pass.f(x)
  {
    double area, temp, part, h;
    double err, minint, maxerr;
    int nmax = 2000;
    double sbin[] = new double[2000];
    double ebin[] = new double[2000];
    double abin[] = new double[2000];
    int    fbin[] = new int[2000];
    int i, j, k, n, nn, done, nf;
    int kmax = 20;

    n=32; // initial number of bins 
    h = (xmax-xmin)/(double)n;
    for(i=0; i<n; i++)
    {
      sbin[i] = xmin+i*h;
      ebin[i] = xmin+(i+1)*h;
      fbin[i] = 0;
    }
    minint = xmax-xmin;
	maxerr = 0.0;
	area = 0.0;
    k = 0;
    done = 0;
    nf = 0; // number of function evaluations 
    nn = n; // next available bin 
    while(done==0)
    {
      done = 1;
      k++;
      if(k>=kmax) break; // quit if more than kmax subdivisions 
      area = 0.0;
      maxerr = 0.0;
      for(i=0; i<n; i++)
      {
        if(fbin[i]==1) // this interval finished 
        {
          area = area + abin[i]; // accumulate total area each pass 
          continue;
        }
        temp = pass.f((sbin[i]+ebin[i])/2.0); // two integration methods 
        part = pass.f((3.0*sbin[i]+ebin[i])/4.0);
        part = part + pass.f((sbin[i]+3.0*ebin[i])/4.0);
        nf = nf + 3;
        abin[i] = (part+2.0*temp)*(ebin[i]-sbin[i])/4.0;
        area = area + abin[i];
        err = Math.abs(temp-part/2.0);
        if(err > maxerr) maxerr = err;
        if(err*1.414 < eps) // heuristic 
        {
          fbin[i] = 1; // this interval finished 
        }
        else
        {
          done = 0; // must keep dividing 
          if(nn>=nmax) // quit, out of space 
          {
            done = 1;
            for(j=i+1;  j<n; j++) area = area + abin[j];
            break;     // result not correct 
	  }
          else // divide interval into two halves 
          {
            if((sbin[i]+ebin[i])/2.0<minint) minint=(sbin[i]+ebin[i])/2.0;
            sbin[nn] = (sbin[i]+ebin[i])/2.0;
            ebin[nn] = ebin[i];
            fbin[nn] = 0;
            ebin[i] = sbin[nn];         
            nn++;
	  }
        }
      } // end for i 
      n = nn;
    } // end while   
    System.out.println(nn+" intervals, "+nf+" funeval, "+k+" divides, small="+
                       minint+", maxerr="+maxerr+" \n");

    return area;
  } // end aquad


  //double f(double x){return 1.0/x;} comes from pass.f(x)
  //                   compile into pass_funct.java
  //double f1(double x){return 1.0/((x-0.3)*(x-0.3)+0.01) + 1.0/((x-0.9)*(x-0.9)+0.04) -6.0;}

  public static pass_funct pass=new pass_funct();
  
  public static void main (String[] args)
  {
    new test_aquad();
  }
} // end test_aquad.java