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.interp ;
with gsl.spline ;

-- 30.7 1D Interpolation Example Programs

procedure Interp is
   use gsl.double_text_io;

   Verbose : boolean := True ;
   Status : Int ;
   procedure Test1 is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 10 ;
      x,y : gsl.double_array(1..N);
      acc : access gsl.interp.gsl_interp_accel := gsl.interp.accel_alloc ;
      spline : access gsl.spline.gsl_spline := gsl.spline.alloc(gsl.interp.gsl_interp_akima,N);
      xi, yi,xdelta : double ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;

      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);

      for i in 1..N
      loop
         x(i) := double(i) + double( 0.5 * Sin(Float(i))) ;
         y(i) := double(i) + double( Cos(Float(i)) * Cos(Float(i))) ;
         Put(x(i)) ; Put(" ; ");
         Put(y(i)) ;
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

      Status := gsl.spline.init(spline,x,y,N);
      Create(logfile,Out_File,myname & ".interp.csv");
      Set_Output(logfile);
      xi := x(x'First) ;
      xdelta := 0.1 ;
      loop
         if xi >= x(x'Last)
         then
            exit ;
         end if ;
         yi := gsl.spline.eval(spline,xi,acc) ;
         Put(xi) ; Put(" ; ");
         Put(yi) ;
         New_Line ;

         xi := xi + xdelta ;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

   end Test1 ;

   procedure Test2 is
   -- a periodic cubic spline with 4 data points
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
      N : constant := 4 ;
      NP : constant := 100 ;
      x : gsl.double_array := (0.00, 0.10,  0.27,  0.30) ;
      y : gsl.double_array := (0.15, 0.70, -0.10,  0.15) ;
      xi, yi : double ;
      xdelta : double := 0.1 ;
      acc : access gsl.interp.gsl_interp_accel := gsl.interp.accel_alloc ;
      spline : access gsl.spline.gsl_spline := 
               gsl.spline.alloc(gsl.interp.gsl_interp_cspline_periodic,N);
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      for i in 1..N
      loop
         Put(x(i-1)) ; Put(" ; ");
         Put(y(i-1)) ;
         New_Line;
      end loop ;

      Set_Output(Standard_Output);
      Close(logfile);

      Status := gsl.spline.init(spline,x,y,N);
      Create(logfile,Out_File,myname & ".interp.csv");
      Set_Output(logfile);
      xi := x(x'First) ;
      xdelta := (x(x'Last) - x(x'First)) / double(NP) ;
      for i in 1..NP
      loop
         if xi >= x(x'Last)
         then
            exit ;
         end if ;
         yi := gsl.spline.eval(spline,xi,acc) ;
         Put(xi) ; Put(" ; ");
         Put(yi) ;
         New_Line ;
         xi := xi + xdelta ;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);
   end Test2 ;

   procedure Test3 is
   -- Different techniques of interpolation on difficult dataset
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;


-- this dataset is taken from
-- J. M. Hyman, Accurate Monotonicity preserving cubic interpolation, 
--         * SIAM J. Sci. Stat. Comput. 4, 4, 1983. 

      x : gsl.double_array := ( 7.99, 8.09, 8.19, 8.7, 9.2,
                           10.0, 12.0, 15.0, 20.0 ) ;
      y : gsl.double_array := ( 0.0, 2.76429e-5, 4.37498e-2,
                           0.169183, 0.469428, 0.943740,
                           0.998636, 0.999919, 0.999994 );
      N : size_t := size_t(x'Length) ;
      acc : access gsl.interp.gsl_interp_accel := gsl.interp.accel_alloc ;
      cubic_spline : access gsl.spline.gsl_spline := 
                  gsl.spline.alloc(gsl.interp.gsl_interp_cspline,N);
      akima_spline : access gsl.spline.gsl_spline := 
                  gsl.spline.alloc(gsl.interp.gsl_interp_akima,N);
      steffen_spline : access gsl.spline.gsl_spline := 
                  gsl.spline.alloc(gsl.interp.gsl_interp_steffen,N);
      NP : constant := 100 ;
      xi, yi : double ;
      xdelta : double ;

   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);

      for i in 1..Integer(N)
      loop
         Put(x(i-1)) ; Put(" ; ");
         Put(y(i-1)) ;
         New_Line;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);
      Status := gsl.spline.init(cubic_spline,x,y,N);
      Status := gsl.spline.init(akima_spline,x,y,N);
      Status := gsl.spline.init(steffen_spline,x,y,N);
      Create(logfile,Out_File,myname & ".interp.csv");
      Set_Output(logfile);
      xi := x(x'First) ;
      xdelta := (x(x'Last) - x(x'First)) / double(NP) ;
      for i in 1..NP
      loop
         if xi >= x(x'Last)
         then
            exit ;
         end if ;

         Put(xi) ; Put(" ; ");
         yi := gsl.spline.eval(cubic_spline,xi,acc) ; Put(yi) ; Put(" ; ");
         yi := gsl.spline.eval(akima_spline,xi,acc) ; Put(yi) ; Put(" ; ");
         yi := gsl.spline.eval(steffen_spline,xi,acc) ; Put(yi) ; Put(" ; ");
         New_Line ;
         xi := xi + xdelta ;
      end loop ;
      Set_Output(Standard_Output);
      Close(logfile);

   end Test3 ;

   procedure Template is
      myname : String := gnat.Source_Info.enclosing_entity ;
      logfile : File_Type ;
   begin
      if Verbose
      then
         Put_Line(myname);
      end if ;
      Create(logfile,Out_File,myname & ".csv");
      Set_Output(logfile);
      Set_Output(Standard_Output);
      Close(logfile);
   end Template ;

begin
   Test1 ;
   Test2 ;
   Test3 ;
end Interp;