with Ada.Text_Io; use Ada.Text_IO ;
with Ada.Long_Float_Text_IO ; use Ada.Long_Float_Text_IO ;
with Ada.Command_Line; use Ada.Command_Line;

with Interfaces.C; use Interfaces.C ;
with Interfaces.C.Strings; use Interfaces.C.Strings ;
with gsl ;
with gsl.rng ;
with gsl.randist ;

procedure bigauss is

   rng : access gsl.rng.gsl_rng ; 

   sigma_x : double := 1.0 ;
   sigma_y : double := 1.0 ;
   rho : double := 0.9 ;

   outfile : File_Type ;

   procedure Create(fn : String ) is
   begin
      Create(outfile, Out_File, fn );
      Set_Output(outfile);
   end Create ;
   procedure Close is
   begin
      Close(outfile);
      Set_Output(Standard_Output);
   end Close;

   procedure TestBivariateGaussian is
      x , y : aliased double ;
   begin 
      Put(Interfaces.C.Strings.Value( gsl.rng.name(rng) ));
      New_Line;
      Create (Value(gsl.rng.name(rng)) & ".bi.txt");
      for i in 1..1024
      loop
         gsl.randist.bivariate_gaussian(rng,sigma_x, sigma_y , rho , x'access , y'access );
         Put(Long_Float(x)); Put(" ; "); Put(Long_Float(y)); 
         New_Line;
      end loop ;
      Close ;
   end TestBivariateGaussian ;

   procedure TestBivariateGaussianPdf is
      x : double := -2.0 ;
      xdelta : double := 2.0 * abs(x) / 1024.0 ;
      y : double := -2.0 ;
      ydelta : double := 2.0 * abs(x) / 1024.0 ;

      pd : double ;
   begin 
      Put(Interfaces.C.Strings.Value( gsl.rng.name(rng) ));
      New_Line;
      Create (Value(gsl.rng.name(rng)) & ".bi.pdf.txt");
      for i in 1..1024
      loop
         y := -2.0 ;
         for j in 1..1024
         loop
            pd := gsl.randist.bivariate_gaussian_pdf(x,y,sigma_x, sigma_y , rho );
            Put(Long_Float(x)); Put(" ; ");
            Put(Long_Float(y)); Put(" ; ");
            Put(Long_Float(pd));
            New_Line;
            y := y + ydelta ;
         end loop ;
         x := x + xdelta ;
      end loop ;
      Close ;
   end TestBivariateGaussianPdf ;

begin
   Put_Line("Default Random Number Generator ");
   Put_Line("Setup from environment");
   rng := gsl.rng.alloc(gsl.rng.env_setup);
   if Argument_Count > 0
   then
      sigma_x := double'Value(Argument(1));
      if Argument_Count > 1
      then
         sigma_y := double'Value(Argument(2));
         if Argument_Count > 2
         then
            rho := double'Value(Argument(3));
         end if ;
      end if ;
   end if;

   TestBivariateGaussian;
   TestBivariateGaussianPdf;
   gsl.rng.free(rng);
end bigauss ;