/*************************************************************************/
/* Program file name: Ex2_6.java                                         */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 9, 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.     */
/*                                                                       */
/*************************************************************************/


// Least-squares approximation of the Bessel function.

import java.lang.*;
public class Ex2_6 {
  public static void main(String argv[]) {
    double xi[] = {0, 0.2, 0.4, 0.6, 0.8, 0.9, 1};
    double fi[] = {1, 0.912005, 0.671133, 0.339986, 0.002508,
      -0.142449, -0.260052};
    int n = xi.length-1;
    int m = 21;
    double p[] = new double[n+1];
    double u[][] = orthogonalPolynomialFit(m, xi, fi);

    for (int i=0; i<=n; ++i)
      for (int j=0; j<=m; ++j)
        p[i] += u[j][n+1]*u[j][i];

    for (int i=0; i<=n; ++i) {
      double df = Math.abs(f(xi[i]*xi[i])-p[i]);
      System.out.println("J_(3x) = " + p[i] + " +- " +df + " for x = " + xi[i]);
    }
  }

// 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;
  }

// Method to approximate the Bessel function J_0(3x).

  public static double f(double x) {
    double x2 = x*x;
    double x3 = x2*x;
    double x4 = x2*x2;
    double x5 = x3*x2;
    double x6  = x3*x3;

    return 1 - 2.2499997*x + 1.2656208*x2 - 0.3163866*x3
      + 0.0444479*x4 - 0.0039444*x5 + 0.0002100*x6;
  }
}