/* 
 *  Laplace.java
 *  An applet to compute the solution of a (built-in) 2nd order PDE.
 *  The solution method is Jacobi iteration, and the partial solutions
 *  are displayed as they are computed, until the errot tolerance is reached.
 */

import java.awt.*;
import java.awt.event.*;

public class Laplace extends java.applet.Applet implements Runnable 
//				    implements ActionListener
{
  Thread thread;
  int R = 64;		// number of rows in the solution matrix
  int C = 128;		// number of columns
  double[][] phi0 = new double [R][C];	
  double[][] phi1 = new double [R][C];	
			// solution and initial matrices

  double tolerance = .001;  // error tolerance
  int n = 10;		// number of iterations between error computation
			// and displays of solution

  int npix;		// number of pixels displayed per matrix element
			//     in each direction
  int minz;  	        // min value of solution in range for color map
  int maxz;	        // max value of solution in range

  Image buffer;		// off screen buffer
  Graphics og;		// off screen graphics
  

  public void init () 
  {
    int i,j;
    //initialize initial matrix to 0
    for (i=0; i<R; i++)
    for (j=0; j<C; j++)
      phi0[i][j] = 0.0;

    inittestcase ();

    buffer = createImage( getSize().width, getSize().height );
    og = buffer.getGraphics(); 
  }

  public void initpoly ( )
  {  int i, j;

    // initialize data bounds for display
     npix = 4;		// number of pixels displayed per matrix element
     minz = -32;	// min value of solution in range for color map
     maxz = 67200;	// max value of solution in range

    // initialize boundary values of initial matrix
    i=0; for(j=0; j<C; j++)
		phi0[i][j] = -1.4 + 0.12*i - 0.23*j 
			+0.001 * Math.pow(i,2) * Math.pow(j,2);
    i=R-1; for(j=0; j<C; j++)
		phi0[i][j] = -1.4 + 0.12*i - 0.23*j 
			+0.001 * Math.pow(i,2) * Math.pow(j,2);
    j=0; for(i=0; i<R; i++)
		phi0[i][j] = -1.4 + 0.12*i - 0.23*j 
			+0.001 * Math.pow(i,2) * Math.pow(j,2);
    j=C-1; for(i=0; i<R; i++)
		phi0[i][j] = -1.4 + 0.12*i - 0.23*j 
			+0.001 * Math.pow(i,2) * Math.pow(j,2);
  }     

  public void inittestcase ( )
  {  int i, j;

    // initialize data bounds for display
     npix = 4;		// number of pixels displayed per matrix element
     minz = 0;	// min value of solution in range for color map
     maxz = 1;	// max value of solution in range

    // initialize boundary values of initial matrix
    i=0; for(j=0; j<C; j++)
		phi0[i][j] = 1.0;
    i=R-1; for(j=0; j<C; j++)
		phi0[i][j] = 1.0;
    j=0; for(i=0; i<R; i++)
		phi0[i][j] = 1.0;
    j=C-1; for(i=0; i<R; i++)
		phi0[i][j] = 1.0;
  }     

  public void start() 
  {
    if ( thread == null ) 
    {
      thread = new Thread( this );
      thread.start();
    }
  }


  public void stop() 
  {
    if ( thread != null ) 
    {
      thread.stop();
      thread = null ;
    }
  }


  public void run() 
  { 
    double eps = 1.0;	// computed error between iterations,
                        // initialized to some large value
    while ( eps > tolerance ) 
    {
      // phi0 has current solution, compute phi1 as new solution
      // after n iterations
      Jacobi(n);  	

      // compare phi0 and phi1 for error
      eps = 0.0;
      for (int i=1; i<(R-1); i++)
      for (int j=1; j<(C-1); j++)
          eps = Math.max(eps, Math.abs(phi1[i][j] - phi0[i][j]));

      // move phi1 to be new current solution phi0
      for (int i=0; i<R; i++)
      for (int j=0; j<C; j++)
    	  phi0[i][j] = phi1[i][j];  

      /*
      System.out.println ( "phi[1][C/2]" + phi0[1][C/2]);
      System.out.println ( "phi[R/2][C/2]" + phi0[R/2][C/2]+"\n");
      */
      repaint();  	// displays phi0
    }
    System.out.println ("Done with epsilon = " + eps);
  }

  public void update( Graphics g ) 
  {
    paint(g);
  }

  public void paint( Graphics g ) 
  { int colorno;
    double diffz, scalez;
    // scale array values to color numbers between 0 and 255
    // paint array in shades of blue
    diffz = maxz - minz;
    scalez = ((double)255)/diffz;
    for (int i=0; i<R; i++)
    for (int j=0; j<C; j++)
    { colorno = (int)Math.floor((phi0[i][j]-minz) * scalez);
      og.setColor(new Color(0, 0, colorno));
      og.fillRect(j*npix, i*npix, npix, npix);
    }
    g.drawImage( buffer, 0, 0, this );
  }

  void Jacobi(int n)
  { // starts with phi0 and computes new phi1 after n iterations,
    //     leaving phi0 unchanged
    double[][] w = new double [R][C];	

    for (int i=0; i<R; i++)
    for (int j=0; j<C; j++)
      { w[i][j]= phi0[i][j];
        phi1[i][j]= phi0[i][j];
       }

    for(int k=0; k<n; k++)
    {
    // compute new phi1 for interior points of the grid
      for (int i=1; i<(R-1); i++)
      for (int j=1; j<(C-1); j++)
	{
	phi1[i][j] = 0.25 *(  w[i-1][j]
			    + w[i+1][j]
			    + w[i][j-1]
			    + w[i][j+1] );
	}

     // copy new phi1 to w to repeat iteration
      for (int i=0; i<R; i++)
      for (int j=0; j<C; j++)
          w[i][j]= phi1[i][j];
    }
  }
}