Driscoll&Healy integration for Multipole

[Erik Schnetter]

The attached patch implements a more accurate integration over the sphere, using an algorithm by Driscoll & Healy. This algorithm uses Gaussian integration weights, leading (almost) to exponential convergence.

Keyword: Multipole

[Roland Haas]

this seems to be buggy. The results differ by what looks like Pi from the results of test_simpson.par which is otherwise an identical parfile. Testing this with a constant value fr[i]=1.0*sin(th) in utils.cc also yields an answer 3.14 times too large when driscoll&Healy is used as compared to Simpson (which gives the correct 4*Pi result).

I don't know how to fix this since I have no idea where the extra Pi creeps in.

Do not apply in this form. :-)

[Erik Schnetter]

I found the missing factor pi. This factor is actually present in the (unnumbered!) equations in the Driscoll&Healy paper. In my tests, I left it out because I implicitly assumed that theta ranges [0...1], whereas of course it ranger [0...pi]. This made the results too large by a factor of pi.

There was also an error in the index bounds. Driscoll&Healy want a point on the equator, which corresponds to even ntheta, not odd ntheta as previously in the code.

Applied.