with System;
with Ada.Text_IO;            use Ada.Text_IO;
with Ada.Text_IO.C_Streams;
with Ada.Long_Float_Text_IO; use Ada.Long_Float_Text_IO;
with Ada.Integer_Text_IO;    use Ada.Integer_Text_IO;
with Interfaces.C;           use Interfaces.C;
with Interfaces.C.Strings;   use Interfaces.C.Strings;
with Ada.Strings.Unbounded;  use Ada.Strings.Unbounded;
with Ada.Command_Line;       use Ada.Command_Line;
with Ada.Numerics.Elementary_Functions ; use Ada.Numerics.Elementary_Functions ;

with GNAT.Source_Info; use GNAT.Source_Info;
with gsl;
with gsl.math ;
with gsl.integration ;
with gsl.sf.gamma ;

with fns ;

procedure integration is
   use gsl.double_text_Io;
   use gsl.double_elementary_functions ;
   Verbose : boolean := True ;
   Status : Int ;

-- 17.16.1 Adaptive integration example
   procedure Test1 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 1000 ;
      ws : access gsl.integration.gsl_integration_workspace 
            := gsl.integration.workspace_alloc(N) ;
      
      result, error : aliased double ;
      expected : constant double := -4.0;
      alpha : aliased double := 1.0;

      fn : aliased gsl.math.gsl_function ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      fn.c_function := fns.LogSqrt'Access ;
      fn.params := alpha'Address ;
      Status := gsl.integration.qags( fn'Access , 0.0 , 1.0 , 0.0 , 1.0e-7 , 1000, ws, result'access , error'access);
      Put("Result "); Put(result, aft => 6 ) ; Put(" ");
      Put("Rel Error "); Put(error) ; Put(" ");
      Put("Intervals "); Put(Integer(ws.size));
      New_Line;
      -- Create(logfile,Out_File,myname & ".csv");
      -- Set_Output(logfile);
      -- Set_Output(Standard_Output);
      -- Close(logfile);
      gsl.integration.workspace_free(ws);
   end Test1 ;

-- 17.16.2 Fixed-point quadrature example
   procedure Test2 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 6 ;
      M : aliased int := 10 ;
      ws : access gsl.integration.gsl_integration_fixed_workspace ;
      fn : aliased gsl.math.gsl_function ;
      result : aliased double ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      ws := gsl.integration.fixed_alloc(gsl.integration.gsl_integration_fixed_hermite, 
                     N , 0.0 , 1.0 , 0.0 , 0.0 );
      fn.c_function := fns.Power'access ;
      fn.params := M'Address ;
      Status := gsl.integration.fixed(fn'access, result'access, ws ) ;
      Put("Result "); Put(result); New_Line ;
      Put("Expected ");
      if m mod 2 = 0
      then
         Put(gsl.math.M_SQRTPI + gsl.sf.gamma.gamma(0.5*(1.0 + double(m))));
      else
         Put(double(gsl.math.M_SQRTPI));
      end if ;
      New_Line;
      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      Set_Output(Standard_Output);
      Close(logfile);
      gsl.integration.fixed_free(ws);
   end Test2 ;  

   procedure Test3 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 6 ;
      ws : access gsl.integration.gsl_integration_fixed_workspace ;
      fn : aliased gsl.math.gsl_function ;
      result : aliased double ;
      x : double ;
      xdelta : constant double := 2.0 * Ada.Numerics.Pi/256.0 ;
      i : Integer := 1 ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;


      fn.c_function := fns.Sin'access ;
      fn.params := System.Null_Address ;

      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      x := xdelta;
      i := 0 ;
      while i < 256
      loop
         i := i + 1 ;
         ws := gsl.integration.fixed_alloc(gsl.integration.gsl_integration_fixed_legendre, 
                     N , 0.0 , x , 0.0 , 0.0 );
         Put(i) ; Put(" ; ");
         Put(x) ; Put(" ; ");
         Put(fns.Sin(x,System.Null_Address)); Put(" ; ");
         Status := gsl.integration.fixed(fn'access , result'access , ws);
         Put(result) ; Put(" ; ");
         New_Line ;
         x := x + xdelta ;
         gsl.integration.fixed_free(ws);
      end loop ;

      Set_Output(Standard_Output);
      Close(logfile);
  end Test3 ;

begin
   if Verbose
   then
      Put_line(Gnat.Source_Info.enclosing_entity);
   end if;
   gsl.double_text_Io.Default_exp := 0;

   Test1 ;
   Test2 ;
   Test3 ;

end integration;