GetMultipoleIntegral(k,f1,f2,N) calculates the k-th multipole moment $\int f_1(r)r^kf_2(r)\mathrm{d}r$
The Integral is calculated over the full domain of the interpolation function f1, where a Gaussian quadrature rule with N points between each two knots of the interpolated function is used.