A two-domain MATLAB implementation for efficient computation of the Voigt/complex error function

In this work we develop a new algorithm for the efficient computation of the Voigt/complex error function. In particular, in this approach we propose a two-domain scheme where the number of the interpolation grid-points is dependent on the input parameter $y$. The error analysis we performed shows that the MATLAB implementation meets the requirements for radiative transfer applications involving the HITRAN molecular spectroscopic database. The run-time test shows that this MATLAB implementation provides rapid computation, especially at smaller ranges of the parameter $x$.


Introduction
The complex error function, also known as the Faddeeva function, is defined as [1][2][3][4] w (z) = e −z 2 where z = x + iy is a complex argument. The complex error function w (z) is closely related to the complex probability function [2] The complex probability function can be written in terms of its real and imaginary parts [2] W (z) = K (x, y) + iL (x, y) such that and Both functions w (z) and W (z) are equal to each other on the upper half of the complex plane, when y = Im [z] > 0 [2]. Consequently, it follows that w (z) = K (x, y) + iL (x, y) , y > 0.
Further we will imply that the parameter y = Im [z] is always greater than zero. The real part of the complex error function K (x, y) is known as the Voigt function [1][2][3][4] that is widely used in Atmospheric Science to describe emission and absorption of the photons by atmospheric molecules [5][6][7][8][9][10][11][12]. Specifically, the Voigt function is used to compute wavelength-dependent absorption coefficients by using the HITRAN molecular spectroscopic database [13]. The imaginary part of the complex error function L (x, y) is also used in many applications. For example, it can describe the spectral behavior of the index of refraction in various materials [14,15].
In our previous publication we proposed an algorithm for the efficient computation of the complex error function based on a single-domain implementation with vectorized interpolation [33]. However, despite rapid performance it has some limitations, including a limited range for the parameter x as well as restrictions for the input array type. As a further development, in this work we present a new MATLAB implementation without those drawbacks. The numerical analysis and computational tests we performed show that the proposed algorithmic implementation meets all the requirements in terms of accuracy and run-time performance for the efficient computation of the Voigt/complex error function in the radiative transfer applications.

Sampling Based Approximation
In our work [29] we have proposed a new sampling method based on incomplete cosine expansion of the sinc function. In particular, it is shown that, using a new sampling method based on incomplete cosine expansion of the sinc function, we can obtain the following approximation where the expansion coefficients are given by with parameters N , M , h and ς that can be taken as 23, 23, 0.25 and 2.75, respectively.
For more rapid performance in algorithmic implementation, it is reasonable to define the function (4) can be represented as [30] w (z) ≈ Ω (z + iς/2) .
Overall, approximation (5) is highly accurate. However, its accuracy deteriorates with decreasing parameter y. In order to resolve this problem we can use the following approximation [30] w where the expansion coefficients are It is interesting to note that this approximation of the complex error function is obtained by substituting approximation (4) into the right side of the identity (see [30] for details in derivation) For |z| > 8 one of the best choices is the approximation based on the Laplace continued fraction [1,3]. In particular, in our algorithm we used This algorithm is implemented in MATLAB as a script file fadsamp.m that utilizes three approximations (5), (6) and (7) as follows [30] As we have shown in our publication [30], approximation (8) provides highly accurate and rapid computation of the complex error function without poles and can be used to cover the entire complex plane.

