(*  Takes as input the table fcurrent[nx,ny] and the *)
(*  boundary defined by the table boundary[nx,ny]    *)
(*  and creates the table fnew[nx,ny] that is the    *)
(*  result of one step of the relaxation technique   *)
BeginPackage["RelaxMethod`"]

RelaxIteration::usage = 
	"RelaxIteration[f, boundary, nx, ny] performs one iteration of the 
	relaxation method to solve Laplace's equation for the function f
	expressed as a table with x indices 1->nx and y index 1->ny.  The
	{nx,ny} table boundary defines the boundary of the function.  
	value 0 implies a [i,j] is in the interior, while a value 1 implies
	the [i,j] is a boundary point."

RelaxSolve::usage = 
	"RelaxSolve[phi, boundary, nx, ny, maxsteps, tolerance] 
	uses the relaxation technique to solve for
	Laplaces equation for the function phi[[i,j]] 
	defined over the points i=1 to nx and j=1 to ny.  The boundary function
	boundary[[i,j]] = 0 indicates [i,j] is in the interior of the 
	region while boundary[[i,j]] ^= 0 indicates [i,j] is on the 
	boundary or exterior of the region.  The solution is iterated
	until the last fractional change in the function is less than 
	tolerance of its previous value."

Begin["`Private`"]

RelaxIteration[fcurrent_, boundary_, nx_, ny_] :=
	Module[ {i,j,fnew},

	fnew = Table[0, {i,nx}, {j,ny}];
	fnew = fcurrent;   (*  Initialize the new value for f *)

	Do[
	    If[boundary[[i,j]]==0,

		(*  In the interior of the region  *)
		fnew[[i,j]] = 0.125 * (
			fcurrent[[i-1,j-1]] +
			fcurrent[[i-1,j]] +
			fcurrent[[i-1,j+1]] +
			fcurrent[[i,j-1]] +
			fcurrent[[i,j+1]] +
			fcurrent[[i+1,j-1]] +
			fcurrent[[i+1,j]] +
			fcurrent[[i+1,j+1]] ),

		(*  On the boundary  *)
		fnew[[i,j]] = fcurrent[[i,j]] ],

	    {j,1,ny}, {i,1,nx} ];

	fnew 
	]

RelaxSolve[phi_, boundary_, nx_, ny_, maxiterations_, tolerance_] :=
	Module[ {i, j, phiold, phinew, iterations, maxvalue, 
			difference, ntoobig},

	phinew = Table[0, {i,nx}, {j,ny}];
	phiold = Table[0, {i,nx}, {j,ny}];
	iterations = 0;
	phiold = phi;

	While[ iterations<maxiterations,

	    (* Make a step first *)
	    phinew = RelaxIteration[phiold, boundary, nx, ny];
	    iterations++;

	    (* Determine how much each element changed *)
	    difference = phinew - phiold;
	    maxvalue = Max[phinew];
	    ntoobig = 0;
	    Do [
		If[ Abs[ difference[[i,j]] ] > 
			tolerance Abs[ phiold[[i,j]] ],
		    If[ Abs[ difference[[i,j]] ] >
			0.001 maxvalue,
		        ntoobig++ 
			]
		    ], 
		{i,1,nx}, {j,1,ny} 
		];

	    Print["RelaxSolve::Iterating: At iteration ",
		iterations,", ", ntoobig," points out of tolerance"];

	    phiold = phinew;		(* Keep the new step... *)

	    (*  If all within the tolerance, end iterations.  
		Note that we require at least 3 iterations in order 
		to avoid the initial case where the initial potential
		is zero everywhere...  *)
	    If [ ntoobig == 0 && iterations > 2,
		Break[ ] 
		]
	    ];

	phinew 	(*  Return most up-to-date value for f  *)
	]

End[ ]		(*  End of Private stuff... *)
		    
EndPackage[ ]
