/*
**  $Id: gausswts.c,v 1.2 1997/12/22 19:56:53 frank Exp $
**  gausswts.c	1.2	Solaris 2.3 Unix	940906	SIO/ODF	fmd
**
**  Generate a series of Gaussian weights.
*/
#include <stdlib.h>
#include <math.h>
#include <odflib.h>

/*
**  Compute a series of N Gaussian weights with a
**  mean of 1.0 at W[N/2]. The series is malloc'd.
*/
double *				/* <- weights series (malloc'd) */
GaussianWeights(n)
	int	n;			/* -> number of weights to generate */
{
	double	*pW, *pWn, *wts, alpha, wSum, v;
	int	halfN, count;

	if (n <= 0.0)
	    return (0);
	n |= 1;		/* insure n is odd */
	if ((wts = (double *)malloc((unsigned)n*sizeof (double)))==0)
	    return (0);
	halfN = n>>1;
	alpha = ((M_PI*M_PI)/(n*n))/4.5;
	for (wSum=0.0, pW=wts, count=halfN, v = count;
	    count; count--, v--, pW++) {
	    *pW   = exp(-alpha*v*v);
	    wSum += *pW;
	}
	*pW   = 1.0;
	wSum += wSum + 1.0;
	wSum  = 1.0/wSum;
	for (pW=wts, count=halfN+1; count; count--, *pW++ *= wSum);
	for (pWn=pW, pW-=2, count=halfN; count; *pWn++ = *pW--, count--);
	return (wts);
}