///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Millikan.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 applying the least-squares approximation
// to the Millikan data on a straight line q=k*e+dq.

import java.lang.*;
public class Millikan {
  public static void main(String argv[]) {
    double k[] = {4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
      15, 16, 17, 18};
    double q[] = {6.558, 8.206, 9.880, 11.50, 13.14,
      14.81, 16.40, 18.04, 19.68, 21.32, 22.96, 24.60,
      26.24, 27.88, 29.52};
    int n = k.length-1;
    int m = 21;
    double u[][] = orthogonalPolynomialFit(m, k, q);
    double sum = 0;
    for (int i=0; i<=n; ++i) sum += k[i];
    double e = u[1][n+1];
    double dq = u[0][n+1]-u[1][n+1]*sum/(n+1);
    System.out.println("Fundamental charge: " + e);
    System.out.println("Estimated error: " + dq);
  }

// Method to generate the orthogonal polynomials and the
// least-squares fitting coefficients.

  public static double[][] orthogonalPolynomialFit
    (int m, double x[], double f[]) {
    int n = x.length-1;
    double u[][] = new double[m+1][n+2];
    double s[] = new double[n+1];
    double g[] = new double[n+1];
    double h[] = new double[n+1];

 // Check and fix the order of the curve
    if (m>n) {
      m = n;
      System.out.println("The highest power"
        + " is adjusted to: " + n);
    }

 // Set up the zeroth-order polynomial
    for (int i=0; i<=n; ++i) {
       u[0][i] = 1;
       double stmp = u[0][i]*u[0][i];
       s[0] += stmp;
       g[0] += x[i]*stmp;
       u[0][n+1] += u[0][i]*f[i];
    }
    g[0] = g[0]/s[0];
    u[0][n+1] = u[0][n+1]/s[0];

 // Set up the first-order polynomial
    for (int i=0; i<=n; ++i) {
      u[1][i] = x[i]*u[0][i]-g[0]*u[0][i];
      s[1] += u[1][i]*u[1][i];
      g[1] += x[i]*u[1][i]*u[1][i];
      h[1] += x[i]*u[1][i]*u[0][i];
      u[1][n+1] += u[1][i]*f[i];
    }
    g[1] = g[1]/s[1];
    h[1] = h[1]/s[0];
    u[1][n+1] = u[1][n+1]/s[1];

 // Obtain the higher-order polynomials recursively
    if (m >= 2) {
      for (int i=1; i<m; ++i) {
        for (int j=0; j<=n; ++j) {
          u[i+1][j] = x[j]*u[i][j]-g[i]*u[i][j]
            -h[i]*u[i-1][j];
          s[i+1] += u[i+1][j]*u[i+1][j];
          g[i+1] += x[j]*u[i+1][j]*u[i+1][j];
          h[i+1] += x[j]*u[i+1][j]*u[i][j];
          u[i+1][n+1] += u[i+1][j]*f[j];
        }
        g[i+1] = g[i+1]/s[i+1];
        h[i+1] = h[i+1]/s[i];
        u[i+1][n+1] = u[i+1][n+1]/s[i+1];
      }
    }
    return u;
  }
}