/*************************************************************************/
/* Program file name: millikan.c                                         */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 7, 2006                                           */
/*                                                                       */
/* (1) This C 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.     */
/*                                                                       */
/*************************************************************************/

#include <stdio.h>
#include <math.h>

/* An example of applying the least-squares approximation
   to the Millikan data on a straight line q=k*e+dq. */

main() {
  int i, j, n=14, m=21;
  double e, dq, sum;
  double u[m+1][n+2];
  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};
  void orthpoly();

  orthpoly(n, m, k, q, u);
  sum = 0;
  for (i=0; i<=n; ++i) sum += k[i];
  e = u[1][n+1];
  dq = u[0][n+1]-u[1][n+1]*sum/(n+1);
  printf("Fundamental charge: %8.4lf\n", e);
  printf("Estimated error: %8.4lf\n", dq);
}

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

void orthpoly(n, m, x, f, u)

int n, m;
double x[], f[];
double u[m+1][n+2];

{
  int i, j;
  double stmp;
  double s[n+1], g[n+1], h[n+1];

  if (m>n) {
    m = n;
    printf ("The highest power is adjusted to %d\n", n);
  }

/* Set up the zeroth-order polynomial */

  s[0] = 0;
  g[0] = 0;
  u[0][n+1] = 0;
  h[0] = 0;
  for (i=0; i<=n; ++i) {
    u[0][i] = 1;
    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] /= s[0];
  u[0][n+1] /= s[0];

/* Set up the first order polynomial */

  s[1] = 0;
  g[1] = 0;
  u[1][n+1] = 0;
  h[1] = 0;
  for (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] /= s[1];
  u[1][n+1] /= s[1];
  h[1] /= s[0];

/* Obtain the higher-order polynomials recursively */

  if (m >= 2) {
    for (i=1; i<m; ++i) {
      s[i+1] = 0;
      g[i+1] = 0;
      u[i+1][n+1] = 0;
      for (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] /= s[i+1];
      u[i+1][n+1] /= s[i+1];
      h[i+1] /= s[i];
    }
  }
}