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

The complex error function, also known as the Faddeeva function, is defined as [1, 2, 3, 4]

where $z=x+iy$ is a complex argument. The complex error function $w\left(z\right)$ 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]

such that

and

Both functions $w\left(z\right)$ and $W\left(z\right)$ are equal to each other on the upper half of the complex plane, when $y=\operatorname{Im}\left[z\right]>0$ [2]. Consequently, it follows that

Further we will imply that the parameter $y=\operatorname{Im}\left[z\right]$ is always greater than zero.

The real part of the complex error function $K\left(x,y\right)$ 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\left(x,y\right)$ is also used in many applications. For example, it can describe the spectral behavior of the index of refraction in various materials [14, 15].

Despite simple representations, the integrals (1), (2) and (3) do not have analytical solutions and must be computed numerically. Many approximations are available in the scientific literature [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Although this problem has been known for many decades, the derivation of new approximations for the complex error function $w\left(z\right)$ and the development of their efficient algorithms still remain an interesting topic [30, 31, 32, 33, 34, 35, 36, 37].

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.

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\geq 10^{-8}$ for the HITRAN [13] applications. However, this implementation has several limitations. In particular, there is a limitation $\left|x\right|\leq 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]

where the coefficients are [20]

such that

Consequently, the complex error function can also be approximated as [33]

Although the algorithms for the Voigt and complex error functions, built on Equations (3) and (3) 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 $\delta$ is the offset that was found experimentally to be $5\times 10^{3}$. The corresponding interpolation grid-points range is bounded by $\left[-r,r\right]$, where $r=35$ is the radius that was also found experimentally. The interpolation grid-points 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\geq 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 $\sqrt{y}$. The parameter $\delta$ in Equation (18) is to provide a sufficient number of the interpolation grid-points $N_{gp}$ within the internal domain $|x+iy|\leq 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

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 $\left|x+iy\right|=r$. The internal domain is situated within the circle while the external domain is situated outside of it.

All points within internal domain $\left|x+iy\right|\leq r$ are computed by the MATLAB 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 $\left|x+iy\right|>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\times 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.

In order to exclude the rounding and truncation errors in computation by using the most recent HITRAN database [13], the values of the $K\left(x,y\right)$ and $L\left(x,y\right)$ 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}$.

Figure 1 shows the absolute difference $\left|K_{ref.}\left(x,y\right)-K\left(x,y\right)\right|$, where $K_{ref.}\left(x,y\right)$ is a highly accurate reference, in the range $-5\leq x\leq 5$ at $y=10^{-8}$. As we can see, the absolute difference does not exceed $2.5\times 10^{-13}$.

Figure 2 shows the absolute difference $\left|L_{ref.}\left(x,y\right)-L\left(x,y\right)\right|$, where $L_{ref.}\left(x,y\right)$ is a highly accurate reference, in the range $-5\leq x\leq 5$ also at $y=10^{-8}$. We can see that for the imaginary part the absolute difference also does not exceed $2.5\times 10^{-13}$.

For more rigorous error analysis, we can apply the following relative errors

and

for the real and imaginary parts of the complex error function $w\left(z\right)$, respectively.

Figure 3 depicts the relative error for the real part of the complex error function $\operatorname{Re}\left[z\right]=K\left(x,y\right)$ in the domain $0\leq x\leq 50$ and $0\leq y\leq 10^{1.699}$, where $10^{1.699}\approx 50$. Inset in Figure 3 shows the smaller domain near the origin $0\leq x\leq 5.5$ and $0\leq y\leq 10^{0.74}$, where $10^{0.74}\approx 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 $\sim$$10^{-10}$.

Figure 4 illustrates the relative error for the imaginary part of the complex error function $\operatorname{Im}\left[z\right]=L\left(x,y\right)$ in the domain $0\leq x\leq 50$ and $0\leq y\leq 10^{1.699}$. Inset in Figure 4 shows the smaller domain near origin $0\leq x\leq 5.5$ and $0\leq y\leq 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 $\sim 10^{-11}$.

It is commonly known that the accuracy of $K\left(x,y\right)$ and $L\left(x,y\right)$ 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\left(z\right)$.

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.

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 $\sim$$10^{-15}$ and $\sim$$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.

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 [37, 30], respectively.

A new algorithm for the efficient computation of the Voigt/complex error function is shown. We propose a two-domain scheme where the number of interpolation grid-points $N_{gp}$ is dependent on parameter $y$. The error analysis shows that our MATLAB implementation meets the requirements for radiative transfer applications utilizing the HITRAN molecular spectroscopic database. The run-time test we performed shows that this MATLAB implementation can provide rapid computation especially at smaller ranges of parameter $x$.

