/***************************************************************
* EqHeat.java: Solve heat equation using with finite differences
* copyright RH Landau, Oregon State Univ, 8/99 *  
* comment: Output data saved in file in 3D grid format of gnuplot
* modified 8/01 DH McIntyre to create ouput file for Visad spreadsheet
***************************************************************/ 
import java.io.*;                      // Import input-output Library 
public class EqHeat 
{ 
// Class constants  
public static final int NX = 101;                        // grid size  
public static final int NTIMES = 30000;              // No iterations  
public static final double THC = 0.12;         // thermal conductivity
public static final double SPH = 0.113;             // specific heat  
public static final double RHO = 15.;                    // density  
// method must throw exceptions if does not handle them. 
public static void main(String[] argv) throws IOException, 
					FileNotFoundException  
	{  
	int ix, t; 
	double[] [] T = new double[NX][2]; // 2-D array declared, allocated  
	double cons; 
	//Two options for output files are given, comment one out.
	// Open file EqHeat.csv for  VisAd output data, use JDK 1.1 methods  
	//PrintWriter q = new PrintWriter( new FileOutputStream("EqHeat.csv"), true); 
	/*Print two header lines for VisAd csv (comma separated values) file
	* First header gives mapping of domain(x,y) to range (z)
	* Second header tells which data appears
	* We only ouput z here, so x and y values are assume to go from 
	* 0 to (N-1) where N is the number of values (per line for x, # lines for y)
	*/
	//q.println("(position,time)->(temperature)");
	//q.println("temperature");

	// OR Open file EqHeat.dat for  gnuplot output data, use JDK 1.1 methods 
	PrintWriter q = new PrintWriter( new FileOutputStream("EqHeat.dat"), true); 

	for(ix=1; ix < NX-1; ix++) T[ix][0]=100.; // t=0 all points at 100 C  
	for(t=0; t < 2; t++) // except the endpoints  
		{ 
		T[0][t] = 0.; 
		T[NX-1][t] = 0.; 
		}  
	cons=THC/(SPH*RHO); // material constants  
	for(t=1; t <= NTIMES; t++) // loop over max timesteps  
		{ 
		for(ix=1; ix < NX-1; ix++) // loop over space  
			{ 
			T[ix][1] = T[ix][0]+cons*(T[ix+1][0]+T[ix-1][0]-2.0*T[ix][0]); 
			} 
		if((t%1000==0) || (t==1)) // save every 1000 time steps 
			{ 
			for(ix=0 ; ix < NX; ix++) 
				{
				//q.print(T[ix][1]+" , "); //VisAd: output data with commas
				q.println(T[ix][1]+" ");    //gnuplot: output data with spaces
				}
				q.println();	// end line of output (same for both)

			} 
		for(ix=1; ix < NX-1; ix++) T[ix][0]=T[ix][1]; //shift new values to old 
		}  //end time loop
	//System.out.println("data stored in EqHeat.csv");
	System.out.println("data stored in EqHeat.dat");
	}  //end main
}//end class