Modified Trapezoidal Rule
In 1945, English mathematician and cryptanalyst Alan Turing, who succeeded to decrypt sophisticated machine codes of the Enigma during the Second World War [38], published an interesting paper where he proposed an elegant method of numerical integration for some class of integrals [39]. Nowadays, his method of the numerical integration that involves the residue calculus is regarded as the modified trapezoidal rule [40] or the generalized trapezoidal rule [37]. The comprehensive and detailed description of Turing's method of integration may be found in literature [40].
In 1949, Goodwin showed how to implement Turing's idea to the integrals of kind [41] ∞ Using the method described by Goodwin, Chiarella and Reichel in their work [42] derived the series expansion for the following integral where Ω (x, t) = (1 − ix) / 2t 1/2 . In particular, they showed that that the function Ψ (x, t) can be approximated as a series where h is a small fitting parameter and H (t) is the Heaviside step function defined as Thus, due to Heaviside step function, Equation (9) can be separated into three parts Equation (11) deals only with a single point h 2 /π 2 for the parameter t and does not represent any practical interest. Therefore, further we will consider only two equations, Equations (10) and (12).
Mata and Reichel [43] showed the relations that link both functions Ψ (x, t) and w (x, y) with each other. Consequently, using these relations the series expansions (10) and (12) can be reformulated as and respectively, where t k = (k + 1/2) h and h can be chosen to be equal to π/ (N + 1) [37].
Consider an expansion series for the complementary error function that was reported by Hunter and Regan [44] (see also [37]) Using the identity relating complex error function and complementary error function [3] w where τ k = kh [37]. Equations (13) and (14) have poles at z = t k while Equation (15) has poles at z = τ k . Furthermore, each of these equations can cover with high accuracy only in the corresponding domain. However, Al Azah and Chandler-Wilde [37] recently proposed the following approximation where ϕ (t) = t − t ∈ [0, 1), that appeared to be very efficient since it can be used for rapid and highly accurate computation without poles at N = 11. This is possible to achieve since, according to approximation (16), Equations (13), (14) and (15) are used interchangeably depending on the domain over the entire complex plain.

Algorithmic Implementation
Previously we have reported a new algorithm based on a vectorized interpolation over a single-domain [33]. Such an approach provides accuracy better than 10 −6 at y ≥ 10 −8 for the HITRAN [13] applications. However, this implementation has several limitations. In particular, there is a limitation |x| ≤ 10 5 . Although it is possible to increase the range for more than 10 5 , it requires to introduce more interpolation grid-points for precomputation. Furthermore, this MATLAB implementation accepts an input only as a vector (one dimensional array) or as a scalar.
One of the ways to implement an efficient algorithm is to use an approximation based on the two-domain scheme that we proposed in our earlier publication [23] K where the coefficients are [20] Consequently, the complex error function can also be approximated as [33] w Although the algorithms for the Voigt and complex error functions, built on Equations (17a) and (17b) can provide rapid computations, their accuracies deteriorate at y < 10 −6 .
In order to resolve this problem we developed a new algorithm that utilizes the internal MATLAB built-up features. Unlike interpolation algorithms shown in [23,33], the proposed approach implies that the number of the interpolation grid-points required for precomputation is not constant and dependent on the input parameter y such that where δ is the offset that was found experimentally to be 5 × 10 3 . The corresponding interpolation grid-points range is bounded by [−r, r], where r = 35 is the radius that was also found experimentally. The interpolation gridpoints are distributed non-equidistantly in logarithmic scale to increase their density near the origin along the x axis. Equation (18) is entirely empirical. It is reasonable to suggest that the number of the interpolation grid-points should be increased with decreasing y. However, taking the reciprocal of y n , such that n ≥ 1, results in a very rapid increase of the interpolation grid-points N gp with decreasing y. Our experimental results show that a more or less optimal value N gp can be found by taking the reciprocal of √ y. The parameter δ in Equation (18) is to provide a sufficient number of the interpolation grid-points N gp within the internal domain |x + iy| ≤ r as parameter y increases.
The use of Equation (18) is given as a basic for applications. If higher accuracy is required, the number of the interpolation grid-points N gp can be further increased by using, for example, a modified form of the equation above N gp = 2 √ y + 3δ.
Such a modification improves accuracy by an order of the magnitude. The algorithm utilizes two domains, internal and external, that are bounded by a circle of radius |x + iy| = r. The internal domain is situated within the circle while the external domain is situated outside of it.
All points within internal domain |x + iy| ≤ r are computed by the MAT-LAB built-in interpolation function interp1 though interpolation grid-points that are computed by using the function file fadsamp.m provided in our article [30]. The spline method was found to be the best for interpolation.
All points outside external domain |x + iy| > r are computed by the following approximation that represents a simplified version of Equation (7) above. The parameter y is dependent on temperature and pressure. Therefore, it is related to the atmospheric height. In a radiative transfer model, the atmospheric column is sliced into layers [5,6]. The layers can be considered homogeneous if their thicknesses are small enough. Typically, a thickness of 50 m (or 20 layers per km) is sufficient to consider the layer homogeneous for atmospheric modeling. This implies that, within a given layer, the temperature and pressure are constants. Since the density of air in the mesosphere is extremely low, its contribution to the radiance is practically negligible. Consequently, taking into consideration that only the troposphere and stratosphere, with an atmospheric column corresponding to a height of up to 40-50 km, actually contribute to the radiance, we can conclude that the total number of layers as well as elements for parameter y does not exceed 1000 per gas in atmosphere. However, in many practical tasks, for a single value y at a given atmospheric temperature and pressure, we may need to cover a wide spectral range with high resolution. This signifies that the number of elements for parameter x can exceed a million. Since the number of elements for parameter y is much smaller than the number of elements for parameter x, it is reasonable to implement a configuration in which parameter x is an array and parameter y is a scalar. This approach is more rapid in the MATLAB implementation since the application of a 2D array M × N , where N and M are numbers of elements for parameters x and y, respectively, requires considerably more computer memory.
As it has been mentioned above, the computation of the interpolation grid-points in its original version of the function file w2dom.m is performed by external function fadsamp.m. However, any other MATLAB function file that can provide highly accurate computation of the complex error function may also be used for computation of interpolation grid-points. For example, the function files such as fadf.m [27] and fadfunc.m [31] can also be used as alternatives. The script of the function file w2dom.m is given in Appendix A.

