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

#define MAX_N 1000000
#define MAX_BINS 100


// function prototipes
double ExpDist(void);
int fillData(double data[], 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 writeExp(int nbins, double min, double max, double lambda);


int main()
{

  /* Generate numbers distributed with exponential distribution with mean life lambda.
   Make histogram.  
  */


  unsigned long int seed;
  double data[MAX_N];
  int hist[MAX_BINS]={0};
  int N,i,nbins;
  double lambda;
  double min,max;

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


  do { 
    printf("Inserire la vita media \n");
    scanf("%lf",&lambda);
  } while (lambda <= 0);


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

  min = 0.;
  max = 15.*lambda;
 

 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);
  writeExp(nbins,min,max,lambda); // write analytical distribution function


  return 0;
}


int fillData(double data[], double lambda) {
  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=ExpDist()/lambda; // exponential distribution with "Mean life" lambda
    data[k]=g;
  }
  return N;
}



double ExpDist(void) {
  
  double v;	
  double r;
	
    do{r=lrand48()/(RAND_MAX+1.);
	}while(r==0.);
	v=-log(r);  

  return v;
}


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("histExp.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 writeExp(int nbins, double min, double max, double lambda){

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

  int j;
  FILE *fp;
  double x,y;
  double width;

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

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

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

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

}