#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

#define PI 3.1415926

int main(void)
{
	int i,Ndata,nbin;
	double *x,*bin,*t, *gaus, h;
	double xmin=1.e30,xmax=-1.e30, higher,lower,binsize,tot;
	FILE *fp;

	Ndata=40;
	x=malloc(Ndata*sizeof(double));
	fp=fopen("data.dat","r");
	for(i=0;i<Ndata;i++) 
	{
		fscanf(fp,"%lf",&x[i]);
		xmin=(x[i]<xmin) ? x[i] : xmin;
		xmax=(x[i]>xmax) ? x[i] : xmax;
		//printf("%g\n",x[i]);
	}
	fclose(fp);

	nbin = 100;
	bin=malloc(nbin*sizeof(double));
	t=malloc(nbin*sizeof(double));
	xmin=0.0; xmax=10.0;
	binsize=(xmax-xmin)/(float) nbin;
	gsl_histogram * hist = gsl_histogram_alloc(nbin);
	gsl_histogram_set_ranges_uniform(hist,xmin,xmax);
	for(i=0;i<Ndata;i++)
	{
	        gsl_histogram_increment(hist,x[i]);
	}
	//gsl_histogram_fprintf(stdout, hist, "%f", "%f");
	for(i=0;i<nbin;i++) 
	{
		if(!gsl_histogram_get_range(hist,i,&lower,&higher))
		t[i]=0.5*(lower+higher);
		bin[i]=gsl_histogram_get(hist,i);
		//printf("%g %g\n",t[i],bin[i]);
	}
	gsl_histogram_free(hist);

	// FFT window
	h=0.368;
	gaus=malloc(nbin*sizeof(double));
	for(i=0;i<nbin/2;i++) 
	{
		gaus[i*2]=/*pow(2.0*PI,1.5)*h*h*h*/exp(-0.5*i*i*h*h);
		gaus[i*2+1]=/*pow(2.0*PI,1.5)*h*h*h*/exp(-0.5*i*i*h*h);
	}

	// FFT data
	gsl_fft_real_wavetable * real;
	gsl_fft_halfcomplex_wavetable * hc;
	gsl_fft_real_workspace * work;

	work = gsl_fft_real_workspace_alloc(nbin);
	real = gsl_fft_real_wavetable_alloc(nbin);

	gsl_fft_real_transform(bin, 1, nbin, real, work);
	gsl_fft_real_wavetable_free(real);
	//for(i=0;i<nbin;i++) printf("%g %g\n",t[i],bin[i]);
	
	// product
	for(i=0;i<nbin;i++) bin[i]=bin[i]*gaus[i];

	// IFFT
	hc = gsl_fft_halfcomplex_wavetable_alloc(nbin);
	gsl_fft_halfcomplex_inverse(bin,1,nbin,hc,work);
	gsl_fft_halfcomplex_wavetable_free(hc);
	for(i=0,tot=0;i<nbin;i++) 
	{
		printf("%g %g\n",t[i],bin[i]);
		tot+=bin[i];
	}
	//printf("%g\n",tot);
	
	gsl_fft_real_workspace_free(work);
	return 1;
}