Error Analysis
In order to exclude the rounding and truncation errors in computation by using the most recent HITRAN database [13], the values of the K (x, y) and L (x, y) functions with 6 or more accurate decimal digits in their mantissas are required. Therefore, in radiative transfer applications involving the HITRAN database, the accuracy of computation of these functions has to be better than 10 −6 .  We can see that for the imaginary part the absolute difference also does not exceed 2.5 × 10 −13 . For more rigorous error analysis, we can apply the following relative errors for the real and imaginary parts of the complex error function w (z), respectively. Figure 3 depicts the relative error for the real part of the complex error function Re [z] = K (x, y) in the domain 0 ≤ x ≤ 50 and 0 ≤ y ≤ 10 1.699 , where 10 1.699 ≈ 50. Inset in Figure 3 shows the smaller domain near the origin 0 ≤ x ≤ 5.5 and 0 ≤ y ≤ 10 0.74 , where 10 0.74 ≈ 5.5. Beyond this smaller domain the observed color is blue indicating that the relative error is below 10 −12 . As we can see from Figure 3, the relative error does not exceed ∼10 −10 .    Figure 4 shows the smaller domain near origin 0 ≤ x ≤ 5.5 and 0 ≤ y ≤ 10 0.74 . Beyond this smaller domain the color is blue and the relative error is below 10 −12 . As we can see from Figure 4, the relative error is generally lower in the imaginary part and does not exceed ∼ 10 −11 .
It is commonly known that the accuracy of K (x, y) and L (x, y) functions tend to deteriorate when y decreases [2,16,21,22,28]. However, from Figures 3 and 4 we can see that the decrease of the parameter y does not deteriorate the accuracy. This is possible to achieve since in accordance with Equation (18) the number of the interpolation grid-points N gp increases with decreasing y. Therefore, this technique enables us to resolve efficiently this problem in computation of the complex error function w (z).

Run-Time Test
The run-time test was performed with MATLAB function files wTrap.m [37], fadsamp.m [31] and w2dom.m at equidistantly distributed 10 million points for parameter x at y = 10 −8 . The results of the run-time test are shown in Table 1.

Algorithm
Run It has been reported that both algorithms wTrap.m and fadsamp.m are highly accurate in computation [37]. In particular, the maximum values in relative errors for the algorithms wTrap.m and fadsamp.m are found to be ∼10 −15 and ∼10 −14 , respectively. However, the computational speed of these two algorithms differs depending on the range for the input parameter x. In particular, Table 1 shows that at smaller range for parameter x the function file wTrap.m performs computation faster than the function file fadsamp.m. However, as the range for parameter x increases, our algorithm fadsamp.m becomes faster. The run-time test also reveals that the algorithm w2dom.m is always faster regardless of the chosen range. This is particularly evident when the range for the input parameter x is smaller. Table 2 shows the results of the run-time test at 30 million equidistantly distributed points for parameter x at y = 10 −8 . As we can see, the algorithm w2dom.m remains more rapid proportionally at extended size of the input array as compared to the wTrap.m and fadsamp.m algorithms.

Algorithm
Run  The MATLAB code for the run-time test is provided in Appendix B. The scripts of function files wTrap.m and fadsamp.m can be accessed from the cited literature [30,37], respectively.