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

/*********************************************************************
 This is a simple filtering/convolution/smoothing program for a
 one-dimensional sequence of (x,y) input data, which alternatively uses
   - Gaussian
   - Lorentz
   - Sinc
 curves for the convolution kernel ("instrument funtion").
 The syntax of the UNIX call is
    thisprog {g|l|b} fwhm [range] < input_file > output_file
 or
    thisprog {g|l|b} fwhm xstart xend N [range] < input_file > output_file

 where each line of the input ASCII ("formatted") file contains either
 a comment line that starts with the character '#' and is ignored, or
 a line that contains one x and one y coordinate in floating point format
 as described under sscanf(3s), separated by white space. Any additional
 bytes in a line after those two coordinates is also ignored.

 The second command line argument 'fwhm' must be a floating point number
 larger than zero, which is used as the full width at half maximum of
 the instrument function.

 If the first argument of the command line is 'g', a Gaussian profile
 weight function of the form ~exp{-Delta^2/(2*sigma^2)} with
 sigma=fwhm/(2*sqrt(2*log(2))) is used as the filter profile.
 The convolution sum (see below) is limited to |Delta|<range*fwhm in this
 case.

 If the first argument of the command line is 'l', a Lorentzian profile
 weight function of the form ~1/[Delta^2+fwhm^2/4] is used as the
 filter profile.
 The convolution sum (see below) is limited to |Delta|<range*fwhm in this
 case.

 If the first argument of the command line is 'b', a profile
 weight function of the form ~sin(Delta*fwhm/2)/Delta is used as the
 filter profile, which would become a boxcar filter in the Fourier
 domain.
 The convolution sum (see below) is limited to |Delta|<range*fwhm in this
 case.

 In all of these cases 'Delta' is the distance between the 
 current point x of interest in the output file and the point in the
 original data file. The calculated output value y is calculated
 via a discrete convolution
 y(x) = sum_over_all_i[weight(x,x_i,fwhm)*y_i]
        /sum_over_all_i[weight(x,x_i,fwhm)]

 The last command line parameter 'range' is optional and
 is set to 3. if not provided.

 The standard output (usually redirected to another ASCII file)
 contains those re-calculated pairs x and y. The set of
 x is the same as provided by the input if the first syntax is
 used to call the program, and a set of 'N' equidistant abscissa
 values from 'xstart' to 'xend' if the second syntax is used.
 'N' is an integer >=2 .
 
 The standard output does not contain the comment lines of the input
 or any characters that were following the original pair of input
 data.
.
 Richard J Mathar, 18 Jan 2002
*********************************************************************/

/* could be 0,1,2,...4,5 */
#define DEBUG_LEVEL 0

static void usage(char *argv0)
{
	printf("Usage: %s {g|l|b} fwhm [range] < infile > outfile\n",argv0) ;
	printf("\t%s {g|l|b} fwhm xstart xend N [range] < infile > outfile\n",argv0) ;
}

/* Read standard input with two ASCII doubles separated by white space per line.
Ingnore input lines that start with a #...
Return number of valid non-comment lines.
*/
static int readstdin(double (**xydata)[2])
{
	int lcount=0 ;		/* input line count */
	char line[LINE_MAX] ;	/* input line */

	while( fgets(line,LINE_MAX,stdin) )
	{
		int items ;
		if( line[0] == '#')
			continue ;
		*xydata= (double(*)[2])realloc(*xydata,(lcount+1)*sizeof(double[2])) ;
		if( xydata == NULL)
		{
			fprintf(stderr,"insufficient memory for %d bytes\n",(lcount+1)*sizeof(double[2]) ) ;
			return 0 ;
		}
		items=sscanf(line,"%le%le",&(*xydata)[lcount][0],&(*xydata)[lcount][1]) ;
#if DEBUG_LEVEL >= 4
		fprintf(stderr,"%d input %e %e\n",lcount,(*xydata)[lcount][0],(*xydata)[lcount][1]) ;
#endif
		if( items == 2)
			lcount ++ ;
	}
#if DEBUG_LEVEL >= 3
	fprintf(stderr,"input line count %d\n",lcount) ;
#endif
	return lcount ;
}

static double gaussian(const double x, const double (* const xydata)[2], const int n, const double fwhm, const double range)
{
	const double sigma = fwhm/(2.*sqrt(2.*M_LN2)) ;
	const double denom = 2.*sigma*sigma ;
	double result=0. ;
	double weight=0. ;
	int i ;
#if DEBUG_LEVEL >= 3
	fprintf(stderr,"Gaussian for %e with %d inputs, sigma=%e\n",x,n,sigma) ;
#endif
	for(i=0; i<n; i++)
	{
		const double dist = x- xydata[i][0] ;
#if DEBUG_LEVEL >= 4
		fprintf(stderr,"Computing %d for %e\n",i,x) ;
#endif
		/* limit output to a scan across range * fwhm */
		if( fabs(dist) < range*fwhm)
		{
			const double tmp=exp(-dist*dist/denom) ;
			weight += tmp ;
			result += tmp*xydata[i][1] ;
		}
#if DEBUG_LEVEL >= 4
		fprintf(stderr,"Computing %d for %e: dist %e, weight %e, result %e\n",i,x,dist,weight,result) ;
#endif
	}
	if (weight)
		return result/weight ;
	else
		return 0. ;
}

