--
-- By Nasser M. Abbasi   sept 25, 2010
-- Math 228A, UC DAVIS, Fall 2010
--
-- to help me learn Ada
--
-- gnatmake -P t1.gpr
--
--  $ ./t1.exe
--      h           D+         D_           D0              D3
--    1.0E-01 -4.2939E-02  4.1138E-02 -9.0005E-04  6.8207E-05
--    5.0E-02 -2.1257E-02  2.0807E-02 -2.2510E-04  8.6491E-06
--    1.0E-02 -4.2163E-03  4.1983E-03 -9.0050E-06  6.9941E-08
--    5.0E-03 -2.1059E-03  2.1014E-03 -2.2513E-06  8.7540E-09
--    1.0E-03 -4.2083E-04  4.2065E-04 -9.0050E-08  6.9979E-11

with ada.text_io;              use ada.text_io;
with Ada.Characters.Latin_1;   use Ada.Characters.Latin_1;
with Ada.Numerics.Long_Elementary_Functions; use
  Ada.Numerics.Long_Elementary_Functions;

procedure t1 is
   package my_float_IO is new Ada.Text_IO.float_IO(long_float);
   use  my_float_IO;

   h : constant array(1..5) of long_float := (0.1,0.05,0.01,0.005,0.001);
   result : array(1..4) of long_float;
   exact_value : constant long_float := cos(1.0);

   --------------------------------------------------------------
   function derivative_left(h: in long_float) return long_float is
   begin
      return ( sin(1.0) - sin(1.0 - h)  ) / h;
   end derivative_left;

   --------------------------------------------------------------
   function derivative_right(h: in long_float) return long_float is
   begin
      return ( sin(1.0+h) - sin(1.0) ) / h;
   end derivative_right;

   --------------------------------------------------------------
   function derivative_center(h: in long_float) return long_float is
   begin
      return ( sin(1.0+h) - sin(1.0-h) ) / (2.0*h);
   end derivative_center;

   --------------------------------------------------------------
   function  derivative_3rd_order(h: in long_float) return long_float is
   begin
      return ( 2.0*sin(1.0+h) + 3.0*sin(1.0)-6.0*sin(1.0-h) + sin(1.0-2.0*h) )
        / (6.0*h);
   end derivative_3rd_order;
   --------------------------------------------------------------

begin
   put_line("    h"& HT&HT& "D+" & HT& "   D_" & HT&HT&"D0" & HT&HT&"D3");

   for I in h'range loop
      result(1) := derivative_right(h(I))     -  exact_value;
      result(2) := derivative_left(h(I))      -  exact_value;
      result(3) := derivative_center(h(I))    -  exact_value;
      result(4) := derivative_3rd_order(h(I)) -  exact_value;

      put(item=>h(I),Fore=>3,Aft=>1);
      for J in result'range loop
         put(item=>result(j),Fore=>3,Aft=>4);
      end loop;

      new_line;

   end loop;

end t1;