#include <math.h>
#define EPS 1.0e-3
#define JMAX 14
#define JMAXP (JMAX+1)
#define K 5

float midpnt(float (*func)(float theta, int mode, float *ca, float *fa, int
    *ifa, int nf, float c), float a, float b, int n, int mode, float *ca,
    float *fa, int *ifa, int nf, float c);

float fqromo(float (*func)(float theta, int mode, float *ca, float *fa, int
	*ifa, int nf, float c), float a, float b, int mode, float *ca, 
        float *fa, int *ifa, int nf, float c)
{
	void polint(float xa[], float ya[], int n, float x, float *y, float *dy);
	void nrerror(char error_text[]);
	int j;
	float ss,dss,h[JMAXP+1],s[JMAXP];

	h[1]=1.0;
	for (j=1;j<=JMAX;j++) {
		s[j]=midpnt(func,a,b,j, mode, ca, fa, ifa, nf, c);
		if (j >= K) {
			polint(&h[j-K],&s[j-K],K,0.0,&ss,&dss);
			if (fabs(dss) <= EPS*fabs(ss)) return ss;
		}
		h[j+1]=h[j]/9.0;
	}
/*	nrerror("Too many steps in routing qromo"); */
	return 0.0;
}
#undef EPS
#undef JMAX
#undef JMAXP
#undef K