This work is supported by the National Research Council Canada, Thoth Technology Inc., York University, Epic College of Technology and Epic Climate Green (ECG) Inc.

function FF = w2dom(x,y,opt)

% SYNOPSIS:
%          x   - array or scalar
%          y   - scalar
%          opt - option for the real and imaginary parts
%
% opt = 1 returns value(s) for the Voigt function, K(x,y)-function
% opt = 2 returns value(s) for the L(x,y)-function
% opt = 3 returns value(s) for the complex error function
%
% This code is primarily intended to work in the range y >= 1e-8 for the
% HITRAN applications for accelerated computation of the Voigt/complex
% error function.
%
% -------------------------------------------------------------------------
% Example:
%          x = linspace(-10,10,1e7); y = 1e-8; tic; w2dom(x,y); toc
% -------------------------------------------------------------------------
if ~isscalar(y)
    disp(’Input parameter y must be a scalar’)
    return % terminate computation if y is not a scalar
end

if nargin == 2
    opt = 3;
    disp(’Default value opt = 3 is assigned.’)
end

if opt ~= 1 && opt ~= 2 && opt ~= 3
    disp([’Wrong parameter opt = ’,num2str(opt),’! Use 1, 2 or 3.’])
    return
end

if y < 1e-8
    if opt == 1
        FF = real(fadsamp(x+1i*y));
    elseif opt == 2
        FF = imag(fadsamp(x+1i*y));
    else
        FF = fadsamp(x+1i*y);
    end
    return
end

FF = zeros(size(x));
radius = 35; % define radius
ind = abs(x + 1i*y) <= radius;

gp =[-radius,radius]; % define grid-points
if ~isempty(x(ind))

    offset = 5*1e3; % assign offset
    nump = 1/sqrt(y) + offset; % assign number of points
    % For better accuracy use for example nump = 2/sqrt(y) + 3*offset;

    gp = radius*(logspace(log10(1 + eps),log10(2),nump) ...
    - 1); % notice the log scale

    gp = [flip(-gp),gp];
end

switch opt
    case 1
        FF(ind) = real(internD(x(ind),y,gp));
        FF(~ind) = real(externD(x(~ind) + 1i*y));
    case 2
        FF(ind) = imag(internD(x(ind),y,gp));
        FF(~ind) = imag(externD(x(~ind) + 1i*y));
    otherwise
        FF(ind) = internD(x(ind),y,gp);
        FF(~ind) = externD(x(~ind) + 1i*y);
end

    function IntD = internD(x,y,gp) % internal domain
        IntD = interp1(gp,fadsamp(gp + 1i*y),x, ...
            ’spline’); % interpolated values
    end

    function ExtD = externD(z) % external domain

        num = 1:4; % define a row vector
        num = num/2;

        ExtD = num(end)./z; % start computing from the end
        for m = 1:length(num) - 1
            ExtD = num(end - m)./(z - ExtD);
        end
        ExtD = 1i/sqrt(pi)./(z - ExtD);
    end
end

clear
clc

% Table for run-time in seconds
tab{1,1}=’Algorithm’;
tab{1,2}=’-10 to 10’; % range 1
tab{1,3}=’-10^2 to 10^2’; % range 2
tab{1,4}=’-10^3 to 10^3’; % range~3

tab{2,1}=’wTrap’;
tab{3,1}=’fadsamp’;
tab{4,1}=’w2dom’;

y=1e-8; % smallest value for the parameter y
maxN=10; % max number of~cycles

for k=[1,3] % k is a factor
    for m=1:3 % three ranges
        x0=10^m; x=linspace(-x0,x0,k*1e7); % 1 and 3 are for 1e7 and 3*1e7 ...
                                           % points, respectively

        if k==1
            disp([’10 million points in range ’,num2str(m)])
        else
            disp([’30 million points in range ’,num2str(m)])
        end
        disp(’Computing, please wait!’)

        tic; for n=1:maxN; wTrap(x+1i*y,11); end; tab{2,m+1}=toc/maxN;
        tic; for n=1:maxN; fadsamp(x+1i*y); end; tab{3,m+1}=toc/maxN;
        tic; for n=1:maxN; w2dom(x,y); end; tab{4,m+1}=toc/maxN;

        clc
    end
    if k==1; tab1=tab; else tab2=tab; end % assign tab1 or tab2
end

disp(’Displaing Table 1:’)
disp(tab1)

disp(’Displaing Table 2:’)
disp(tab2)