///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Collide.java                                       //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: January 18, 2006                                       //
//                                                                       //
// (1) This Java program is part of 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 calculating the differential cross section
// of classical scattering on the Yukawa potential.

import java.lang.*;
public class Collide {
  static final int n = 10000, m = 20;
  static final double a = 100, e = 1;
  static double b;
  public static void main(String argv[]) {
    int nc = 20, ne = 2;
    double del = 1e-6, db = 0.5, b0 = 0.01, h = 0.01;
    double g1, g2;
    double theta[] = new double[n+1];
    double fi[] = new double[n+1];
    double sig[] = new double[m+1];
    for (int i=0; i<=m; ++i) {
      b = b0+i*db;

   // Calculate the first term of theta
      for (int j=0; j<=n; ++j) {
        double r = b+h*(j+1);
        fi[j] = 1/(r*r*Math.sqrt(fb(r)));
      }
      g1 = simpson(fi, h);

   // Find r_m from 1-b*b/(r*r)-U/E = 0
      double rm = secant(nc, del, b, h);

   // Calculate the second term of theta
      for (int j=0; j<=n; ++j) {
        double r = rm+h*(j+1);
        fi[j] = 1/(r*r*Math.sqrt(f(r)));
      }
      g2 = simpson(fi, h);
      theta[i] = 2*b*(g1-g2);
    }

 // Calculate d theta/d b
    sig = firstOrderDerivative(db, theta, ne);

 // Put the cross section in log form with the exact
 // result of the Coulomb scattering (ruth)
    for (int i=m; i>=0; --i) {
      b = b0+i*db;
      sig[i] = b/(Math.abs(sig[i])*Math.sin(theta[i]));
      double ruth = 1/Math.pow(Math.sin(theta[i]/2),4);
      ruth /= 16;
      double si = Math.log(sig[i]);
      double ru = Math.log(ruth);
      System.out.println("theta = " + theta[i]);
      System.out.println("ln sigma(theta) = " + si);
      System.out.println("ln sigma_r(theta) = " + ru);
      System.out.println();
    }
  }

// Method to achieve the evenly spaced Simpson rule.

  public static double simpson(double y[], double h) {
    int n = y.length-1;
    double s0 = 0, s1 = 0, s2 = 0;
    for (int i=1; i<n; i+=2) {
      s0 += y[i];
      s1 += y[i-1];
      s2 += y[i+1];
    }
    double s = (s1+4*s0+s2)/3;

 // Add the last slice separately for an even n+1
    if ((n+1)%2 == 0)
      return h*(s+(5*y[n]+8*y[n-1]-y[n-2])/12);
    else
      return h*s;
  }

// Method to carry out the secant search.

  public static double secant(int n, double del,
    double x, double dx) {
    int k = 0;
    double x1 = x+dx;
    while ((Math.abs(dx)>del) && (k<n)) {
      double d = f(x1)-f(x);
      double x2 = x1-f(x1)*(x1-x)/d;
      x  = x1;
      x1 = x2;
      dx = x1-x;
      k++;
    }
    if (k==n) System.out.println("Convergence not" +
      " found after " + n + " iterations");
    return x1;
  }

// Method for the 1st-order derivative with the 3-point
// formula.  Extrapolations are made at the boundaries.

  public static double[] firstOrderDerivative(double h,
    double f[], int m) {
    int n = f.length-1;
    double[] y = new double[n+1];
    double[] xl = new double[m+1];
    double[] fl = new double[m+1];
    double[] fr = new double[m+1];

    for (int i=1; i<n; ++i)
      y[i] = (f[i+1]-f[i-1])/(2*h);

 // Lagrange-extrapolate the boundary points
    for (int i=1; i<=(m+1); ++i) {
      xl[i-1] = h*i;
      fl[i-1] = y[i];
      fr[i-1] = y[n-i];
    }
    y[0] = aitken(0, xl, fl);
    y[n] = aitken(0, xl, fr);
    return y;
  }

// The Aitken method for the Lagrange interpolation.

  public static double aitken(double x, double xi[],
    double fi[]) {
    int n = xi.length-1;
    double ft[] = (double[]) fi.clone();
    for (int i=0; i<n; ++i) {
      for (int j=0; j<n-i; ++j) {
        ft[j] = (x-xi[j])/(xi[i+j+1]-xi[j])*ft[j+1]
               +(x-xi[i+j+1])/(xi[j]-xi[i+j+1])*ft[j];
      }
    }
    return ft[0];
  }

// Method to provide function f(x) for the root search.

  public static double f(double x) {
    return 1-b*b/(x*x)-Math.exp(-x/a)/(x*e);
  }

// Method to provide function 1-b*b/(x*x).

  public static double fb(double x) {
    return 1-b*b/(x*x);
  }
}