#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <unistd.h>
#include <string.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 [-k xidx,yidx] {-g|-l|-b} -f fwhm [range] < input_file > output_file
 or
    thisprog {-g|-l|-b} -f 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 option '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 option '-g' us used, 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 option '-l' is used, 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.

 The option -k followed by a comman-separated pair of integers allows
 to select other columns than the first two for the x and y data. The
 default is xidx=1 and yidx=2, taking the first two data columns of the input. 

 If option 'b' is used, 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, including order, 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.

 There are no Fourier transform methods used to implement the 
 convolutions; the computational time grows as the square of the
 number of recognized input lines. Only the standard mathematical library with
 the sine and exponential functions needs to be linked. The standard
 compilation is
      gcc -std=gnu99 -O -o convol convol.c -lm
 for example.

 Richard J Mathar http://www.strw.leidenuniv.nl/~mathar , 2010-10-06
*********************************************************************/

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

#define CONVOL_MAX_COLS 10

static void usage(char *argv0)
{
	printf("Usage: %s [-k xidx,yidx] {-g|-l|-b} -f fwhm [range] < infile > outfile\n",argv0) ;
	printf("\t%s [-k idx,yidx] {-g|-l|-b} -f 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], const int xyidx[2])
{
	int lcount=0 ;		/* input line count */
	char line[LINE_MAX] ;	/* input line */
	char fmtx[1+4*CONVOL_MAX_COLS] ,	/* space for a maximumof CONVOL_MAX_COLS %le or %*le */
	     fmty[1+4*CONVOL_MAX_COLS] ;

	fmtx[0]='\0' ;
	for(int f =0 ; f < xyidx[0]-1 && f < CONVOL_MAX_COLS; f++)
		strcat(fmtx,"%*le") ;
	strcat(fmtx,"%le") ;
	fmty[0]='\0' ;
	for(int f =0 ; f < xyidx[1]-1 && f < CONVOL_MAX_COLS; f++)
		strcat(fmty,"%*le") ;
	strcat(fmty,"%le") ;

	while( fgets(line,LINE_MAX,stdin) )
	{
		int items =0;
		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,fmtx,&(*xydata)[lcount][0]) ;
		items += sscanf(line,fmty,&(*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> [-k xidx,yidx] {-g|-l|-b} -f fwhm [range] < input_file
   <thisprog> {-g|-l|-b} -f 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 */
	    xyidx[2] = {1,2},
	    doG = 0 ,
	    doL = 0 ,
	    doB = 0 ,
	    oc;
	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.;

	while (  (oc=getopt(argc,argv,"k:f:glb")) != -1 )
        {
                switch(oc)
                {
                case 'g' :
                        doG=1;
			break ;
                case 'l' :
                        doL=1;
			break ;
                case 'b' :
                        doB=1;
			break ;
                case 'f' :
                        fwhm=atof(optarg) ;
			break ;
                case 'k' :
                        xyidx[0] = atoi(optarg) ;
                        xyidx[1] = atoi( 1+index(optarg,',') ) ; /* skip comma */
			break ;
                case '?' :
                        fprintf(stderr,"Invalid command line option %c\n",optopt) ;
                        return 1 ;
                        break ;
                }
        }

	if ( xyidx[0] > CONVOL_MAX_COLS || xyidx[1] > CONVOL_MAX_COLS)
	{
		fprintf(stderr,"Cannot handle more than %d columns per line\n",CONVOL_MAX_COLS) ;
		return 1 ;
	}

	
	if( argc < 2)
	{
		usage(argv[0]) ;
		return 1;
	}
	if( optind == argc -1)
		range=atof(argv[optind]) ;

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

	if( optind < argc - 1 )
	{
		xrange[0] = atof(argv[optind]) ;
		xrange[1] = atof(argv[optind+1]) ;
		N = atoi(argv[optind+5]) ;
		if( N <2)
		{
			fprintf(stderr,"N= %d must be >=2\n",N) ;
			usage(argv[0]) ;
			return 1;
		}
		xladder=1 ;
	}
	else
	{
		N= n ;
		xladder=0 ;
	}

	if (fwhm <=0.)
	{
		fprintf(stderr,"Fwhm %e must be >0\n",fwhm) ;
		free(xydata) ;
		return 1 ;
	}
#if DEBUG_LEVEL >= 1
	fprintf(stderr,"Using FWHM %e\n",fwhm) ;
#endif
	if ( doG)
		/* Gaussian */
		for(int 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) ;
		}
	else if (doL)
		/* Lorentz profile */
		for(int 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) ;
		}
	else if (doB)
		/* Sinc (boxcar) profile */
		for(int 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) ;
		}
	free(xydata) ;
	return 0 ;
}