Add 2PCF multipoles

[lgarrison]

Is your feature request related to a problem? Please describe.
Corrfunc has no feature to compute 2PCF multipoles. Users can use halotools' wrapper around DDsmu, but that requires tuning the number of mu bins. Too few and the discrete approximation to the integral is inaccurate; too many, and the estimator gets noisy and the pair counting slows down.

Describe the solution you'd like
Corrfunc could implement multipoles directly by weighting pairs by the Legendre polynomials in mu. No binning in mu required.

In theory, we could compute multiple multipoles simultaneously, but that would require support for returning multiple weighted results. This is similar to the issue of jackknife covariance, where we'd like to return one result for each round of leave-one-out. But for multipoles, we'd only have to return a few results, rather than hundreds.

It's unclear what the best way to implement this is. Either as a new module, "xi_ell", or as a custom weighting function for DDsmu. Either way, we'd like to support returning multiple weighted values.

Describe alternatives you've considered
Users can use the halotools wrapper for most applications. But for sparse data sets where the binning would be noisy, or large data sets where a large number of bins is costly, this is less than ideal. And all use of the halotools wrapper requires validation of the number of mu bins.

Additional context
plot_wp_xi0_xi2_xi4_reseeded_varmubins.pdf

Plot from @SandyYuan showing the dependence of multipole on number of mu bins. The issue is most noticeable in the hexadecapole, where one would probably need several hundred mu bins to be convinced of convergence.