ENH: use our quadrature routines instead of fixed_quad

Scipy's fixed_quad routine is (as far as I have been able to tell) equivalent to our legendre routine with 5 points.

When I used our own routines in Julia instead of the quadrature that comes in Julia (quadgk) I got huge speedups. I'm hoping we can see a similar performance boost in python, but we would have to do some benchmarks. One place we could really get speedups is when we integrate over the same region inside a for loop, or nested for loops. In this case we could just compute the integration nodes and weights one time and to the quadrature by simply evaluating the integrand at the nodes and having numpy compute a dot product between that result and the weights. Not recomputing nodes/weights each iteration can only help.

Disclaimer: the quadgk routine in Julia is an adaptive routine, so my little performance comparison wasn't exactly fair. However, for all the tests I threw at the two routines they produced effectively equivalent results for the integral .

Locate the uses of SciPy's fixed_quad and the project's legendre routine first. Benchmark repeated and nested integrations, including reuse of nodes and weights, and compare both results and runtime; done means equivalent integrals with documented performance measurements.