/* Numerical Recipies Integration Routines */ double trapzd(double (*func)(double),double a,double b,int n); /* trapezoidal */ double qsimp(double (*func)(double),double a,double b,double eps); /* Simpson's rule */