/* Numerical Recipies Integration Routines */

#include 
#include 
#include "integrat.h"

#define FUNC(x) ((*func)(x))

#define EPS 1.0e-6      /* default fractional accuracy */
#define JMAX 20		    /* 2^(JMAX-1) is max. no. steps */

double trapzd(double (*func)(double),double a,double b,int n)
{			/* trapezoidal rule */
  double x,tnm,sum,del;
  static double s;
  static int it;
  int j;

  if (n==1) {
    it=1;                   /* no. of points to be added on next call */
    return (s=0.5*(b-a)*(FUNC(a)+FUNC(b)));
  }else{
    tnm=it;
    del=(b-a)/tnm;          /* spacing of points to be added */
    x=a+0.5*del;
    for (sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x);
    it *= 2;
    s=0.5*(s+(b-a)*sum/tnm);  /* replaces s by refined value */
    return s;
  }
}

double qsimp(double (*func)(double),double a,double b,double eps)       /* Simpson's rule */
{
  int j;
  double s,st,ost,os,trapzd(double (*func)(double),double a,double b,int n);

  if((eps<=0.0)||(eps>=1.0)) eps=EPS;
  ost=os=-1.0e30;
  for (j=1;j<=JMAX;j++) {
    st=trapzd(func,a,b,j);
    s=(4.0*st-ost)/3.0;
	if (fabs(s-os)<=eps*fabs(os)) return s;
    os=s;
    ost=st;
  }
  fprintf(stderr,"Too many steps in routine QSIMP");
  return(0);
}