double integrate_pp(double (*f)(void *,double),
    void *params, double a, double b, int n)
{
  int i;
  double sum;
  sum = 0.5*((*f)(params,a) + (*f)(params,b));
  for ( i = 1; i <= n-1; i++ )
    sum = sum + (*f)(params,a+i*(b-a)/n);
  return sum*(b-a)/n;
}