-----------------------------------------------------------------------
-- package Runge_Coeffs_pd_5, coefficients for Prince Dormand Runge Kutta
-- Copyright (C) 2008-2009 Jonathan S. Parker.
--
-- This program is free software: you can redistribute it and/or modify
-- it under the terms of the GNU General Public License as published by
-- the Free Software Foundation, either version 3 of the License, or
-- (at your option) any later version.
--
-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-- GNU General Public License for more details.
-- You should have received a copy of the GNU General Public License
-- along with this program.  If not, see  http://www.gnu.org/licenses/
--
-- As a special exception, if other files instantiate generics from
-- this unit, or you link this unit with other files to produce an
-- executable, this  unit  does not  by itself cause  the resulting
-- executable to be covered by the GNU General Public License. This
-- exception does not however invalidate any other reasons why the
-- executable file  might be covered by the  GNU Public License.
-----------------------------------------------------------------------
with Text_IO; use Text_IO;

package body runge_coeffs_pd_5 is

    Print_to_Screen_Desired : constant Boolean := False;
  --Print_to_Screen_Desired : constant Boolean := True;

    -----------------------
    -- Test_Runge_Coeffs --
    -----------------------

    procedure Test_Runge_Coeffs
    is
       Factors : array(Stages, Stages) of Real := (others => (others => 0.0));
       Factorial, Sum, Error, Max_Error : Real := +0.0;
    begin

       for I in Stages loop
          Factors(I,Stages'First) := 1.0;
       end loop;
       --    K(I = 1) = h*F*Y                         J = 1..0
       --
       --    K(I = 2) = h*F*(Y + K(1)*A21)               J = 1..1
       --
       --    K(I = 3) = h*F*(Y + K(1)*A31 + K(2)*A32)          J = 1..2
       --
       --    K(I = 4) = h*F*(Y + K(1)*A41 + K(2)*A42 + K(3)*A43)     J = 1..3

       for I in stages range 1..Stages'Last loop

          for J in Stages range Stages'First..I-1 loop
          for Order in Stages range Stages'First+1..J+1 loop
             Factors(I,Order) := Factors(I,Order) + A(I)(J)*Factors(J,Order-1);
          end loop;
          end loop;

       end loop;

       -- Now should have SUM(B(i) * Factors(i,n)) = 1/n! up to the order of the
       -- the Taylors series.

       for Order in Stages range 1..5 loop
          Sum := 0.0;
          for I in Stages loop
             Sum := Sum + B5(I) * Factors (I, Order-1);
          end loop;

          --  Calculate 1.0 / Order! to get
          Factorial := 1.0;
          for N in Stages range 2..Order loop
             Factorial := Factorial * Real(N);
          end loop;
          Factorial := 1.0 / Factorial;

          if Print_to_Screen_Desired then
             new_line(2);
             put ("           "); put (Real'Image (Sum));
             new_line;
             put ("Should be: "); put (Real'Image (Factorial));
             new_line;
          end if;

          Error := Abs (Sum - Factorial);
          if Error > 1.0e-14 then
             put ("Problem with the Runge Kutta Coeffs in runge_coeffs_5 (1).");
             --raise Program_Error;
          end if;
       
       end loop;

       for Order in Stages range 1..4 loop
          Sum := 0.0;
          for I in Stages loop
            Sum := Sum + B4(I) * Factors (I, Order-1);
          end loop;

          --  Calculate 1.0 / Order! to get
          Factorial := 1.0;
          for N in Stages range 2..Order loop
             Factorial := Factorial * Real(N);
          end loop;
          Factorial := 1.0 / Factorial;

          if Print_to_Screen_Desired then
             new_line(2);
             put ("           "); put (Real'Image (Sum));
             new_line;
             put ("Should be: "); put (Real'Image (Factorial));
             new_line;
          end if;

          Error := Abs (Sum - Factorial);
          if Error > 1.0e-14 then
             put ("Problem with the Runge Kutta Coeffs in runge_coeffs_5 (2).");
             --raise Program_Error;
          end if;
       
       end loop;

    end Test_Runge_Coeffs;

end runge_coeffs_pd_5;