// test_psimeq.java

public class test_psimeq
{
  double rms_err, avg_err, max_err;
  int idata = 0; // reset after each use of the data set
  int n_data = 18; // number of data points 
  double x_data[] = new double[n_data];
  double y_data[] = new double[n_data];
  boolean debug = false;

  public test_psimeq()
  {
    int n;

    System.out.println("test_psimeq.java");
    for(int i=0; i<n_data; i++) x_data[i]=(double)i/10.0;
    y_data[0]=0.0;  
    y_data[1]=6.0;  
    y_data[2]=14.1;  
    y_data[3]=5.01;  
    y_data[4]=4.326;  
    for(int i=5; i<n_data; i++) y_data[i]=y_data[i-1];
    y_data[n_data-1]=0.0;
    if(debug)
      for(int i=0; i<n_data; i++)
	System.out.println("i="+i+", x="+x_data[i]+", y="+y_data[i]);

    // sample polynomial least square fit, nth power 
    for(n=4; n<=n_data; n++) // was n_data = 5 for test
    { 
      double A[][] = new double[n][n];
      double C[] = new double[n];
      double Y[] = new double[n];
      System.out.println("fit data to "+(n-1)+" degree polynomial");
      fit_pn(n, A, Y, C);
      check_pn(n, C);
      for(int i=0; i<n; i++) System.out.println("C["+i+"]="+C[i]);
      System.out.println("max_err="+max_err+", avg_err="+avg_err+
                         ", rms_err="+rms_err);
      if(debug) System.out.println(" ");
    }
  } // end test_psimeq

  int data_set(double yx[]) // returns 1 for data, 0 for end  
  {                                  // sets value of y for value of x 
    if(idata>=n_data) {idata=0; return 0;} // ready for check 
    yx[1] = x_data[idata];
    yx[0] = y_data[idata];
    idata++;  
    return 1;
  } // end data_set 

  void fit_pn(int n, double A[][], double Y[], double C[])
  {
    int i, j, k;
    double x, y, t;
    double yx[] = new double[2];
    double pwr[] = new double[2*n]; // at least 2n 
  
    for(i=0; i<n; i++)
    {
      for(j=0; j<n; j++)
      {
        A[i][j] = 0.0;
      }
      Y[i] = 0.0;
    }
    while(data_set(yx)==1)
    {
      y = yx[0];
      x = yx[1];
      pwr[0] = 1.0;
      for(i=1; i<2*n; i++) pwr[i] = pwr[i-1]*x;
      for(i=0; i<n; i++)
      {
        for(j=0; j<n; j++)
        {
          A[i][j] = A[i][j] + pwr[i]*pwr[j];
        }
        Y[i] = Y[i] + y*pwr[i];
      }
    }
    new psimeq(A, Y, C);
    if(debug) for(i=0; i<n; i++) System.out.println("C["+i+"]="+C[i]);
  } // end fit_pn 

  void check_pn(int n, double C[])
  {
    double x, y, ya, diff;
    double yx[] = new double[2];
    double sumsq = 0.0;
    double sum = 0.0;
    double maxe = 0.0;
    double xmin = 0.0;
    double xmax = 0.0;
    double ymin = 0.0;
    double ymax = 0.0;
    double xbad = 0.0;
    double ybad = 0.0;
    int i, k, imax=0;
  
    k = 0;
    while(data_set(yx)==1)
    {
      y = yx[0];
      x = yx[1];
      if(k==0)
      {
        xmin=x;
        xmax=x;
        ymin=y;
        ymax=y;
        imax=0;
        xbad=x;
        ybad=y;
      }
      if(x>xmax) xmax=x;
      if(x<xmin) xmin=x;
      if(y>ymax) ymax=y;
      if(y<ymin) ymin=y;
      k++;
      ya = C[n-1]*x;
      for(i=n-2; i>0; i--)
      {
        ya = (C[i]+ya)*x;
      }
      ya = ya + C[0];
      diff = Math.abs(y-ya);
      if(diff>maxe)
      {
        maxe=diff;
        imax=k;
        xbad=x;
        ybad=y;
      }
      sum = sum + diff;
      sumsq = sumsq + diff*diff;
    }
    if(debug)
    System.out.println("check_pn k="+k+", xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax);
    max_err = maxe;
    avg_err = sum/(double)k;
    rms_err = Math.sqrt(sumsq/(double)k);
    if(debug)
    System.out.println("max="+max_err+" at i="+imax+", xbad="+xbad+
                       ", ybad="+ybad);
  } // end check_pn 

  public static void main (String[] args)
  {
    new test_psimeq();
  } // end main

} // end class test_psimeq