// lsfit2.java  do an order=order least square fit of data z=F(x,y)
//              provide an interpolation z_new = F1(x_new, y_new)
// data sorted by X primary, Y secondary, Z independent, could be xg,yg
// given nx, ny and nxy = nx * ny,  index  X,Y,Z  [k] = i*ny+j
// needs simeq.java
//
// usage:      lsfit2 mylsfit = new lsfit2(xord, yord, nx, ny, X, Y, Z);
//             for(int i= ...
//               for(int j= ...
//                 your choice of x and y
//                 z = mylsfit.eval2(x,y); // typically includes
//                                         // intermediate values

class lsfit2
{
    int xord, yord, maxord, xyord, nx, ny;
  double C[];

  lsfit2(int xord, int yord, int nx, int ny,
         double X[], double Y[], double Z[])
  {
    double xfct[] = new double[xord+1];
    double yfct[] = new double[yord+1];
    double term[] = new double[(xord+1)*(yord+1)]; // maximum
    double B[][];
    double x, y, z;
    int k;

    this.xord = xord;
    this.yord = yord;
    this.nx = nx;
    this.ny = ny;
    maxord = Math.max(xord,yord); // maximum term order

    xyord = 0;
    for(int i=0; i<=xord; i++)
    {
      for(int j=0; j<=yord; j++)
      {
        if(i+j<=maxord)
	{ 
          xyord++;
          System.out.println("term = x^"+i+"*y^"+j);
	}
      } // end j
    } // end i
    B = new double[xyord][xyord+1];
    C = new double[xyord];
    for(int iii=0; iii<xyord; iii++)
    {
      for(int jjj=0; jjj<xyord+1; jjj++)
      {
	B[iii][jjj] = 0.0;
      } // end jjj
    } // end iii


    for(int ii=0; ii<nx; ii++)
    {
      for(int jj=0; jj<ny; jj++)
      {
        x = X[ii*ny+jj];
	y = Y[ii*ny+jj];
	z = Z[ii*ny+jj];
        xfct[0] = 1.0;
        for(int i=1; i<=xord; i++)
        {
          xfct[i] = xfct[i-1]*x;
        }
        yfct[0] = 1.0;
        for(int j=1; j<=yord; j++)
        {
          yfct[j] = yfct[j-1]*y;
        }
        k = 0;
        for(int i=0; i<=xord; i++)
        {
          for(int j=0; j<=yord; j++)
          {
	    if(i+j<=maxord)
	    {
              term[k++] = xfct[i]*yfct[j];
	    }
          } // end j
        } // end i

        for(int iii=0; iii<xyord; iii++)
	{
          for(int jjj=0; jjj<xyord; jjj++)
	  {
	    B[iii][jjj] += term[iii]*term[jjj];
	  } // end jjj
          B[iii][xyord] += term[iii]*z;
	} // end iii
      } // end jj
    } // end ii
    new simeq(xyord, B, C);
  } // end lsfit2 constructor

  public double eval2(double x, double y)
  {
    double z;
    double xfct[] = new double[xord+1];
    double yfct[] = new double[yord+1];
    int k;

    xfct[0] = 1.0;
    for(int i=1; i<=xord; i++)
    {
      xfct[i] = xfct[i-1]*x;
    }

    yfct[0] = 1.0;
    for(int j=1; j<=yord; j++)
    {
      yfct[j] = yfct[j-1]*y;
    }

    z = 0.0;
    k = 0;
    for(int i=0; i<=xord; i++)
    {
      for(int j=0; j<=yord; j++)
      {
	if(i+j<=maxord)
	{
          z += C[k++]*xfct[i]*yfct[j];
        }
      } // end j
    } // end i
    return z;
  }

    void test(double X[], double Y[], double Z[])
  {
    System.out.println("test xord="+this.xord+", yord="+this.yord+", nx="+this.nx+", ny="+this.ny);

    int nxy = nx*ny;
    double z_true, z_eval, err, maxerr;
    maxerr = 0.0;
    for(int i=0; i<nx;  i++)
    {
      for(int j=0; j<ny; j++)
      {
        z_true = Y[j]*Y[j]+2.0*X[i]*Y[j]+3.0*X[i]*X[i]+
	         4.0*Y[j]+5.0*X[i]+6.0;
        z_eval = eval2(X[i], Y[j]);
        err = Math.abs(z_true-z_eval);
        if(err>maxerr) maxerr = err;
        System.out.println("X["+i+"]="+X[i]+", Y["+j+"]="+Y[j]);
        System.out.println("z_true="+z_true+", z_eval="+z_eval+", err="+err);
      }
    }
    System.out.println("maxerr="+maxerr);
  } // end test

  public static void main(String args[])
  {
    int xord=2;
    int yord=2;
    int nx=3;
    int ny=3;
    double X[] = {1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 4.0, 4.0, 4.0};
    double Y[] = {2.0, 4.0, 5.0, 2.0, 4.0, 5.0, 2.0, 4.0, 5.0};
    double Z[] = new double[nx*ny];

    System.out.println("lsfit2.java running");
    System.out.println("lsfit2 xord="+xord+", yord="+yord+", nx="+nx+", ny="+ny);
    for(int i=0; i<nx;  i++)
    {
      for(int j=0; j<ny; j++)
      {
        Z[i*ny+j] = Y[j]*Y[j]+2.0*X[i]*Y[j]+3.0*X[i]*X[i]+
	    4.0*Y[j]+5.0*X[i]+6.0;
      }
    }
    lsfit2 mylsfit = new lsfit2(xord,yord,nx,ny,X,Y,Z);  
    mylsfit.test(X, Y, Z);
    System.out.println("\n---------------------------------------------\n");
    nx = 4;
    ny = 4;
    double X2[] = {0.5, 0.5, 0.5, 0.5,
                   1.0, 1.0, 1.0, 1.0,
                   2.0, 2.0, 2.0, 2.0,
                   4.0, 4.0, 4.0, 4.0};
    double Y2[] = {0.3, 2.0, 4.0, 5.0,
                   0.3, 2.0, 4.0, 5.0,
                   0.3, 2.0, 4.0, 5.0,
                   0.3, 2.0, 4.0, 5.0};
    double Z2[] = new double[nx*ny];
    for(int i=0; i<nx;  i++)
    {
      for(int j=0; j<ny; j++)
      {
        Z2[i*ny+j] = Y2[j]*Y2[j]+2.0*X2[i]*Y2[j]+3.0*X2[i]*X2[i]+
	    4.0*Y2[j]+5.0*X2[i]+6.0;
      }
    }
    lsfit2 mylsfit2 = new lsfit2(xord,yord,nx,ny,X2,Y2,Z2);  
    mylsfit2.test(X2, Y2, Z2);
    System.out.println("lsfit2.java finished");
  } // end main

} // end class lsfit2