#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>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_math.h>

#define PI 3.1415926

int main(void)
{
	int i,Ndata,Ndata1,Ndata2,nbin;
	double *x,*x1,*x2,*bin,*fftbin,*t, *gaus, h;
	double y,logL;
	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);

	Ndata1=Ndata2=Ndata/2;
	x1=malloc(Ndata/2*sizeof(double));
	x2=malloc(Ndata/2*sizeof(double));
	for(i=0;i<Ndata/2;i++)
	{
		x1[i]=x[i*2];
		x2[i]=x[i*2+1];
		//printf("%g %g\n",x1[i],x2[i]);
	}

	nbin = 100;
	bin=malloc(nbin*sizeof(double));
	fftbin=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<Ndata1;i++)
	{
	        gsl_histogram_increment(hist,x1[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);

	gaus=malloc(nbin*sizeof(double));
	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);
	hc = gsl_fft_halfcomplex_wavetable_alloc(nbin);
	gsl_interp_accel *acc = gsl_interp_accel_alloc();
	gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline,nbin);
	// FFT data
	gsl_fft_real_transform(bin, 1, nbin, real, work);
	for(i=0;i<nbin;i++) fftbin[i]=bin[i];

	h=0.36;
	do{
	
	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);
	}
	
	// product
	for(i=0;i<nbin;i++) bin[i]=fftbin[i]*gaus[i];

	// IFFT
	gsl_fft_halfcomplex_inverse(bin,1,nbin,hc,work);
	for(i=0,tot=0;i<nbin;i++) 
	{
		//printf("%g %g\n",t[i],bin[i]);
		tot+=bin[i];
	}
	//printf("%g\n",tot);
	
	gsl_spline_init(spline,t,bin,nbin);
	for(i=0,logL=0.0;i<Ndata2; i++)
	{
		y = gsl_spline_eval(spline, x2[i], acc);
		//printf("%g %g\n",x2[i],y);
		logL += log10(gsl_max(y,0));
	}
	//printf("\n");
	printf("%g %g\n",h,logL);
	h+=0.0001;
	}
	while(h<=0.38);
	gsl_fft_real_wavetable_free(real);
	gsl_fft_halfcomplex_wavetable_free(hc);
	gsl_fft_real_workspace_free(work);
	gsl_spline_free(spline);
	gsl_interp_accel_free(acc);

	return 1;
}