Lbl4IR – Line-by-line for InfraRed

Optimized implementations of Rational Approximations for the Voigt and Complex Error Function

(Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 112, pp. 1010-1025, 2011)

Abstract

Rational functions are frequently used as efficient yet accurate numerical approximations for real and complex valued functions. For the complex error function w(x+iy), whose real part is the Voigt function K(x,y), code optimizations of rational approximations are investigated. An assessment of requirements for atmospheric radiative transfer modeling indicates a y range over many orders of magnitude and accuracy better than 10^-4. Following a brief survey of complex error function algorithms in general and rational function approximations in particular the problems associated with subdivisions of the x, y plane (i.e., conditional branches in the code) are discussed and practical aspects of Fortran and Python implementations are considered. Benchmark tests of a variety of algorithms demonstrate that programming language, compiler choice, and implementation details influence computational speed and there is no unique ranking of algorithms. A new implementation, based on subdivision of the upper half-plane in only two regions, combining Weideman's rational approximation for small |x|+y < 15 and Humlicek's rational approximation otherwise is shown to be efficient and accurate for all x, y.

(See the weideman.py module for Python implementations of various Humlicek-Weideman combinations.)

Computational aspects of speed-dependent Voigt profiles

(Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 187, pp. 44-53, 2017)

Abstract

The increasing quality of atmospheric spectroscopy observations has indicated the limitations of the Voigt profile routinely used for line-by-line modeling, and physical processes beyond pressure and Doppler broadening have to be considered. The speed-dependent Voigt (SDV) profile can be readily computed as the difference of the real part of two complex error functions (i.e. Voigt functions). Using a highly accurate code as a reference, various implementations of the SDV function based on Humlicek's rational approximations are examined for typical speed dependences of pressure broadening and the range of wavenumber distances and Lorentz to Doppler width ratios encountered in infrared applications. Neither of these implementations appears to be optimal, and a new algorithm based on a combination of the Humlicek (1982) and Weideman (1994) rational approximations is suggested.

Supplement Material

demo.py (called mmc1.txt in the JQSRT supplement): A (I)Python interactive session demonstrating the use of the humlicek, weideman, and sdv modules and how to generate the contour plots.
sdv.py (called mmc4.txt in the JQSRT supplement): The 'main functions' of this module, sdvProfile and voigtProfile, evaluate the corresponding profiles (normalized to one) for given pressure & temperature, and line position & widths.
Further functions defined here provide implementations of the SDV function based on various codes (and combinations therof) for the complex error function.
humlicek.py (called mmc2.txt in the JQSRT supplement): Various implementations of Humlicek's cpf12 and w4 rational approximations
weideman.py (called mmc5.txt in the JQSRT supplement): The rational approximations with N=24 and N=32 terms and the Humlicek-Weideman combinations.
ir.py (called mmc3.txt in the JQSRT supplement): Mathematical and physical constants

Comments on the Voigt function implementation in the Astropy and SpectraPlot.com packages

Journal of Quantitative Spectroscopy and Radiative Transfer, 213, 13-16, 2018;
Preprint available at arXiv:1806.10338

Abstract

The Voigt profile is important for spectroscopy, astrophysics, and many other fields of physics, but is notoriously difficult to compute. McLean et al. [J. Electron Spectrosc. & Relat. Phenom., 1994] have proposed an approximation using a sum of Lorentzians. Our assessment indicates that this algorithm has significant errors for small arguments. After a brief survey of the requirements for spectroscopy we give a short list of both efficient and accurate codes and recommend implementations based on rational approximations.

Supplement Material and Further Links

mcLean.py: The Python script defining the Voigt function approximations and generating the figures.
astropy: The Astropy Collaboration, T.P. Robitaille, E.J. Tollerud, et al.
Astropy: A community Python package for astronomy
Astronomy & Astrophysics, Vol. 558, Article A33, October 2013
spectraPlot.com: C.S. Goldenstein, V.A. Miller, R.M. Spearrin, C.L. Strand:
SpectraPlot.com: Integrated spectroscopic modeling of atomic and molecular gases
J. Quant. Spectrosc. & Radiat. Transfer, 200, pp. 249-257, 2017

The Voigt and complex error function: Humlicek's rational approximation generalized

