[SciPy-Dev] Proposal for a new scipy.special function: Voigt

[Warren Weckesser]

On 9/23/18, Muhammad Firmansyah Kasim <firman.kasim at gmail.com> wrote:
> I am proposing to implement a new special function under `scipy.special`
> module: Voigt function.
>
> Voigt function/profile is a very common line profile in atomic transitions
> and widely used in the spectroscopy field. It is a convolution between
> Gaussian and Cauchy distribution, thus it is not trivial to calculate it
> efficiently.
>
> There are several approximations to calculate it efficiently, but I am
> proposing to use the one in the "Stellar atmosphere" by Mihalas (1978). It
> is the one used in FLYCHK by NIST and IAEA (I work closely with the author
> of the code).
>
> Please let me know if you have any suggestions or objections to the
> proposal.
> Thanks!
>
> Muhammad
>


Thanks, Muhammad.  I've seen questions about implementing the Voigt
function more than once over on stackoverflow, so I think it would be
a worthwhile addition to scipy.

The Voigt function can be expressed in terms of the complex Faddeeva
function w(z), which is in scipy.special.  Using w(z) allows the Voigt
function to be implemented in just a few lines.

I experimented with the Voigt function when working on an answer to a
question at stackoverflow:
https://stackoverflow.com/questions/40536467/scipy-integrate-quad-not-returning-the-expected-values

Below is a modified version of the script in my stackoverflow answer.
It generates the attached plot, which reproduces the plot at
https://dlmf.nist.gov/7.19.  It would be good to know how the accuracy
and performance of these implementations compare to the code that you
are proposing.

Warren

-----
from __future__ import division

import numpy as np
from scipy.special import wofz, erfc
from scipy.integrate import quad


_SQRTPI = np.sqrt(np.pi)
_SQRTPI2 = _SQRTPI/2


def voigtuv(x, t):
    """
    Voigt functions U(x,t) and V(x,t).

    The return value is U(x,t) + 1j*V(x,t).
    """
    sqrtt = np.sqrt(t)
    z = (1j + x)/(2*sqrtt)
    w = wofz(z) * _SQRTPI2 / sqrtt
    return w


def voigth(a, u):
    """
    Voigt function H(a, u).
    """
    x = u/a
    t = 1/(4*a**2)
    voigtU = voigtuv(x, t).real
    h = voigtU/(a*_SQRTPI)
    return h


if __name__ == "__main__":
    # Create plots of U(x,t) and V(x,t) like those
    # in http://dlmf.nist.gov/7.19

    import matplotlib.pyplot as plt

    xmax = 15
    x = np.linspace(0, xmax, 100)
    t = np.array([0.1, 2.5, 5, 10])
    y = voigtuv(x, t.reshape(-1, 1))

    colors = ['g', 'r', 'b', 'y']

    plt.subplot(2, 1, 1)
    for tval, ya, clr in zip(t, y, colors):
        plt.plot(x, ya.real, clr, label="t=%4.1f" % tval)
    plt.grid()
    plt.xlim(0, xmax)
    plt.ylim(0, 1)
    plt.legend(loc='best')
    plt.title('U(x,t)')

    plt.subplot(2, 1, 2)
    for tval, ya, clr in zip(t, y, colors):
        plt.plot(x, ya.imag, clr, label="t=%4.1f" % tval)
    plt.grid()
    plt.xlim(0, xmax)
    plt.ylim(0, 0.5)
    plt.legend(loc='best')
    plt.title('V(x,t)')
    plt.xlabel('x')

    plt.tight_layout()
    plt.show()

-----
-------------- next part --------------
A non-text attachment was scrubbed...
Name: voigtuv.png
Type: image/png
Size: 114375 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20180923/b6e41081/attachment-0001.png>