/*
 * Heat.java by Richard J. Davies
 * from `Introductory Java for Scientists and Engineers'
 * chapter: `Physical Modelling'
 * section: `Example - Heat Flow in a Solid'
 *
 * This example simulates a square of thermally conducting material which is
 * placed in an environment in which various thermal conditions are applied
 * to each side. We use Gauss-Seidel iteration to find the stable state of
 * the system.
 */
import uk.co.jscieng.*;
import java.awt.*;

public class Heat
{
  // The number of grid squares along each side of
  // the piece of material.
  public static final int SIZE = 100;
  
  // The number of Gauss-Seidel iterations to
  // perform between redraws.
  public static final int ITERPERSTEP = 100;
  
  // The number of redraws to perform.
  public static final int STEPS=10;
  
  
  // A method that applies the boundary
  // conditions to the grid.
  public static void setBCs(double[][] state)
  {
    int i;
    
    // The left hand side is held at 100.
    for (i=0; i<SIZE; i++)
      state[0][i] = 100;
    
    // The bottom is held at 50.
    for (i=0; i<SIZE; i++)
      state[i][0] = 50;
    
    // Half of the right is held at 50.
    for (i=0; i<SIZE/2; i++)
      state[SIZE-1][i] = 50;
    
    // The rest of the right is held at 0.
    for (i=SIZE/2; i<SIZE; i++)
      state[SIZE-1][i] = 0;
  }
  
  
  // A method to perform Gauss-Seidel iteration.
  // Note that we overcorrect by a factor of
  // 1.8 because we are using SOR.
  public static void iterate(double[][] state)
  {
    double newstate;
    int i, j;
    
    // Iterate each of the squares in the middle
    // of the material. These are the average of
    // four neighbours.
    for (i=1; i<SIZE-1; i++)
      for (j=1; j<SIZE-1; j++)
      {
        newstate = (state[i-1][j] + state[i][j-1]
                    + state[i][j+1] + state[i+1][j])/4;
        state[i][j] += 1.8 * (newstate - state[i][j]);
      }
    
    // Iterate squares on the top edge. These are
    // the average of three neighbours.
    for (i=1; i<SIZE-1; i++)
    {
      newstate = (state[i-1][SIZE-1] + state[i][SIZE-2]
                  + state[i+1][SIZE-1])/3;
      state[i][SIZE-1] += 1.8 * (newstate - state[i][SIZE-1]);
    }
  }
  
  
  // A method used by the graphical output.
  // We return a suitable colour for a temperature t.
  public static Color colorFor(double t)
  {
    if (t>86)
      return Color.magenta;
    else if (t>71)
      return Color.red;
    else if (t>57)
      return Color.orange;
    else if (t>43)
      return Color.yellow;
    else if (t>29)
      return Color.green;
    else if (t>14)
      return Color.cyan;
    else
      return Color.blue;
  }
  
  
  // The main method.
  public static void main(String[] argv)
  {
    int i, j, k;
    
    // Create grid and initialize all squares to 50
    // since most will end up near here.
    double[][] state = new double[SIZE][SIZE];

    for (i=0; i<SIZE; i++)
      for (j=0; j<SIZE; j++)
        state[i][j] = 50;

    // Apply the actual BCs to the edge squares.
    setBCs(state);
    
    // Show the graphics window.
    SciGraph.showGraph(0, SIZE, 0, SIZE, SIZE, SIZE);

    // For each redrawing step that we will do.
    for (i=0; i<STEPS; i++)
    {
      // Perform the correct number of iterations.
      for (j=0; j<ITERPERSTEP; j++)
        iterate(state);
      
      // Plot the temperature distribution,
      // displaying this as a colour gradient.
      for (j=0; j<SIZE; j++)
        for (k=0; k<SIZE; k++)
          SciGraph.point(j, k, colorFor(state[j][k]));
    }
  }
}