Voigt profile functions wrong in Sherpa

[olaurino]

From a CXC help desk user:

I am pretty sure I have found two quite serious mistakes in the functions defined in Sherpa as "absorptionvoigt" and emissionvoigt". These are supposed to be (copying-and-pasting from the documentation):
"emission(/absorption) profile, where a Gaussian model of the core is combined with a Lorentzian model of the wings. The approximation presented in Astrophysical Formulae (K. R. Lang, 1980, 2nd ed., p. 220) is used. "

So, I went to check the reference, where the approximation is given and compared it to the lines of Ciao code at ~/ciao/lib/python/site-packages/sherpa/astro/optical/int.py.
Here is the relevant (for the emissionvoigt model) part of that file:

# This model computes emission as a Voigt function -- i.e., with
# a Gaussian core and Lorentzian wings.
class EmissionVoigt(ArithmeticModel):
"""Emission Voigt function."""

def __init__(self, name='emissionvoigt'):
self.center = Parameter(name, 'center', 5000., tinyval, hard_min=tinyval, frozen=True, units="angstroms")
self.flux = Parameter(name, 'flux', 1.)
self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s")
self.lg = Parameter(name, 'lg', 1., tinyval, hard_min=tinyval, frozen=True)

# Create core and wings from Gaussian and Lorentz
self._core = EmissionGaussian()
self._wings = EmissionLorentz()

ArithmeticModel.__init__(self, name, (self.center, self.flux, self.fwhm, self.lg))

#@modelCacher1d
def calc(self, p, x, xhi=None, **kwargs):
# Combining two emission components means
# adding them, it appears (at least according
# to addEmission in NarrowBandFunction.java)

core_pars = numpy.array([ p[2], p[0], p[1], 1.0])
wing_pars = numpy.array([ p[3]*p[2], p[0], p[1], 2.0])
return self._core.calc(core_pars, x) + self._wings.calc(wing_pars, x)

So, this function in Sherpa computes a voigt profile as the sum of a Gaussian and a Lorentizan: THIS IS WRONG!
A Voigt profile is the convolution product (very different from a sum, mathematically) of a Gaussian and a Lorentzian, i.e. the weighted integral of the Lorentzian over the Gaussian profile. And it is the same thing for emission and absorption (the only difference is that in one case the profile is added to the continuum and in the other it is subtracted from it).

The approximation in Lang (1980) is not the sum of the Gaussian and the Lorentzian, but the sum of the Gaussian and an 1/(x-x0)^2 like function function, weighted with the product of the FWHMs of the two original functions (the Gaussian and the Lorentzian).

The piece of code where the "absorptionvoigt" function is defined (see below) is even more seriously incorrect, I think: because there the Gaussian and the Lorentzian are multiplied by each other!
Again, there is no difference between absorption and emission: the profile is defined in exactly the same way: it only depends on whether the user "adds" it to the continuum (emission) or subtract it from it (absorption).

# This model computes absorption as a Voigt function -- i.e., with
# a Gaussian core and Lorentzian wings.
class AbsorptionVoigt(ArithmeticModel):
"""Absorption Voigt function expressed in equivalent width."""

def __init__(self, name='absorptionvoigt'):
self.center = Parameter(name, 'center', 5000., tinyval, hard_min=tinyval, frozen=True, units="angstroms")
self.ew = Parameter(name, 'ew', 1., tinyval, hard_min=tinyval, units="angstroms")
self.fwhm = Parameter(name, 'fwhm', 100., tinyval, hard_min=tinyval, units="km/s")
self.lg = Parameter(name, 'lg', 1., tinyval, hard_min=tinyval)

# Create core and wings from Gaussian and Lorentz
self._core = AbsorptionGaussian()
self._wings = AbsorptionLorentz()

ArithmeticModel.__init__(self, name, (self.center, self.ew, self.fwhm, self.lg))

#@modelCacher1d
def calc(self, p, x, xhi=None, **kwargs):
# Combining two absorption components means
# multiplying them, it appears (at least according
# to addAbsorption in NarrowBandFunction.java)

core_pars = numpy.array([ p[2], p[0], p[1] / 2.0, 4.0 ])
wing_pars = numpy.array([ p[2] * p[3], p[0], p[1] / 2.0 ])
return self._core.calc(core_pars, x) * self._wings.calc(wing_pars, x)

[olaurino]

@DougBurke @anetasie on a more general note it is not the first time we find issues with the models we got from specview a few years ago, maybe we should review these models?

[olaurino]

I added this to the 4.11.1 milestone as this seems pretty bad and I want it to be on my radar.

[anetasie]

@olaurino yes, the review of all the models is needed.

[anetasie]

Here is the link to the python implementation of voigt profile.
https://scipython.com/book/chapter-8-scipy/examples/the-voigt-profile/