// test_eqn3.java   compile    javac -cp . test_eqn3.java
//                  execute    java  -cp . test_eqn3
// solve 2*U(x,y)+3*V(x,y)+4*W(x,y) = C(x,y)
// solve 4*U(x,y)+2*V(x,y)+3*W(x,y) = D(x,y)
// solve 3*U(x,y)+4*V(x,y)+2*W(x,y) = E(x,y)
// for values of U(x,y), V(x,y), W(x,y) at points in a grid
// U(x,y) = 1+x*x*x-y   V(x,y) = x*x+y*y  W(x,y) = 2-x+2y  unknown
// C, D, E values given
// needs set of 3 equations for 3 variables at each grid
// needs simeq.java compiled


public class test_eqn3
{
  int nx = 2; // grid
  int ny = 3;
  int k, kk, ii, jj;
  int row = 0; // typically S(i,j)
  int grid = nx*ny; // points in grid
  int nvar = 3; // U(x,y)  V(x,y)  W(x,y)
  int neqn = nvar*grid; // nvar unknowns at each grid point
  double xmin = 1.0;
  double xmax = 3.0;
  double ymin = 2.0;
  double ymax = 4.0;
  double x, y, dx, dy, val;
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double A[][] = new double[neqn][neqn]; // A*X=Y  
  double X[] = new double[neqn];  // unknown  
  double Y[] = new double[neqn];  // known from C(x,y),D(x,y),E(x,y)  
  double err, max_err;
    
  public test_eqn3()
  {
    System.out.println("test_eqn3.java running");
    System.out.println("grid nx="+nx+"  ny="+ny);
    System.out.println("xmin="+xmin+"  xmax="+xmax);
    System.out.println("ymin="+ymin+"  ymax="+ymax);
    System.out.println("x grid xg[]");
    xg[0] = xmin;
    dx = (xmax-xmin)/(double)(nx-1);
    System.out.println("dx="+dx+"  xg[0]="+xmin);
    for(int i=1; i<nx; i++)
    {
      xg[i] = xg[i-1]+dx; 
      System.out.println("  xg["+i+"]="+xg[i]);
    } // end i
    
    System.out.println(" ");
    System.out.println("y grid yg[]");
    yg[0] = ymin;
    dy = (ymax-ymin)/(double)(ny-1);
    System.out.println("dy="+dy+"  yg[0]="+yg[0]);
    for(int j=1; j<ny; j++)
    {
      yg[j] = yg[j-1]+dy; 
      System.out.println("  yg["+j+"]="+yg[j]);
    } // end i
    System.out.println(" ");
    System.out.println("build A and Y ");
    for(int i=0; i<neqn; i++)
    {
      for(int j=0; j<neqn; j++)
      {
	A[i][j] = 0.0;
      } // end j
    } // end i
    k = 0;
    for(int i=0; i<nx; i++)
    {
      x = xg[i]; 
      for(int j=0; j<ny; j++)
      {
        y = yg[j];
	row = S(i,j);
	Y[row] = C(x,y);
	Y[row+grid] = D(x,y);
	Y[row+2*grid] = E(x,y);
        System.out.println("U("+x+","+y+")="+U(x,y)+"  V(x,y)="+V(x,y)+
			   "  W(x,y)="+W(x,y));
	
        System.out.println(" at i="+i+",j="+j+"	Y["+row+"]="+Y[row]);
        System.out.println("Y["+(row+grid)+"]="+Y[row+grid]);
        System.out.println("Y["+(row+2*grid)+"]="+Y[row+2*grid]);
 	A[row][k] = 2.0;   // 2*U(x,y)
	A[row][k+1] = 3.0; // 3*V(x,y)
	A[row][k+2] = 4.0; // 4*W(x,y)
 	A[row+grid][k] = 4.0;   // 3*U(x,y)  sum is C(x,y)
        A[row+grid][k+1] = 2.0; // 2*V(x,y)  sum is D(x,y)
        A[row+grid][k+2] = 3.0; // 2*V(x,y)  sum is E(x,y)
	A[row+2*grid][k] = 3.0;   // 3*U(x,y)  sum is C(x,y)
        A[row+2*grid][k+1] = 4.0; // 2*V(x,y)  sum is D(x,y)
        A[row+2*grid][k+2] = 2.0; // 2*V(x,y)  sum is E(x,y)
	k = k+nvar; 
      } // end j
    } // end i
    System.out.println(" ");
    printA();
    System.out.println("run simeq ");
    new simeq(A,Y,X);
    System.out.println("X results from simeq");
    max_err = 0.0;
    k = 0; // k==0 first variable  k==2 second variable
    for(int row=0; row<neqn; row++)
    {
      kk = row/nvar;
      ii = kk/ny;
      jj = kk-ii*ny;
      System.out.println("ii="+ii+"  jj="+jj+"   k="+k);
      if(k==0) { val = U(xg[ii],yg[jj]); }
      if(k==1) { val = V(xg[ii],yg[jj]); }
      if(k==2) { val = W(xg[ii],yg[jj]); }
      err = Math.abs(X[row]-val);
      if(err>max_err) { max_err=err; }
      System.out.println("X["+row+"]="+X[row]+"  err="+err);
      k = k + 1;
      if(k>=nvar) { k=0; }
    } // end row    
    System.out.println(" ");
    System.out.println("max_err="+max_err);
    System.out.println("test_eqn3.java finished");
  } // end test_eqn3

  double U(double x, double y) // for testing and generating matrix
  {
    return 1.0+x*x*x-y;
  } // end U

  double V(double x, double y)
  {
    return x*x+y*y;
  } // end V

  double W(double x, double y)
  {
    return 2.0-x+2.0*y;
  } // end W

  double C(double x, double y)
  {
    return 2.0*U(x,y)+3.0*V(x,y)+4.0*W(x,y);
  } // end C

  double D(double x, double y)
  {
    return 4.0*U(x,y)+2.0*V(x,y)+3.0*W(x,y);
  } // end D

  double E(double x, double y)
  {
    return 3.0*U(x,y)+4.0*V(x,y)+2.0*W(x,y);
  } // end E
      
  int S(int i, int j) // index to one grid
  {
    return i*ny+j;
  } // end S

  void  printA() // non zero values of A
  {
    System.out.println("non zero A");
    for(int i=0; i<neqn; i++)
    {
      for(int j=0; j<neqn; j++)
      {
	if(Math.abs(A[i][j]) > 0.0001)
	{
          System.out.println("A["+i+"]["+j+"]="+A[i][j]);
	} // end if
      } // end j
    } // end i
  } // end printA

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