Multipole works by first interpolating (e.g., psi4) data onto a uniform grid in theta and phi, such that
theta = [0,pi], where theta = 0,dth,2 dth, …, N dth=pi
phi = [0, 2 pi], where phi = 0, dph, 2 dph, …, N dph=2 pi
Clearly these are not midpoints. Further the integration for midpoint method proceeds as (https://bitbucket.org/einsteintoolkit/einsteinanalysis/src/26d17d9edfecafd0c516e38fe8343bdc1d11eb5d/Multipole/src/integrate.cc#lines-51):
for (iy = 0; iy <= ny; iy++) for (ix = 0; ix <= nx; ix++) integrand_sum += f[idx(ix,iy)]; return hx * hy * integrand_sum;
which double counts points on the periodic boundaries. This might not be a problem for theta since the integrands naturally include a multiplication by sin(theta), and sin(0)=sin(pi)=0. However it does pose a major problem with phi, since there is no guarantee that the integrand at phi=0 or 2pi vanishes.
There is even a warning in code comments about this adding a first-order error to a method that’s “already first order”, while the trapezoidal rule is claimed to be “2nd order convergent” (param.ccl). When properly implemented the error order on the midpoint method matches that of the trapezoidal rule and in fact midpoint produces about half the error (see https://en.wikipedia.org/wiki/Riemann_sum#Midpoint_rule ; cf https://en.wikipedia.org/wiki/Riemann_sum#Trapezoidal_rule).
Despite the above problems, for some reason midpoint integration is chosen as the default.
I propose that we either error out if the user attempts to use the midpoint rule in its current form, or properly implement the midpoint method in the thorn.
Empirical effects of the bug: I find that certain odd-m modes extracted from QC-0 (e.g., l=2,m=1 or l=m=3) are of order 1% the amplitude the dominant (l=m=2) mode with the bug, but drop to 0.001% the amplitude when trapezoidal rule (a properly implemented integration method) is chosen.