A single-domain implementation of the Voigt/complex error function by vectorized interpolation
AA single-domain implementation of the Voigt/complex error functionby vectorized interpolation
S. M. Abrarov
1, 2 , B. M. Quine
1, 3 , R. Siddiqui
1, 2, 3 , and R. K. Jagpal
2, 31Dept. Earth and Space Science and Engineering, York University, 4700 Keele St., Canada, M3J 1P32Epic College of Technology, 5670 McAdam Rd., Mississauga, Canada, L4Z 1T23Dept. Physics and Astronomy, York University, 4700 Keele St., Toronto, Canada, M3J 1P3
August 31, 2019
Abstract
In this work we show how to perform a rapid computation ofthe Voigt/complex error over a single domain by vectorized inter-polation. This approach enables us to cover the entire set of theparameters x, y ∈ R required for the HITRAN-based spectroscopicapplications. The computational test reveals that within domains x ∈ [0 , ∩ y ∈ (cid:2) − , (cid:3) and x ∈ [0 , ∩ y ≥ − our al-gorithmic implementation is faster in computation by factors of about8 and 3, respectively, as compared to the fastest known C/C++ codefor the Voigt/complex error function. A rapid MATLAB code is pre-sented. Keywords: complex error function; Faddeeva function; Voigt func-tion; interpolation a r X i v : . [ m a t h . G M ] A ug Introduction
The complex error function, also commonly known as the Faddeeva function,can be defined as [1, 2, 3, 4] w ( z ) = e − z i √ π z (cid:90) e t dt , (1)where z = x + iy is the complex argument. The real part of the complexerror function (1), known as the Voigt function, can be represented as givenby [2, 5, 6] K ( x, y ) = 1 √ π ∞ (cid:90) e − t / e − yt cos ( xt ) dt (2a)that is proportional to the Voigt profile describing the spectral broadeningof atmospheric gas absorption or emission g V ( ν − ν , α L , α D ) = (cid:112) ln 2 /πα D K ( x, y ) , where ν is the frequency, ν is the frequency of the line center, α L and α D are the Lorentz and Doppler half widths at half maximum (HWHM),respectively, and x = √ ln 2 ν − ν α D , y = √ ln 2 α L α D . The imaginary part of the complex error function (1) has no specific name.Historically, it is denoted by L ( x, y ) and can be written in form [5, 6] L ( x, y ) = 1 √ π ∞ (cid:90) e − t / e − yt sin ( xt ) dt. (2b)None of the integrals above are analytically integrable in closed form.Therefore, these equations must be solved numerically. There are two mostimportant aspects that have to be taken into consideration for efficient algo-rithmic implementation of the integrals above in spectroscopic applicationsbased on HITRAN database [7]. Specifically, the latest versions of the HI-TRAN database provide spectroscopic data with 4 and more digits in floating2ormat. In order to exclude the truncation errors that inevitably occur inany line-by-line atmospheric modeling [8, 9, 10, 11, 12, 13], the algorithmicimplementation of the Voigt/complex error function should provide accuracyat least as close as possible to 10 − . This accuracy requirement is particularlyrelevant for the set of input numbers { x, y ∈ R : | x | + y ≤ } while it is notso much critical for the set { x, y ∈ R : | x | + y > } according to literature[14, 15]. Furthermore, the algorithm should be rapid as it may deal withmany millions of input numbers z in computation for the various radiativetransfer applications [15, 16].In order to satisfy these two criteria in computation of the integrals above,the complex plane is segmented into several domains. For example, Huml´ıˇcekproposed a rapid algorithm for computation of the Voigt function (2a) basedon rational functions by segmenting the complex plane into several domainssuch that each of them is computed by corresponding rational approxima-tion [17]. Kuntz proposed some modifications and showed efficiency of theHuml´ıˇcek algorithm with 4 domains [18]. The Huml´ıˇcek’s algorithm is rapidand provides accuracy 10 − . Later, Wells developed a FORTRAN code wherehe succeeded to improve accuracy by an order of the magnitude by makingsome modifications to the original Huml´ıˇcek’s algorithm [19]. Another in-teresting variation of the Huml´ıˇcek’s algorithm with improved accuracy waspresented by Imai et al. [20].It is very desirable to reduce wherever possible the number of the domainsto accelerate an algorithm since each area segmentation of the complex planeleads to run-time increase due to additional logical operations and subsequentsorting of the input numbers x and y .In our earlier publication [21] we have shown how to reduce the numberof domains to 2 by interpolation. In particular, we suggested two-domainscheme for rapid computation of the Voigt function that can be representedas K ( x, y ) ≈ interpolation, x + y (cid:54) a + b x a + b x + x , x + y > , (3) a = y/ (2 √ π ) + y / √ π ≈ . y + 0 . y b = y/ √ π ≈ . ya = 0 .
25 + y + y b = − y and a + b x a + b x + x = Re (cid:26) iz/ √ πz − / (cid:27) . Thus, the complex plane is segmented into two domains by an ellipse cen-tered at the origin with semi-major and semi-minor axises equal to 27 and15, respectively (see inset in Fig. 2 from our paper [21]). Inside the ellipse(computationally difficult internal domain) we apply interpolation while out-side the ellipse (computationally simple external domain) we apply a simplerational approximation of low order. This scheme shows high efficiency par-ticularly when x = { x , x , x , . . . } is a vector and y is a scalar. In fact,vectorized x = { x , x , x , . . . } and scalar y is a quite common technique inradiative transfer applications [22, 23, 24].Recently Schreier reported a two-domain scheme (see equation (12) in[15]) w ( z ) = K ( x, y ) + iL ( x, y ) ≈ (cid:80) n − k =0 α k z k (cid:80) n(cid:96) =0 β (cid:96) z (cid:96) , | x | + y ≤ iz/ √ πz − / , | x | + y > , (4) where n is assumed to be an even integer, α k and β (cid:96) are the expansion co-efficients that can be readily generated by Computer Algebra System (CAS)supporting symbolic programming. It is suggested that n = 20 in the sin-gle quotient is sufficient for line-by-line calculations in atmospheric modeling[15].In general, for arbitrary n (even or odd) representation of the Huml´ıˇcek’sapproximation as a single quotient can be made by using notation (cid:100) . . . (cid:101) forthe ceiling function as w ( z ) ≈ n (cid:88) k =1 (cid:18) γ k + iθ k z − x k + iδ − γ k − iθ k z + x k + iδ (cid:19) = (cid:80) (cid:100) n/ (cid:101)− k =0 α k z k (cid:80) (cid:100) n/ (cid:101) (cid:96) =0 β (cid:96) z (cid:96) , (5)4here γ k = − π ω k e δ sin (2 x k δ ) ,θ k = 1 π ω k e δ cos (2 x k δ ) ,x k are roots of the Hermite polynomial H n ( x ) of degree n that can be definedby recurrence relations [25] H n +1 = 2 xH n ( x ) − nH n − ( x ) , H ( x ) = 1 , H ( x ) = 2 x,δ is a fitting parameter that at n = 20 can be taken as 1 .
55 [15] and ω k = 2 n − n ! √ πn H n − ( x k )are weights of the Hermite polynomial H n ( x ) of degree n [26].It is interesting to note that if integer n is even and roots are given inascending order such that x k − < x k , then number of the summation termsin the Huml´ıˇcek’s approximation can be reduced by a factor of two (compareequations (27a) and (26b) from [8]) w ( z ) ≈ n/ (cid:88) k =1 (cid:18) γ k + iθ k z − x k + iδ − γ k − iθ k z + x k + iδ (cid:19) . The single quotient reformulation shown in equations (4) and (5) is in-teresting. However, by performing computational test we found empiricallythat deterioration of accuracy with decreasing y is especially inherent to thesingle quotient reformulation of the Huml´ıˇcek’s approximation (5) (deterio-ration of accuracy with decreasing y is a common problem in computation ofthe Voigt/complex error function [2, 27, 28]). Although a multiple precisionarithmetic may be used to resolve this problem, it needs a special package(see for example [29]) that may affect computational speed and makes theMATLAB code inconvenient in practical applications.In plasma physics of rarefied gases or at low atmospheric pressure thattakes place in stratosphere, mesosphere and thermoshpere of the Earth,where the Doppler broadening considerably predominates over the Lorentzbroadening, the value of y dependent on the pressure and temperature may berelatively close to zero. This is particularly important as the latest versionsof the HITRAN supply parameters for high temperatures almost reaching50000 K and there is a tendency that it will be increased in future. There-fore, it would be very desirable to develop a rapid algorithm that can sustainthe required accuracy for input parameter y ≥ − [15].One of the possible ways to overcome these problems is to segment thecomplex plane as a narrow band along x -axis and to use appropriate approx-imation for smaller Im [ z ] = y (see for example C/C++ code [30]). However,this decelerates computation as a result of additional segmentation.Consider a complete version of the two-domain scheme (3) that includesboth, the real and imaginary parts, as follows w ( z ) = K ( x, y ) + iL ( x, y ) ≈ interpolation , x + y (cid:54) iz/ √ πz − / , x + y > . (6) The run-time test we performed shows that for the set of input numbers { x, y ∈ R : | x | + y ≤ } this two-domain scheme is faster in performance bya factor about 2 as compared to that of reported in [15] (we used the MAT-LAB codes built on approximations (4) and (6)). This is possible to achievesince interpolation utilizes a simple cubic spline instead of a rational functionof high order. Therefore, this fact strongly motivated us to develop furtheran algorithm based on a vectorized interpolation for rapid computation ofthe Voigt/complex error function with accuracy that meets the requirementfor the HITRAN spectroscopic applications [7].In this work we propose a new method of algorithmic implementation forrapid computation of the integrals (1), (2a) and (2b) by vectorized interpola-tion that enables us to employ just a single domain. This approach providesboth, the rapid computation and accuracy that meets the requirement forthe HITRAN spectroscopic applications. To the best of our knowledge, thismethod of computation of the Voigt/complex error function is new and hasnever been reported in scientific literature. Computation of the Voigt function (1) by interpolation with help of thelookup table was first reported in the work [31] (see also [32]). However, thelookup table involved in this method requires two dimensional interpolatinggrid-points that need extra memory and time to handle large-size data during6omputation to obtain a reasonable accuracy. In contrast, the vectorized (onedimensional) approach, where the interpolating grid points are computeddynamically at x = { x , x , x , . . . } and given fixed value y [21] (rather thanpicked up from the lookup table [31, 32]) is considerably advantageous ininterpolation due to significantly smaller quantity of the interpolating grid-points stored in computer memory.Let us consider two-domain scheme shown by equation (6) more closely.We note that semi-major and semi-minor axises may be chosen arbitrarilydepending on algorithmic implementations (see for example [33, 34]). Sincesemi-major and semi-minor axises may be flexible in magnitude, we maysuggest to increase them to cover the entire set x, y ∈ R for the HITRAN ap-plications such that we could exclude completely the rational approximationof low order in equation (6) from the consideration.Previously it was reported that the HITRAN spectroscopic database re-quires a domain 0 ≤ | x | < − < y < [19]. Consequently,within the I-st quadrant a single domain should be large enough to coverrectangular area 40000 × y into considerationwould be preferable for practical applications to account for low pressure andhigh temperature of the HITRAN gases. In our approach, the single domainrepresents an area | x | (cid:54) y sinceit is a scalar. Formally stating, this single domain is inclosed by a reshapedellipse in such a way that x + lim n →∞ y n (cid:54) . In our work [21] we noted that the interpolating grid-points may not benecessarily spaced equidistantly. This provides a significant advantage sincewe can minimize the number of the interpolating grid-points by putting themdenser in more curved subintervals and sparser in more linear subintervals.Particularly, we separated interval [ − , x -axis into subin-tervals as it is shown in the Table 1. Initially each subinterval contains anumber of the interpolating grid-points multiple to 100 accumulating 1600interpolating grid-points in total (see command lines below < function MF= mainF(x,y,opt) > in the MATLAB code shown in Appendix A). However,after exclusion of repeating values the total number of the interpolating grid-points decreases from 1600 to 197 + 398 + 4 ×
198 + 200 = 1587 (see Table7). This small quantity of the interpolating grid-points needs a very mi-nor amount of computer memory. Therefore, it can be easily handled bypractically any modern computer. We cannot infer that such a distribu-tion provides the smallest number of the interpolating grid-points. Perhaps,this number of the interpolating grid-points can be reduced significantly byfurther optimization.
Subintervals Number of IGP ( − . , .
5) 1 ×
197 = 197( − . , − . ∪ [2 . , .
5) 2 ×
199 = 398( − , − . ∪ [5 . ,
15) 2 ×
99 = 198( − , − ∪ [15 , ×
99 = 198( − , − ∪ [100 , ×
99 = 198( − , − ∪ [1000 , ×
99 = 198[ − , − ∪ [10000 , ×
100 = 200
Table 1.
Number of interpolating grid-points (IGP) in subintervals.
If hypothetically there could be some input values | x | greater than 50000,then it is very easy to extend the region along x -axis if required. For exam-ple, it is sufficient to include only 300 additional interpolating grid-points toextend the range for | x | up to 100000; it is not necessary to include manyinterpolating grid-points since the functions K ( x, y ) and L ( x, y ) becomenearly linear as | x | increases.For proper interpolation the accuracy of the computation should be twoorders of the magnitude better than 10 − . Generally, any highly accurate al-gorithm can be used to compute 1587 interpolating grid-points. We appliedthe highly accurate MATLAB function fadsamp.m shown in our recent paper[35] as a subroutine for this purpose (see a brief description of this method inAppendix B). Alternatively, a MATLAB function file that is highly accurateand suitable for computation of the interpolating grid-points can be found in[28] (these two codes can also be downloaded from the Matlab Central web-sites [36] and [37], respectively). Highly accurate C/C++ implementationby Johnson [30] with MEX plugins for MATLAB [38] can also be used for asubroutine that can be invoked from the MATLAB environment to computethese 1587 interpolating grid-points.It should be noted that this methodology can also be generalized further8o other functions (including spectral line profiles). For example, our prelim-inary results demonstrate that this methodology is also applicable for rapidand quite accurate computation of the Spectrally Integrated Voigt Function(SIVF) that enables us to perform line-by-line atmospheric modeling at re-duced spectral resolution [13]. In order to perform error analysis define the relative errors for the real∆ Re = (cid:12)(cid:12)(cid:12)(cid:12) K ref. ( x, y ) − K ( x, y ) K ref. ( x, y ) (cid:12)(cid:12)(cid:12)(cid:12) and imaginary parts ∆ Im = (cid:12)(cid:12)(cid:12)(cid:12) L ref. ( x, y ) − L ( x, y ) L ref. ( x, y ) (cid:12)(cid:12)(cid:12)(cid:12) of the complex error function (1), where K ref. ( x, y ) and L ref. ( x, y ) are thehighly accurate reference values that can be readily obtained by using theCAS. Fig. 1.
Logarithms of relative error for the real a) and imaginary b)parts of the complex error function over the area x ∈ [0 , ∩ y ∈ [10 − , − ] computed by vectorized interpolation. Insets show thesubareas with worst accuracies for the real a) and imaginary b) parts. Figures 1a and 1b show the relative errors for the real and imaginaryparts over the area x ∈ [0 ,
15] and y ∈ [10 − , − ] computed by vectorizedinterpolation. The largest relative errors over this area for the real and9maginary parts are 1 . × − and 7 . × − , respectively. Thus,our algorithm satisfies the accuracy criterion 10 − for the required range y ≥ − . In comparison, the computational test we perform reveals thatapproximation (4) sustains the required accuracy only at y ≥ − .Figures 2a and 2b show the relative errors for the real and imaginary partsover the area x ∈ [0 ,
15] and y ∈ (cid:2) − , (cid:3) . The largest relative errors overthis area for the real and imaginary parts are 2 . × − and 7 . × − ,respectively.Thus, from Figs. 1a, 1b and 2a, 2b we can conclude that our algorithmsatisfies accuracy requirement for the HITRAN-based applications and effec-tively resolves the problem that is inherent to the single quotient shown inequations (4) and (5) as well as many other approximations [2, 27].In order to obtain most objective results for the run-time test we usedan independently written code for the Voigt/complex error function. Specifi-cally, the run-time test has been performed by comparing our MATLAB codeshown in Appendix A with C/C++ implementation that was developed byJohnson [30, 38]. This implementation is known to be the fastest C/C++program. It represents a modified Algorithm 680 [33, 34] with inclusion ofthe Salzer’s approximations [39] (see also [40]). As the computation com-plexity prevails at smaller values of the parameter y , we imply that it is closeto zero, say y = 10 − . The detailed description of how to run the programsand perform the time execution test is shown in Appendix C.The run-time test shows that for 10 million input values within the mostimportant area such that { x, y ∈ R : | x + iy | ≤ } , the MATLAB code isalmost 8 times faster than the C/C++ implementation while for the en-tire domain that is required for the HITRAN spectroscopic applications theMATLAB code is faster than the C/C++ implementation by a factor about3. The MATLAB is one of the fastest scientific languages in computation.However, it is generally slower than C/C++. The more rapid computationhas been achieved because of two main reasons. The function Faddeeva.cc is unnecessarily complicated as it utilizes a large number of domains thatdue to multiple logical operations and sorting of the input numbers x and y decelerate computation. Furthermore, our approach de facto performs1587 actual computations only as the remaining is just an interpolation. Bychoosing different options for interpolation, we found experimentally that theMATLAB built-in method ’spline’ provides the best performance. All theseresults can be readily confirmed by running the codes provided in Appendices10 ig. 2. Logarithms of relative error for the real a) and imaginary b)parts of the complex error function over the area x ∈ [0 , ∩ y ∈ [10 − , computed by vectorized interpolation. A and C.
We propose a new single-domain vectorized interpolation method for rapidcomputation of the Voigt/complex error function (1) that enables us toachieve required accuracy for the HITRAN-based spectroscopic applications.The computational test we performed reveals that within intervals x ∈ [0 , ∩ y ∈ [10 − ,
15] and x ∈ [0 , ∩ y ≥ − our algorithmic imple-mentation is faster in computation by factors of about 8 and 3, respectively,as compared to the fastest known C/C++ code for the Voigt/complex errorfunction. Acknowledgments
This work is supported by National Research Council Canada, Thoth Tech-nology Inc., York University, Epic College of Technology and Epic ClimateGreen (ECG) Inc. The authors wish to thank principal developer of theMODTRAN Dr. Alexander Berk for constructive discussions.11 ppendix A function vecFF = vecfadf(x,y,opt)% This function file computes the Voigt/complex error function, also known% as the Faddeeva function, providing rapid computation at required% accuracy for the HITRAN-based radiative transfer application.%% SYNOPSIS:% x - row or column vector% y - scalar% opt - option for the real and imaginary parts%% The code is written by Sanjar M. Abrarov, Brendan M. Quine, Rehan% Siddiqui and Rajinder K. Jagpal, York University, Canada, May 2019.bound = 5*1e4; % default bound to cover the HITRAN domainnum = 1e2; % common number for grid-points in interpolationif max(size(y)) ~= 1disp( ' Parameter y must be a scalar! ' )returnelseif ~isvector(x)disp( ' Parameter x is not a vector! ' )endif max(abs(x)) > bound || abs(y) < 1e-8disp( ' x or y is beyond HITRAN range! Computation is terminated. ' )returnendif nargin == 2opt = 3;disp( ' Default value opt = 3 is assigned. ' )endif opt ~= 1 && opt ~= 2 && opt ~= 3disp([ ' Wrong parameter opt = ' ,num2str(opt), ' ! Use either 1, 2 or 3. ' ])returnendif y >= 0vecFF = mainF(x,y,opt); % upper half-planeelse ecFF = mainF(x,-y,opt); % lower half-planevecFF = conj(2*exp(-(x + 1i*y).^2) - vecFF);end function MF = mainF(x,y,opt)% Forming non-equidistantly spaced interpolating grid-points (IGP)IGP = linspace(0,2.5,num);IGP = [IGP,linspace(2.5,5.5,100 + num)]; % 100 more grid-pointsIGP = [IGP,linspace(5.5,15,num)];IGP = [IGP,linspace(15,100,num)];IGP = [IGP,linspace(100,1000,num)];IGP = [IGP,linspace(1000,10000,num)];IGP = [IGP,linspace(10000,bound,num)];IGP = [-(flip(IGP)),IGP];IGP = unique(IGP); % exclude repeated valuesMF = interp1(IGP,fadsamp(IGP + 1i*y),x, ' spline ' ); % call ...% external function
75 and the expansion coefficients A m = √ π ( m − / M h N (cid:88) n = − N e ς / − n h sin (cid:18) π ( m − /
2) ( nh + ς/ M h (cid:19) ,B m = − iM √ π N (cid:88) n = − N e ς / − n h cos (cid:18) π ( m − /
2) ( nh + ς/ M h (cid:19) C m = π ( m − / M h .
Similar to the Huml´ıˇcek’s approximation (5), the series expansion (A.2)also deteriorates in accuracy with decreasing y . We have shown how toovercome this problem by transformation of equation (A.2) into followingform w ( z ) ≈ e − z + z M +2 (cid:88) m =1 κ m − λ m z µ m − ν m z + z , (A.3)where κ m = B m (cid:34) C m − (cid:18) ς (cid:19) (cid:35) + iA m ς = B m (cid:34)(cid:18) π ( m − / M h (cid:19) − (cid:16) ς (cid:17) (cid:35) + iA m ς,λ m = B m ,µ m = C m + C m ς ς
16 = (cid:34)(cid:18) π ( m − / M h (cid:19) + (cid:16) ς (cid:17) (cid:35) and ν m = 2 C m − ς (cid:18) π ( m − / M h (cid:19) − ς . The third equation that is used in the function file fadsamp.m is theLaplace continued fraction [4, 3, 33] given by w ( z ) ≈ ( i/ √ π ) z − / z − z − / z − z − / z − z − / z − z − / z − z − / z , (A.4)Equations (A.2) and (A.3) cover the domains { x, y ∈ R : | x + iy | ≤ } \ { y < . | x |} and { x, y ∈ R : | x + iy | ≤ } ∩ { y < . | x |} , respectively, while equation (A.3) covers the domain { x, y ∈ R : | x + iy | > } (see Fig. 1 in [35]). This three-domain scheme excludes all poles in compu-tation and provides largest relative error ∼ − only.15 ppendix C The C/C++ function file
Faddeeva.cc and its header
Faddeeva.hh can bedownloaded from the website [30]. To run the program the following linescan be added to the end of the
Faddeeva.cc function file /**These command lines can be added at the end of the sourcefile ' Faddeeva.cc ' .*/ ' minNum ' to ' maxNum ' .*/double fRand(double minNum, double maxNum){double fMin = minNum, fMax = maxNum;double f = (double)rand() / RAND_MAX;return fMin + f * (fMax - fMin);};int main(void){double(epsVal) = 1E-6; // epsilon value for accuracyfor(int k = 0; k < 1e7; k++){ // execute the computation ...// 10 million times with random numbers 0 < x < 15 and positive y << 1.FADDEEVA(w)(C(fRand(0,15),1E-5),epsVal);};return 0;};
16s we can see from the body of the function < int main(void) > , this pro-gram computes 10 million random digits of x at fixed value of y = 10 − .The epsilon value that determines accuracy of computation is taken to be10 − . The execution time on a typical laptop computer (we used Intel(R)Core(TM) i3-7020U CPU @ 2.30GHz, 8GB RAM) takes about 14 seconds.It is interesting to note that with < double(espVal) = 1E-12; > the execu-tions time increases to 17 seconds. As a subsequent step, at < double(espVal)= 0; > we should expect the highest accuracy and, therefore, the longest exe-cution time. Instead, however, the program becomes significantly faster andit takes about 8 . x and y at < double(espVal) = 0; > . Therefore, we did not take this option as areference.The following is the MATLAB code that also computes 10 million randomnumbers x ∈ [0 ,
15] at fixed y = 10 − . The execution time is about 1 . % ******************************************************************x = 15*rand(1e7,1); % this generates 10 million random numbers ...% in the interval 0 < x < 15y = 1e-5; % y is a scalartic;vecfadf(x,y,3);toc % this shows the run-time% ****************************************************************** In order to compare the execution times for the interval x ∈ [0 , FADDEEVA(w)(C(fRand(0,50000),1E-5),epsVal); and x = 50000*rand(1e7,1); % this generates 10 million random numbers ...% in the interval 0 < x < 50000 respectively. The execution times for the C/C++ and MATLAB imple-mentations for this case become about 3 and about 1 seconds, respectively.Therefore, the MATLAB code is nearly 3 times faster.Execution times can be decreased by more than an order of the magnitudeon a more powerful computer of the latest generation. However, we shouldanticipate that these ratios 8 and 3 will remain same.17 eferences [1] V.N. Faddeyeva, and N.M. Terent’ev, Tables of the probability integral w ( z ) = e − z (cid:16) i √ π (cid:82) z e t dt (cid:17) for complex argument. Pergamon Press,Oxford, 1961.[2] B.H. Armstrong, Spectrum line profiles: The Voigt function, J. Quant.Spectrosc. Radiat. Transfer, 7 (1) (1967) 61-88. https://doi.org/10.1016/0022-4073(67)90057-X [3] W. Gautschi, Efficient computation of the complex error function. SIAMJ. Numer. Anal., 7 (1) (1970) 187-198. https://doi.org/10.1137/0707012 [4] M. Abramowitz and I.A. Stegun. Error function and Fresnel integrals.Handbook of mathematical functions with formulas, graphs, and mathe-matical tables. 9th ed. New York 1972, 297-309.[5] H.M. Srivastava and E.A. Miller, A unified presentation of the Voigtfunctions. Astrophys. Space Sci., 135 (1) (1987) 111118. https://doi.org/10.1007/bf00644466 [6] H.M. Srivastava and M.P. Chen, Some unified presentations of the Voigtfunctions. Astrophys. Space Sci., 192 (1) (1992) 63-74. https://doi.org/10.1007/BF00653260 [7] C. Hill, I.E. Gordon, R.V. Kochanov, L. Barrett, J.S. Wilzewski andL.S. Rothman, HITRAN online : An online interface and the flexible rep-resentation of spectroscopic data in the HITRAN database, J. Quant.Spectrosc. Radiat. Transfer, 177 (2016) 4-14. https://doi.org/10.1016/j.jqsrt.2015.12.012 [8] A. Berk and F. Hawes, Validation of MODTRAN ® https://doi.org/10.1016/j.jqsrt.2017.03.004 [9] D. Pliutau and K. Roslyakov, Bytran - | - spectral calculations for portabledevices using the HITRAN database, Earth Sci. Inform., 10 (3) (2017)395-404. https://doi.org/10.1007/s12145-017-0288-4 https://doi.org/10.1504/IJSPACESE.2015.075911 [11] R. Siddiqui, R. Jagpal and B.M. Quine, Short wave upwelling radia-tive flux (SWupRF) within near infrared (NIR) wavelength bands of O , H O , CO , and CH by Argus 1000 along with GENSPECT line by lineradiative transfer model, Canad. J. Remote Sens., 43 (4) (2017) 330-344. https://doi.org/10.1080/07038992.2017.1346467 [12] R.K. Jagpal, B.M. Quine, H. Chesser, S. Abrarov and R. Lee, Cali-bration and in-orbit performance of the Argus 1000 spectrometer - theCanadian pollution monitor, J. Appl. Remote Sens., 4 (1) (2010) 049501. https://doi.org/10.1117/1.3302405 [13] B.M. Quine and S.M. Abrarov, Application of the spectrally integratedVoigt function to line-by-line radiative transfer modelling, J. Quant. Spec-trosc. Radiat. Transfer, 127 (2013) 37-48. http://dx.doi.org/10.1016/j.jqsrt.2013.04.020 [14] F. Schreier, Optimized implementations of rational approximations forthe Voigt and complex error function, J. Quant. Spectrosc. Radiat. Trans-fer, 112 (6) (2011) 1010-1025. https://doi.org/10.1016/j.jqsrt.2010.12.010 [15] F. Schreier, The Voigt and complex error function: Huml´ıˇcek’s rationalapproximation generalized, Mon. Not. Roy. Astron. Soc., 479 (2018) 3068-3075. https://doi.org/10.1093/mnras/sty1680 [16] S.L. Grimm and K. Heng, HELIOS-K: An ultrafast, open-source opacitycalculator for radiative transfer, Astrophys. J., 808:182 (2015). https://doi.org/10.1088/0004-637X/808/2/182 [17] J. Huml´ıˇcek, Optimized computation of the Voigt and complex prob-ability functions, J. Quant. Spectrosc. Radiat. Transfer, 27 (4) (1982)437-444. https://doi.org/10.1016/0022-4073(82)90078-4 https://doi.org/10.1016/S0022-4073(96)00162-8 [19] R.J. Wells, Rapid approximation to the Voigt/Faddeeva function and itsderivatives, J. Quant. Spectrosc. Radiat. Transfer, 62 (1) (1999) 29-48. https://doi.org/10.1016/S0022-4073(97)00231-8 [20] K. Imai, M. Suzuki and C. Takahashi, Evaluation of Voigt algorithmsfor the ISS/JEM/SMILES L2 data processing system, Adv. Space Res.,45 (5) (2010) 669-675. https://doi.org/10.1016/j.asr.2009.11.005 [21] S.M. Abrarov, B.M. Quine and R.K. Jagpal, A simple interpolatingalgorithm for the rapid and accurate calculation of the Voigt function, J.Quant. Spectrosc. Radiat. Transfer, 110 (67) (2009) 376-383. https://doi.org/10.1016/j.jqsrt.2009.01.003 [22] A.E. Lynas-Gray, VOIGTL a fast subroutine for Voigt function eval-uation on vector processors, Comput. Phys. Commun., 75 (1-2) (1993)135-142. https://doi.org/10.1016/0010-4655(93)90171-8 [23] K.L. Letchworth and D.C. Benner, Rapid and accurate calculation ofthe Voigt function, J. Quant. Spectrosc. Radiat. Transfer, 107 (1) (2007)173-192. https://doi.org/10.1016/j.jqsrt.2007.01.052 [24] F. Schreier and D. Kohlert, Optimized implementations of rational ap-proximations a case study on the Voigt and complex error function,Comput. Phys. Commun., 179 (7) (2008) 457-465. https://doi.org/10.1016/j.cpc.2008.04.012 [25] E.W. Weisstein, Hermite polynomial. http://mathworld.wolfram.com/HermitePolynomial.html [26] E.W. Weisstein, Hermite–Gauss quadrature. http://mathworld.wolfram.com/Hermite-GaussQuadrature.html a where the calculationis notoriously difficult, Amer. J. Anal. Chem., 4 (12) (2013) 725-731. https://doi.org/10.4236/ajac.2013.412087 [28] S.M. Abrarov and B.M. Quine, Efficient algorithmic implementation ofthe Voigt/complex error function based on exponential series approxima-tion, Appl. Math. Comput., 218 (5) (2011) 1894-1902. https://doi.org/10.1016/j.amc.2011.06.072 [29] D. Tsarapkina and D.J. Jeffrey, Exploring rounding errors in Matlabusing extended precision, Proc. Comput. Sci., 29 (2014) 1423-1432. https://doi.org/10.1016/j.procs.2014.05.129 [30] S.G. Johnson, Faddeeva package. http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package [31] J.R. Drummond and M. Steckner, Voigt-function evaluation using a two-dimensional interpolation scheme, J. Quant. Spectrosc. Radiat. Transfer,34 (6) (1985) 517-521. https://doi.org/10.1016/0022-4073(85)90145-1 [32] L. Sparks, Efficient line-by-line calculation of absorption coefficients tohigh numerical accuracy, J. Quant. Spectrosc. Radiat. Transfer, 57 (5)(1997) 631-650. https://doi.org/10.1016/S0022-4073(96)00154-9 [33] G.P.M. Poppe and C.M.J. Wijers, More efficient computation of thecomplex error function. ACM Transact. Math. Software, 16 (1) (1990)38-46. https://doi.org/10.1145/77626.77629 [34] G.P.M. Poppe and C.M.J. Wijers, Algorithm 680: evaluation of thecomplex error function. ACM Transact. Math. Software, 16 (1) (1990)47. https://doi.org/10.1145/77626.77630 [35] S.M. Abrarov, B.M. Quine and R.K. Jagpal, A sampling-based approx-imation of the complex error function and its implementation withoutpoles, Appl. Num. Math., 129 (2018) 181-191. https://doi.org/10.1016/j.apnum.2018.03.009 https://doi.org/10.1016/j.amc.2017.10.032 [40] H.E. Salzer, Formulas for calculating the error function of a complexvariable, Math. Tables Aids Comput. 5 (34) (1951) 67-70. https://doi.org/10.2307/2002163 [41] S.M. Abrarov and B.M. Quine, Sampling by incomplete cosine expansionof the sinc function: Application to the Voigt/complex error function,Appl. Math. Comput., 258 (2015) 425-435. https://doi.org/10.1016/j.amc.2015.01.072 [42] S.M. Abrarov and B.M. Quine, A rational approximation for efficientcomputation of the Voigt function in quantitative spectroscopy, J. Math.Research, 7 (2) (2015) 163-174. https://doi.org/10.5539/jmr.v7n2p163 [43] W.B. Gearhart and H.S. Shultz, The function sin xx , College Math. J.,21 (2) (1990) 90-99. http://dx.doi.org/10.1080/07468342.1990.11973290 [44] E.W. Weisstein, Sinc function. http://mathworld.wolfram.com/SincFunction.htmlhttp://mathworld.wolfram.com/SincFunction.html