/* 	imhf.java
	
	Implicit heat flow equation, finite differenced.
	Adapted from Jon Claerbout's Imaging the Earth's Interior (1985), p. 98,
	by John N. Louie, 28 November 1998.
*/

/* Include the usual classes */
import java.util.*;
import java.io.*;
import java.lang.Math;

/* Declare a class with a main() method to contain the finite-difference code. */
/* When you change the code in this file, make a copy to a new file with
    a new name like "ex6imhf.java", and be sure to make the same change
    to the class name below - e.g. "class ex6imhf" */
class imhf
	{
	/* The main() method compiles into a file named "imhf.class" that
	    you run in the Java virtual machine as a stand-alone application
	    with a command such as "java imhf"; or however you have
	    renamed it. */
	public static void main(String args[]) /* No runtime arguments needed */
		{
		System.out.println("Starting imhf...");
		/* Declare variables */
		int nx, it, ix;
		float alpha, a, q[], d[], e[], f[];
		/* Dimension arrays */
		nx = 12;
		q = new float[nx];
		d = new float[nx];
		e = new float[nx];
		f = new float[nx];
		
		a = 8F;
		System.out.println("a = " + a);
		alpha = 0.5F*a;
	
		/* Initial temperature is a step in x */
		for (ix=0; ix<6; ix++)
			q[ix] = 0F;
		for (ix=6; ix<12; ix++)
			q[ix] = 1F;
	
		/* Iterate in time */
		for (it=0; it<4; it++)
			{
			/* Print current state of array */
			for (ix=0; ix<nx; ix++)
				System.out.print(" " + q[ix]);
			System.out.println("");
	
			/* Set sides of implicit array */
			d[0] = 0F;
			d[nx-1] = 0F;
	
			/* Apply top half of differencing star to array */
			for (ix=1; ix<nx-1; ix++)
				d[ix] = q[ix] + alpha*(q[ix-1] - 2F*q[ix] + q[ix+1]);
				
			/* Solve the tridiagonal system */
			rtris(nx, alpha, -alpha, (1.0F+2.0F*alpha), -alpha, alpha,
				d, q, e, f);
			}
		
		/* End of imhf's main() application method. */
		}

	/*	Real tridiagonal equation solver for the system:

		|l c 0 . . . . . 0| | q[0] |   | d[0] |
		|a b c 0 . . . . .| |  .   |   |  .   |
		|0 a b c 0 . . . .| |  .   |   |  .   |
		|. .` ` `  . . . .| |  .   |   |  .   |
		|. . .` ` `  . . .| |  .   | = |  .   |
		|. . . .` ` `  . .| |  .   |   |  .   |
		|. . . . 0 a b c 0| |  .   |   |  .   | 
		|. . . . . 0 a b c| |  .   |   |  .   |
		|0 . . . . . 0 a r| |q[n-1]|   |d[n-1]|

		Adapted from Claerbout (1985) p. 98-99,
		by John N. Louie, 28 November 1998.

		Arguments:	n	dimension of arrays
				endl	the corner value l (above)
				a,b,c	diagonals of the coefficient matrix
				endr	the corner value r (above)
				d	float[] data array of n elements
				q	float[] solution array of n elements
				e,f	float[] temporary iteration arrays
					of n elements each */
	static void rtris(int n, float endl, float a, float b, float c, float endr,
		float[] d, float[] q, float[] e, float[] f)
		{
		int i;
		float den;
		
		/* Initialize unknowns e and f */
		e[0] = -a/endl;
		f[0] = d[0]/endl;

		/* Recursions of equations 19a and 19b */
		for (i=1; i<n-1; i++)
			{
			den = b + c*e[i-1];
			e[i] = -a/den;
			f[i] = (d[i] - c*f[i-1])/den;
			}

		/* Right-hand boundary condition, equation 20 */
		q[n-1] = (d[n-1] - c*f[n-2])/(endr + c*e[n-2]);

		/* Recursion for equation 21 */
		for (i=n-2; i>=0; i--)
			q[i] = e[i]*q[i+1] + f[i];
		
		/* End of rtris static method of class imhf */
		}

	/* End of methods in class imhf */
	}