/*************************    Program 2.3    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define NMAX 21

main()
/* Main program for a linear fit of the Millikan experimental
   data on the fundamental charge e_0 from e_n = e_0*n + de.
   Copyright (c) Tao Pang 1997. */
{
int i,j,n,m;
double e0,de,sum1,sum2;
double x[NMAX],a[NMAX],u[NMAX][NMAX];
double f[]={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 pfit();

n = 15; m = 2;
for (i = 0; i < n; ++i)
  {
  x[i] = 4+i;
  }
for (i = 0; i < m; ++i)
  {
  a[i] = 0;
  }
for (i = 0; i < m; ++i)
  {
  for (j = 0; j < n; ++j)
    {
    u[i][j] = 0;
    }
  }

pfit(n,m,x,f,a,u);
sum1 = 0; sum2 = 0;
for (i = 0; i < n; ++i)
  {
  sum1 = sum1+u[0][i]*u[0][i];
  sum2 = sum2+x[i]*u[0][i]*u[0][i];
  }
e0 = a[1];
de = a[0]-a[1]*sum2/sum1;
printf("%8.4lf %8.4lf\n", e0,de);
}

void pfit (n,m,x,f,a,u)
/* Subroutine generating orthonormal polynomials U(M,N) up to
   (M-1)th order and coefficients A(M), for the least squares
   approximation of the function F(N) at X(N). Other variables
   used: G(K) for g_k, H(K) for h_k, S(K) for <u_k|u_k>.
   Copyright (c) Tao Pang 1997. */

int n,m;
double x[],f[];
double a[],u[NMAX][NMAX];
{
int i,j;
double stmp;
double s[NMAX],g[NMAX],h[NMAX];

if (n > NMAX)
  {
  printf ("Too many points\n");
  exit (1);
  }
if (m > n)
  {
  printf ("Order too high\n");
  exit (2);
  }

/* Set up the zeroth order polynomial u_0 */

s[0] = 0;  g[0] = 0;  a[0] = 0;

for (i = 0; i < n; ++i)
  {
  u[0][i] = 1;
  stmp   = u[0][i]*u[0][i];
  s[0]   = s[0]+stmp;
  g[0]   = g[0]+x[i]*stmp;
  a[0]   = a[0]+u[0][i]*f[i];
  }
g[0] = g[0]/s[0];
a[0] = a[0]/s[0];
h[0] = 0;

/* Set up the first order polynomial u_1 */

for (i = 0; i < n; ++i)
  {
  u[1][i] = x[i]*u[0][i]-g[0]*u[0][i];
  s[1]   = s[1]+u[1][i]*u[1][i];
  g[1]   = g[1]+x[i]*u[1][i]*u[1][i];
  h[1]   = h[1]+x[i]*u[1][i]*u[0][i];
  a[1]   = a[1]+u[1][i]*f[i];
  }
g[1] = g[1]/s[1];
a[1] = a[1]/s[1];
h[1] = h[1]/s[0];

/* Higher order polynomials u_k from the recursive relation */

if (m >= 3)
  {
  for (i = 1; i < m-1; ++i)
    {
    s[i+1] = 0;  g[i+1] = 0;  a[i+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]   = s[i+1]+u[i+1][j]*u[i+1][j];
      g[i+1]   = g[i+1]+x[j]*u[i+1][j]*u[i+1][j];
      h[i+1]   = h[i+1]+x[j]*u[i+1][j]*u[i][j];
      a[i+1]   = a[i+1]+u[i+1][j]*f[j];
      }
    g[i+1] = g[i+1]/s[i+1];
    a[i+1] = a[i+1]/s[i+1];
    h[i+1] = h[i+1]/s[i];
    }
  }
}