----------------------------------------------------------------------
-- Aufgabe: partielle Differentialgleichung wie zB Laplace,
--                 diskretisiert ueber Gitter: 4*P(i,j) =   P(i-1,j) + P(i+1,j)
--                                                        + P(i,j-1) + P(i,j+1)
--                                                        - F(i,j)
-- Methode nur zum Vergleich: Standard Iteration(Gauss Seidel Typ)
----------------------------------------------------------------------

with Stringpack; use Stringpack;
procedure Gridsynchro is
   N: constant := 5;
   subtype Full_Grid is Integer range 0 .. N;
   subtype Grid is Full_Grid range 1 .. N-1;   -- without border
   type Real is digits 7;  -- declare our own type
   type Matrix is array(Integer range <>, Integer range <>) of Real;
   P: Matrix(Full_Grid, Full_Grid)  := (others=>(others=>0.0));
   Delta_P: Matrix(Grid, Grid) := (others=>(others=>1.0));
   Tolerance: constant Real := 0.0001;
   Error_Limit: constant Real := Tolerance * Real(N-1)**2;
   Converged: Boolean := False;
   Error_Sum: Real;
   New_P: Real;


   function F(I, J: Grid) return Real is   -- keep it simple
   begin return 0.0; end F;


begin            -- main subprogram
   
   for I in Full_Grid loop  --init P;
      P(I,0) := 1.0; 
   end loop;

   loop
      Error_Sum := 0.0;
      for I in Grid loop
	 for J in Grid loop
	    New_P := 0.25 * (P(I-1, J) + P(I+1, J) + P(I, J-1)
			     + P(I, J+1) -  F(I, J));
	    Delta_P(I, J) := New_P - P(I, J);
	    P(I, J) := New_P;
	    Error_Sum := Error_Sum + Delta_P(I, J)**2;
	 end loop;
      end loop;
      Converged := Error_Sum < Error_Limit;  
      exit when Converged;
   end loop;

   for I in Full_Grid loop        -- output results
      Outbuffer := Varnull;
      for J in Full_Grid loop
	 Outbuffer := Outbuffer & Cvfs(Float(P(I,J)),1,3) & " ";
      end loop;
      Print;
   end loop;


end Gridsynchro;