(Monthly Notices of the Royal Astronomical Society, 479(3), 3068-3075, 2018;
Preprint available at arXiv:1806.11560

Abstract

Accurate yet efficient computation of the Voigt and complex error function is a challenge since decades in astrophysics and other areas of physics. Rational approximations have attracted considerable attention and are used in many codes, often in combination with other techniques. The 12-term code "cpf12" of Humlicek (1979) achieves an accuracy of five to six significant digits throughout the entire complex plane. Here we generalize this algorithm to a larger (even) number of terms. The n = 16 approximation has a relative accuracy better than 10^-5 for almost the entire complex plane except for very small imaginary values of the argument even without the correction term required for the cpf12 algorithm. With 20 terms the accuracy is better than 10^-6. In addition to the accuracy assessment we discuss methods for optimization and propose a combination of the 16-term approximation with the asymptotic approximation of Humlicek (1982) for high efficiency.

Supplement Material

cpfX.py --- Various implementations of the Humlicek (1979) algorithm:

Python/NumPy implementations of the original Fortran 77 source code (with two versions to combine the two regions);
A generalization of the region I rational approximation to arbitrary (even) number of terms;
An optimized implementation of the region I approximation for 16 terms using a single fraction;
The combination Eq. (12) with the asymptotic approximation of Humlicek (1982) for NumPy.
NOTE: line 127 does not work in Python 3 because integers are required for array indices! This line correctly reads as
for t,a,b in zip(roots[n//2:],alfa[n//2:],beta[n//2:]):
UPDATE Nov 2022: new implementation of the zpf16 algorithm, including results of recent timeit benchmarks

fig124.py

Execution of this script inside the IPython interpreter should produce the Figures 1, 2, and 4 of the manuscript.

weideman.py (originally provided as supplementary material of the JQSRT 2017 paper)

Weideman's rational approximations with N=24 and N=32 terms and the Humlicek-Weideman combinations.

cpf12.f90

Shigenobu Hirose (Japan Agency for Marine-Earth Science and Technology) has identified several errors (typos) in the Fortran code providing a modified version of Humlicek's cpf.f subroutine.

lines 76-77 (== 109-110) are incorrect, the index of alfa, beta, and t should be 6 instead of 5.

in lines 118-127 the t(1) in the second nominator term is incorrect and should be t(2), t(3), ... etc.

Computational aspects of speed-dependent (Voigt and) Rautian profiles

(co-authord with Philipp Hochstaffl;
Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 258, 107385, 2021)
Preprint available at arXiv:2010.09804

Abstract

For accurate line-by-line modeling of molecular cross sections several physical processes "beyond Voigt" have to be considered. For the speed-dependent Voigt and Rautian profiles (SDV, SDR) and the Hartmann-Tran profile the difference w(iz₋)-w(iz₊) of two complex error functions (essentially Voigt functions) has to be evaluated where the function arguments z± are given by the sum and difference of two square roots. These two terms describing z± can be huge and the default implementation of the difference can lead to large cancellation errors. First we demonstrate that these problems can be avoided by a simple reformulation of z₋. Furthermore we show that a single rational approximation of the complex error function valid in the whole complex plane (e.g. by Humlicek, 1979 or Weideman, 1994) enables computation of the SDV and SDR with four significant digits or better. Our benchmarks indicate that the SDV and SDR functions are about a factor 2.2 slower compared to the Voigt function, but for evaluation of molecular cross sections this time lag does not significantly prolong the overall program execution because speed-dependent parameters are available only for a fraction of strong lines.

Supplement Material

sdv.py (called mmc1.txt in the JQSRT supplement): Various implementations of the speed-dependent Voigt function, essentially the versions discussed in subsection 3.1 and 3.2 investigating different ways to compute the two complex square roots and their difference.
NOTE: The arguments defined here are slightly different from those defined in the old 2017 SDV paper!!! (Here the third argument is essentially the ratio of the second order broadening parameter and the Gauss width γ₂/γ_G, whereas in the old paper we used α = γ_L/γ₂ -3/2.)
sdr.py (called mmc2.txt in the JQSRT supplement)

Py4CAtS – Python for Computational ATmospheric Spectroscopy

Abstract

Py4CAtS is a Python re-implementation of the modern Fortran General Purpose Atmospheric Radiative Transfer Line-by-Line Infrared-Microwave Code (GARLIC), where compute-intensive code sections utilize the Numeric/Scientific Python modules (Numpy / Scipy) for highly optimized array-processing. The individual steps of an infrared or microwave radiative transfer computation are implemented in separate scripts to extract lines of relevant molecules in the spectral range of interest, to compute line-by-line cross sections for given pressure(s) and temperature(s), to combine cross sections to absorption coefficients and optical depths, and to integrate along the line-of-sight to transmission and radiance/intensity. Visualization of all these quantities is done with matplotlib.

Some notes on Voigt and complex error function (cerf) implementations

The "old" scripts (e.g. for the SDV 2017 paper) were written for Python 2, whereas the newer scripts use Python 3
There is some inconsistency w.r.t. to the arguments of the various "cerf" (w) functions (here the second "function" refers to the Python code): For some codes (e.g. in the context of the SDV and SDR papers) we use one complex argument z=x+iy; whereas in other cases two real arguments x, y are used. For lbl number crunching of molecular cross sections with the Voigt profile the two argument version is preferable, because each individual spectral line is computed for an array of wavenumber grid points (corresponding to an x grid array) and a single scalar Voigt parameter y characterizing the broadening (calling the cerf function with a complex argument would almost double the size of the array to be passed).
For (some of) the "either-or" algorithms combining the Humlicek (1982) asymptotic approximation (the region I (or II) of the w4 code) with the rational approximation of Weideman (1994) or Humlicek (1979) (the generalized cpf code) there are two versions:
- a version using numpy's mask function (for time-critical applications, e.g. hum1zpf16m)
- a version using numpy's where function (for contour plots, e.g. hum1zpf16)
Code for the Humlicek (1979 and 1982) and Weideman (1994) rational approximations can be imported from the files given above. Furthermore, Py4CAtS has a module racef.py of "Rational Approximations for the Complex Error Function". For convenience we also provide a module test_racef.py here that can be used to generate contour plots as in the papers above. (See test_racef.html for a typical output.)

Documentation

The Voigt and complex error function: Humliceks rational approximation generalized. MNRAS, 479(3), 3068-3075, 2018 (arxiv: physics.comp-ph/1806.11560)
Computational Aspects of Speed-dependent Voigt Profiles. JQSRT, Vol. 187, pp. 44-53, 2017
Optimized Implementations of Rational Approximations for the Voigt and Complex Error Function. JQSRT, Vol. 112(6), pp. 110-125, 2011
D. Kohlert and FS: Line--by--Line Computation of Atmospheric Infrared Spectra with Field Programmable Gate Arrays. IEEE-JSTARS, 4(3), 701-709, 2011
Comments on "A Common Misunderstanding about the Voigt Line Profile." J. Atmos. Science 66(12), 1860-1864, 2009 (Erratum)
FS and D. Kohlert: Optimized Implementations of Rational Approximations --- A Case Study on the Voigt and Complex Error Function. CPC 179(7), 457-465, 2008
The Voigt and complex error function: A comparison of computational methods. JQSRT, Vol. 48, pp. 743-762, 1992
GARLIC --- A General Purpose Atmospheric Radiative Transfer Line-by-Line Infrared-Microwave Code: Implementation and Evaluation; JQSRT, Vol. 127, pp. 29-50, 2014
Algorithmic vs. finite difference Jacobians for infrared atmospheric radiative transfer; JQSRT, Vol. 164, pp. 147-160, 2016
Py4CAtS --- Python for Computational ATmospheric Spectroscopy; Atmosphere, Vol. 10(5), pp. 262, 2019
An assessment of some closed-form expressions for the Voigt function: III: Combinations of the Lorentz and Gauss functions. JQSRT, 226, 87-91, 2019 (arxiv: physics.comp-ph/1901.08366)
Comments on the Voigt function implementation in the Astropy and SpectraPlot.com packages. JQSRT, 213, 13-16, 2018 (arxiv: physics.comp-ph/1806.10338)
An assessment of some closed-form expressions for the Voigt function: II: Utilizing rational approximations for the Gauss function. JQSRT, 202, 81-89, 2017
An assessment of some closed-form expressions for the Voigt function JQSRT, 176, 1-5, 2016
J.Humlicek: An efficient method for evaluation of the complex probability function: the Voigt function and its derivatives.; JQSRT, Vol. 21, 309-313, 1979
J.Humlicek: Optimized computation of the Voigt and complex probability function.; JQSRT, Vol. 27, 437-444, 1982
J.A.C. Weideman: Computation of the Complex Error Function.; SIAM-J. Num. Anal., Vol. 31, 1497-1518, 1994

For more applications of Lbl-IR see also the listing at

What's new?

October 2016: Initial page setup
April 2018: JQSRT 2018 paper added
June 2018: MNRAS 2018 paper added
August 2019: "Assessment of some closed-form expressions ..." paper series added
December 2020: Speed-dependent Rautian paper (JQSRT) added
February 2021: Broken link corrected, some more papers
August 2021: Corrections of the supplementary files cpfX.py and cpf12.f90 of the MNRAS 2018 paper on Humlicek's rational approximation generalized
November 2022: Supplementary file cpfX.py of the "Generalized Humlicek" paper (MNRAS 2018) updated