/*************************************************************************/
/* Program file name: Ex3_9.java                                         */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 15, 2006                                          */
/*                                                                       */
/* (1) This Java program is created for the book, "An Introduction to    */
/*     Computational Physics, 2nd Edition," written by Tao Pang and      */
/*     published by Cambridge University Press on January 19, 2006.      */
/*                                                                       */
/* (2) No warranties, express or implied, are made for this program.     */
/*                                                                       */
/*************************************************************************/

// An example of evaluating an integral with the adaptive
// trapezoid rule.

import java.lang.*;
public class Ex3_9 {
  static final int n = 100;
  public static void main(String argv[]) {
    double del = 1e-6;
    int k = 0;
    double a = 0;
    double b = Math.PI;
    double s = trapezoid(a, b, del, k, n);
    System.out.println ("S = " + s);
  }

// Method to integrate over f(x) with the adaptive
// Simpson rule.
  
  public static double trapezoid(double a, double b,
    double del, int step, int maxstep) {
    double h = b-a;
    double c = (b+a)/2;
    double fa = f(a);
    double fc = f(c);
    double fb = f(b);
    double s0 = h*(fa+fb)/2;
    double s1 = h*(fa+2*fc+fb)/4;
    step++;
    if (step >= maxstep) {
      System.out.println ("Not converged after "
        + step + " recursions");
      return s1;
    }
    else {
      if (Math.abs(s1-s0) < 3*del) return s1;
      else return trapezoid(a, c, del/2, step, maxstep)
        + trapezoid(c, b, del/2, step, maxstep);
    }
  }

// Method to provide the integrand f(x).

  public static double f(double x) {
    double a0 = 5;
    double s = Math.sin(x);
    double c = Math.cos(x);
    double f = (1+a0*(1-c)*(1-c))
      /((1+a0*s*s)*Math.sqrt(1+2*a0*(1-c)));
    return f;
  }
}