Add a numerical integration method to use with line shapes and plasma models

[vsnever]

To obtain an accurate spectrum, the spectral emission must be averaged over the spectral bins. Currently, only Gaussian-based lineshape models accurately perform this averaging because they are implemented using the Erf functions. In all other cases, the spectral emission is calculated as a half-sum of the values ​​at the boundaries of the spectral bin. It would be nice to have a numerical integration method for integrating emission over spectral bins.

I think fixed-tolerance Gaussian quadrature will do. We should probably limit the maximum quadrature order to some reasonable value, because then we can use lookup tables for the quadrature coefficients.

[Mateasek]

Do you mean that this would be applied for every bin separately? If yes, shouldn't we also add some way for users to specify the maximum size of spectral step, the lineshape evaluation should done at before recalculating to the ray spectrum bins? An example could be a parameter specifying the step as a fraction of line width. For the cases, when ray spectrum has a smaller step, this could be skipped.

[vsnever]

Yes, this would be applied to every bin separately. The idea is to replace, for example, this: https://github.com/cherab/core/blob/0eb7454137b17bfbeca6d5e3eae4281f9cd4c157/cherab/core/model/lineshape.pyx#L416-L430 with an integration method, so the code would look like this:

lower_wavelength = raw_lineshape.min_wavelength + start * raw_lineshape.delta_wavelength
for i in range(start, end):
    upper_wavelength = raw_lineshape.min_wavelength + raw_lineshape.delta_wavelength * (i + 1)
    raw_lineshape.samples_mv[i] += integration_method(self._lineshape_func, lower_wavelength, upper_wavelength, self._rtol, self._maxiter) / (upper_wavelength - lower_wavelength)
    lower_wavelength = upper_wavelength

Regarding the size of spectral step, if we implement the Gaussian quadrature method, the spectral bin is divided according to the roots of Legendre polynomial of a given order, so the method just iterates over the polynomials until the desired relative tolerance is achieved or until the maximum order of the polynomial is reached. Or did I not understand the question?

Basically, I just want to implement a slightly simplified Cython version of this scipy function https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quadrature.html.

[Mateasek]

Thanks for explanation @vsnever, this looks great!

[vsnever]

Implemented in #378.