static double lorentz(const double x, const double (* const xydata)[2], const int n, const double fwhm,const double range)
{
	const double eps2 = 0.25*fwhm*fwhm ;
	double result=0. ;
	double weight=0. ;
	int i ;
#if DEBUG_LEVEL >= 3
	fprintf(stderr,"Lorentzian for %e with %d inputs, fwhm=%e\n",x,n,fwhm) ;
#endif
	for(i=0; i<n; i++)
	{
		const double dist = x- xydata[i][0] ;
#if DEBUG_LEVEL >= 4
		fprintf(stderr,"Computing %d for %e\n",i,x) ;
#endif
		/* limit output to scanning range * fwhm */
		if( fabs(dist) < range*fwhm)
		{
			const double tmp=1./(dist*dist+eps2) ;
			weight += tmp ;
			result += tmp*xydata[i][1] ;
		}
#if DEBUG_LEVEL >= 4
		fprintf(stderr,"Computing %d for %e: dist %e, weight %e, result %e\n",i,x,dist,weight,result) ;
#endif
	}
	if (weight)
		return result/weight ;
	else
		return 0. ;
}

static double sinc(const double x, const double (* const xydata)[2], const int n, const double fwhm, const double range)
{
	double result=0. ;
	double weight=0. ;
	int i ;
#if DEBUG_LEVEL >= 3
	fprintf(stderr,"Boxcar for %e with %d inputs, fwhm=%e\n",x,n,fwhm) ;
#endif
	for(i=0; i<n; i++)
	{
		const double dist = x- xydata[i][0] ;
#if DEBUG_LEVEL >= 4
		fprintf(stderr,"Computing %d for %e\n",i,x) ;
#endif
		/* limit output to scanning range * fwhm */
		if( fabs(dist) < range*fwhm)
		{
			const double tmp=(x ) ? sin(0.5*fwhm*x)/x : 0.5*fwhm ;
			weight += tmp ;
			result += tmp*xydata[i][1] ;
		}
#if DEBUG_LEVEL >= 4
		fprintf(stderr,"Computing %d for %e: dist %e, weight %e, result %e\n",i,x,dist,weight,result) ;
#endif
	}
	if (weight)
		return result/weight ;
	else
		return 0. ;
}

/* UNIX synopsis:
   <thisprog> {g|l|b} fwhm [range] < input_file
   <thisprog> {g|l|b} fwhm xstart xend N [range] < input_file
*/
int main(int argc, char *argv[])
{
	double (*xydata)[2] = NULL ;		/* xydata[i][0] contains the x and xydata[i][1] the y input value */
	int n=0 ,				/* number of input values and size of xydata[i=0..n-1][] */
	    N =2,				/* subdivisions for 2nd command line format */
	    xladder=1,				/* nonzero for 2nd command line format, zero for first one */
	    i;
	double fwhm = -1.,			/* Full width at half maximum of weight function */
	       xrange[2] ,			/* output range of x values for 2nd command line format */
	       range = 3.;
	
	if( argc < 3)
	{
		usage(argv[0]) ;
		exit(EXIT_FAILURE) ;
	}
	if( argc == 4 || argc == 7)
		range=atof(argv[argc-1]) ;

	n=readstdin(&xydata) ;
	if( n==0)
	{
		fprintf(stderr,"No valid lines in stdin\n") ;
		free(xydata) ;
		exit(EXIT_FAILURE) ;
	}

	if( argc >= 6 )
	{
		xrange[0] = atof(argv[3]) ;
		xrange[1] = atof(argv[4]) ;
		N = atoi(argv[5]) ;
		if( N <2)
		{
			fprintf(stderr,"N= %d must be >=2\n",N) ;
			usage(argv[0]) ;
			exit(EXIT_FAILURE) ;
		}
		xladder=1 ;
	}
	else
	{
		N= n ;
		xladder=0 ;
	}

	fwhm=atof(argv[2]) ;
	if (fwhm <=0.)
	{
		fprintf(stderr,"Fwhm %e must be >0\n",fwhm) ;
		free(xydata) ;
		exit(EXIT_FAILURE) ;
	}
#if DEBUG_LEVEL >= 1
	fprintf(stderr,"Using FWHM %e\n",fwhm) ;
#endif
	switch (argv[1][0])
	{
	int i ;
	case 'g':		/* Gaussian */
		for(i=0 ; i<N ; i++)
		{
			const double x= (xladder) ? (xrange[0]+i*(xrange[1]-xrange[0])/(double)N)  : xydata[i][0] ;
			const double resul =gaussian(x,xydata,n,fwhm,range) ;
			printf("%e %e\n",x, resul) ;
		}
		break ;
	case 'l':		/* Lorentz profile */
		for(i=0 ; i<N ; i++)
		{
			const double x= (xladder) ? (xrange[0]+i*(xrange[1]-xrange[0])/(double)N)  : xydata[i][0] ;
			const double resul =lorentz(x,xydata,n,fwhm,range) ;
			printf("%e %e\n",x, resul) ;
		}
		break ;
	case 'b':		/* Sinc (boxcar) profile */
		for(i=0 ; i<N ; i++)
		{
			const double x= (xladder) ? (xrange[0]+i*(xrange[1]-xrange[0])/(double)N)  : xydata[i][0] ;
			const double resul =sinc(x,xydata,n,fwhm,range) ;
			printf("%e %e\n",x, resul) ;
		}
		break ;
	default :		/* Gaussian */
		fprintf(stderr,"Unrecognized filter type %c\n",argv[1][0]) ;
		break ;
	}
	free(xydata) ;
	return 0 ;
}