Faddeeva function (FFT based)

Faddeeva function, aka plasma dispersion function

Description

W = FADDEEVA(Z) is the Faddeeva function, aka the plasma dispersion function, for each element of Z. The Faddeeva function is defined as:

     w(z) = exp(-z^2) erfc(-jz)

where erfc(x) is the complex complementary error function.

This code is a modified version of the code listed in:

J.A.C. Weideman, "Computation of the Complex Error Function," SIAM J. Numerical Analysis, pp. 1497-1518, No. 5, Vol. 31, Oct., 1994 Available Online: http://www.jstor.org/stable/2158232

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
Complex scaled complementary error function

• MATLAB release: MATLAB 7.7 (R2008b)

[Ulf Carlberg]

Very useful and efficient function. I think I found two small bugs though:

Line 25:
Change "z" to "z(idx)" in the exp-function.
w(idx) = exp(-z(idx).^2).*erfc(imag(z(idx)));

Line 54 (last line):
Take the conjugate of the result.
w(idx1) = conj(2*exp(-z(idx1).^2) - w(idx1));

I used this to visualize the function:

% Plot of the faddeeva function.
x = linspace(-10,10,1000); [xx,yy] = meshgrid(x,x);
ww = faddeeva(xx+1j*yy);
figure; imagesc(x,x,real(ww)); caxis([-1 1]); title('re(w)'); xlabel('re(z)'); ylabel('im(z)');
figure; imagesc(x,x,imag(ww)); caxis([-1 1]); title('im(w)'); xlabel('re(z)'); ylabel('im(z)');

[Tom Irvine]

This function helped me solve an acoustic problem for atmospheric sound propagation from a point source via both a direct path and ground reflection.

[Jean-Philippe Demers]

I used this function to simulate a Voigt profile. The function was fast and accurate.

I particularly liked the possibility to reduce the number of terms of the expansion, which further improve execution time, without affecting the accuracy of the solution.

[Woldemar Smith]

Overall the code is somehow OK. However it fails in accuracy at smaller y. For example, at z = x + 1i*y where x = 5.5 and y = 1e-4, the real part of this code gives 1.937525964550158e-006 while Fortran Algorithm 680 returns 1.966263443047126e-006. So the accuracy of this algorithm at y < 1e-4 is unacceptable - only one decimal digit matches. The code works at y > 1e-3. This is important because the real part of the Faddeeva function is the Voigt function used everywhere in science.

[Kesh Ikuma]

@Woldemar - The accuracy of the function is controlled by the (optional) second argument N, which defines the number of expansion terms. The default is N = 12, and for your example it is too low. If you increase N to say 100, you would get 1.96626344303036e-006, which is much closer to your Fortan solution.