/* DSAQ_:
 * A double-precision adaptive quadrature routine, using Simpson's rule as the
 * elemental integration technique.
 *
 * dsaq_ is written to be called from a FORTRAN routine.  To this end,
 * interface arguments are passed by reference and an underscore is appended to
 * the end of the function name in the C code.  This underscore should not
 * appear in the FORTRAN calling routine.  (So, the calling routine does not
 * even know that it is calling a C routine.)
 *
 * This routine is written in C since it uses recursion, which F77 does not support.
 *
 * John Kerl
 * 10/12/93
 */

#define FABS(x) ((x)>0.0?(x):(-x))

/**********************************************************************/
double dsaq_(double (*f) (double *x), double *a, double *b, double *epsilon)
{
	double d_simpson (double (*)(double *), double, double);
	double c, s1, s2, half_eps, diff;
	static int first_call = 1;
	/* static int counts = 0; */

	if (first_call) {
		*epsilon *= 10.;
		first_call = 0;
	}
	/* printf("\ndsaq %5d: a=%10.7e, b=%10.7e, e=%10.7e\n", counts++, *a, *b,
	 *epsilon);
	 */
	c = (*a + *b) / 2.;
	half_eps = *epsilon / 2.;

	s1 = d_simpson (f, *a, *b);
	s2 = d_simpson (f, *a, c) + d_simpson (f, c, *b);
	/*
		 printf("s1=%10.7e, s2=%10.7e, diff=%10.7e,\n\t e=%10.7e, h=%10.7e\n",
		 s1, s2, FABS(s1-s2), *epsilon, half_eps);
	 */
	/* if (counts > 4) { */
	/* exit (1); */
	/* } */
	diff = (FABS (s2 - s1));
	if (diff <= *epsilon) {
		return s2;
	}
	else {
		return dsaq_ (f, a, &c, &half_eps) + dsaq_ (f, &c, b, &half_eps);
	}
}

double d_simpson (double (*f) (double *x), double a, double b)
{
	double c, h, *ap, *bp, *cp;

	h = (b - a) / 2.;
	c = a + h;
	ap = &a;
	bp = &b;
	cp = &c;

	return (h / 3. * ((*f) (ap) + 4. * (*f) (cp) + (*f) (bp)));
}

/**********************************************************************/
float ssaq_ (float (*f) (float *x), float *a, float *b, float *epsilon)
{
	float s_simpson (float (*)(float *), float, float);
	float c, s1, s2, half_eps, diff;
	static int first_call = 1;
	/* static int counts = 0; */

	if (first_call) {
		*epsilon *= 10.;
		first_call = 0;
	}
	/* printf("\nssaq %5d: a=%10.7e, b=%10.7e, e=%10.7e\n", counts++, *a, *b,
	 *epsilon);
	 */
	c = (*a + *b) / 2.;
	half_eps = *epsilon / 2.;

	s1 = s_simpson (f, *a, *b);
	s2 = s_simpson (f, *a, c) + s_simpson (f, c, *b);
	/*
		 printf("s1=%10.7e, s2=%10.7e, diff=%10.7e,\n\t e=%10.7e, h=%10.7e\n",
		 s1, s2, FABS(s1-s2), *epsilon, half_eps);
	 */
	/* if (counts > 4) { */
	/* exit (1); */
	/* } */
	diff = (FABS (s2 - s1));
	if (diff <= *epsilon) {
		return s2;
	}
	else {
		return ssaq_ (f, a, &c, &half_eps) + ssaq_ (f, &c, b, &half_eps);
	}
}

float s_simpson (float (*f) (float *x), float a, float b)
{
	float c, h, *ap, *bp, *cp;

	h = (b - a) / 2.;
	c = a + h;
	ap = &a;
	bp = &b;
	cp = &c;

	return (h / 3. * ((*f) (ap) + 4. * (*f) (cp) + (*f) (bp)));
}