#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define MAX_N 100000
#define MAX_BINS 100


// function prototipes
double rgauss(void);
int fillData(double data[], double, double);
void makeHist(double *data, int *hist, double min, double max, int nbins, int nData);
void writeHist(int *hist, int nbins, double min, double max);
void writeGauss(int nbins, double min, double max, double mu, double sigma);


int main()
{

  /* Generate numbers distributed with Gaussian distribution with mean mu and variance sigma. 
     Make histogram.  
  */


  unsigned long int seed;
  double data[MAX_N];
  int hist[MAX_BINS]={0};
  int N,i,nbins;
  double mu,sigma;
  double min,max;
  FILE *fp;

  seed = time(0);
  srand48(seed);

  printf("Insert mean of Gaussian\n");
  scanf("%lf",&mu);
  do { 
    printf("Insert sigma of Gaussian\n");
    scanf("%lf",&sigma);
  } while (sigma <= 0);


  N=fillData(data,mu,sigma); // Fill vector data with random numbers 

  min = mu-4.*sigma;
  max = mu+4.*sigma;
 

 do { 
    printf("Number of bins for histogram?\n");
    scanf("%d",&nbins);
  } while (nbins <= 0 || nbins > MAX_BINS);

  makeHist(data, hist, min, max, nbins, N); 
  writeHist(hist,nbins,min,max);
  writeGauss(nbins,min,max,mu,sigma); // write analytical gauss. function


  return 0;
}


int fillData(double data[], double mu, double sigma) {
  int k,N;
  double g;
  do {
    printf("Number of data points?\n");
    scanf("%d",&N);
  } while(N<0 || N>MAX_N);
  for (k=0;k<N;k++) {
    g=rgauss()*sigma+mu; // Gaussian distriuted number with mean mu and std sigma
    data[k]=g;
  }
  return N;
}



double rgauss(void) {
  /* How to generate a Gaussian dist. 
     1) generate two random numbers r1 and r2 with uniform dist. between 0 and 1
     2) compute gamma=sin(2.*M_PI * r1) * sqrt(-2.*log(r2)); ; gamma is Gaussian with mu=0 and sigma=1
  */
  double g;
  double r1,r2;
  r1=(double)lrand48()/(RAND_MAX+1.);
	do{r2=(double)lrand48()/(RAND_MAX+1.);
	}while(r2==0.);
  g= sin(2.*M_PI * r1) * sqrt(-2.*log(r2));
  return g;
}


void makeHist(double *data, int *hist, double min, double max, int nbins, int nData) {   
  int i,j;
  double x;
  double invWidth = (double)nbins / (max-min);
  for (i=0; i < nData; i++) {
    x = *(data+i);
    j =  (x - min) * invWidth;
    if (j < nbins ) (*(hist+j))++;
  }
}
   

void writeHist(int *hist, int nbins, double min, double max){

  int j;
  FILE *fp;
  double x,y;
  double norm=0.;
  double width;

  width=(max-min) / (double)nbins;

  // Normalization for histogram
  for (j=0; j < nbins; j++) {
    norm = norm + *(hist+j);
  }
  norm*=width;
  norm=1./norm;

  if((fp = fopen("hist.dat", "w")) == NULL) {
    printf("Cannot open file.\n");
    exit(1);
  } 

  for (j=0; j < nbins; j++) {
    x = (j+0.5)*width + min;    
    y = (double)hist[j];
    fprintf(fp,"%.4f %f\n",x,y*norm);    
  }

  if( fclose( fp )) 
    printf("File close error.\n");
}      


void writeGauss(int nbins, double min, double max, double mu, double sigma){

  // Analitycal expression for Gaussian computed in the center of bins

  int j;
  FILE *fp;
  double x,y;
  double sigma2=2.*pow(sigma,2);
  double norm = 1.0/sqrt(M_PI * sigma2);
  double width;

  width=(max-min) / (double)nbins;

  if((fp = fopen("gauss.dat", "w")) == NULL) {
    printf("Cannot open file.\n");
    exit(1);
  }  

  for (j=0; j < nbins; j++) {
    x = (j+0.5)*width + min;
    y = norm*exp(-pow(x-mu,2)/sigma2);    
    fprintf(fp,"%.4f %f \n",x,y);    
  }

  if( fclose( fp )) 
    printf("File close error.\n");

}