#include <math.h>
#define NRANSI
#include "nrutil.h"
#define M 64
#define NDFT 1024
#define MPOL 6
#define TWOPI (2.0*3.14159265)

void dftint(float (*func)(float), float a, float b, float w, float *cosint,
	float *sinint)
{
	void dftcor(float w, float delta, float a, float b, float endpts[],
		float *corre, float *corim, float *corfac);
	void polint(float xa[], float ya[], int n, float x, float *y, float *dy);
	void realft(float data[], unsigned long n, int isign);
	static int init=0;
	int j,nn;
	static float aold = -1.e30,bold = -1.e30,delta,(*funcold)(float);
	static float data[NDFT+1],endpts[9];
	float c,cdft,cerr,corfac,corim,corre,en,s;
	float sdft,serr,*cpol,*spol,*xpol;

	cpol=vector(1,MPOL);
	spol=vector(1,MPOL);
	xpol=vector(1,MPOL);
	if (init != 1 || a != aold || b != bold || func != funcold) {
		init=1;
		aold=a;
		bold=b;
		funcold=func;
		delta=(b-a)/M;
		for (j=1;j<=M+1;j++)
			data[j]=(*func)(a+(j-1)*delta);
		for (j=M+2;j<=NDFT;j++)
			data[j]=0.0;
		for (j=1;j<=4;j++) {
			endpts[j]=data[j];
			endpts[j+4]=data[M-3+j];
		}
		realft(data,NDFT,1);
		data[2]=0.0;
	}
	en=w*delta*NDFT/TWOPI+1.0;
	nn=IMIN(IMAX((int)(en-0.5*MPOL+1.0),1),NDFT/2-MPOL+1);
	for (j=1;j<=MPOL;j++,nn++) {
		cpol[j]=data[2*nn-1];
		spol[j]=data[2*nn];
		xpol[j]=nn;
	}
	polint(xpol,cpol,MPOL,en,&cdft,&cerr);
	polint(xpol,spol,MPOL,en,&sdft,&serr);
	dftcor(w,delta,a,b,endpts,&corre,&corim,&corfac);
	cdft *= corfac;
	sdft *= corfac;
	cdft += corre;
	sdft += corim;
	c=delta*cos(w*a);
	s=delta*sin(w*a);
	*cosint=c*cdft-s*sdft;
	*sinint=s*cdft+c*sdft;
	free_vector(cpol,1,MPOL);
	free_vector(spol,1,MPOL);
	free_vector(xpol,1,MPOL);
}
#undef M
#undef NDFT
#undef MPOL
#undef TWOPI
#undef NRANSI