// test_eqn1.java   compile    javac -cp . test_eqn1.java
//                  execute    java  -cp . test_eqn1
// solve 2*U(x,y) = C(x,y)  for values of U(x,y) at points in a grid
// U(x,y) = 1+x-y   V(x,y) = x+y
// thus C(x,y) = 2 + 2x - 2y  assume you are given values of C(x,y)
//                            at points on grid and you want to find
//                            values of U(x,y) for any x,y
// needs simeq.java compiled


public class test_eqn1
{
  int nx = 3; // grid
  int ny = 4; // grid
  int k, ii, jj ;
  int row = 0; // typically S(i,j)
  int grid = nx*ny; // points in grid
  int neqn = grid; // two 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;
  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)  
  double err, max_err;
    
  public test_eqn1()
  {
    System.out.println("test_eqn1.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 ");
    // initialize A to all zero
    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);
        System.out.println(" at i="+i+",j="+j+"	Y["+row+"]="+Y[row]);
        // System.out.println("U("+x+","+y+")="+U(x,y)+"  V(x,y)="+V(x,y));
 	A[row][k] = 2.0;   // 2*U(x,y)
	k = k+1;
      } // 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;
    for(int row=0; row<neqn; row++)
    {
      ii = row/ny;
      jj = row-ii*ny;
      err = Math.abs(X[row]-U(xg[ii],yg[jj]));
      if(err>max_err) { max_err = err; }
      System.out.println("at i="+ii+" j="+jj+" X["+row+"]="+X[row]+"  err="+err);
    }    
    System.out.println(" ");
    System.out.println("max_err="+max_err);
    System.out.println("test_eqn1.java finished");
  } // end test_eqn1

  // U(x,y) = 1+x-y   V(x,y) = x+y
  // thus C(x,y)  = 2 + 2*x - 2*y
  double C(double x, double y)
  {
    return 2.0+2.0*x-2.0*y;
  } // end C

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

  int S(int i, int j) // index to one grid
  {
    return ny*i+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_eqn1();
  }
} // end test_eqn1.java