///////////////////////////////////////////////////////////////////////////
//
// PURPOSE: Test the Gaussian random number generator.
// AUTHORS: 1. Rob McEwen, 25 Aug 2016
//          2.
//
///////////////////////////////////////////////////////////////////////////
//
#include"MathP.h"
#include <stdio.h>

void main()
{
   const int Nsamp = 1000;
   int i;
   long idum[1] = {-100};
   float gaussian[Nsamp];
   float uniform[Nsamp];
   double meanG = 0.;
   double ssqG  = 0.;
   double meanU = 0.;
   double ssqU  = 0.;
   double maxU  = 0.;
   
   for( i=0; i<Nsamp; i++ )
   {
      gaussian[i] = Math::gasdev( idum );
      meanG += gaussian[i];
      //
      // 1/sqrt(12) is the s.d. for a unit interval, for example, (-.5,.5).
      // 1/sqrt(3)  is the s.d. for 2*unit interval, for example, (-1,1).
      uniform[i] = sqrt(12.)*(Math::ran1( idum ) -.5) ;
      meanU += uniform[i];
      if( fabs( uniform[i] ) > maxU ) maxU = fabs( uniform[i] );

      printf(" %15.4f   %15.6f\n", gaussian[i], uniform[i] );
   }

   meanG = meanG/Nsamp;
   printf(" Gaussian Mean = %15.6f\n", meanG );

   meanU = meanU/Nsamp;
   printf(" Uniform Mean  = %15.6f\n", meanU );

   for( i=0; i<Nsamp; i++ )
   {
      ssqG += pow( ((double)gaussian[i] - meanG), 2.);
      ssqU += pow( ((double)uniform[i]  - meanU), 2.);
   }
   printf(" Gaussian Standard deviation = %15.6f\n", sqrt( (double)1/(Nsamp-1) * ssqG) );
   printf(" Uniform Standard deviation  = %15.6f\n", sqrt( (double)1/(Nsamp-1) * ssqU) );
   printf(" Max of Uniform R.V.         = %15.6f\n", maxU/sqrt(12.) );
   //printf(" One over square root of 12  = %15.6f\n", 1./sqrt(12));
   
}