faddeeva matlab,Faddeeva函数

Faddeva函数是用在计算等离子体色散函数中一个比较常用的函数，在上半平面有定义，下半平面解析延拓过去。与朗道阻尼联系密切，贴一个matlab的朗道阻尼的代码。

在matlab里 如果用fsolve 可以指定初始值为复数 like 1 + 0.1j 但是python的scipy没有这项功能 sympy 解不了色散函数中的erf ，numpy自然也很烦，这个问题有待于解决。还是赞美一下matlab并且赞美一下谢华生老师的计算等离子体物理。

直接贴github上找的faddeeva的matlab 代码了

function w = faddeeva(z,N)

% FADDEEVA Faddeeva function

% 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(-j*z)

%

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

%

% W = FADDEEVA(Z,N) can be used to explicitly specify the number of terms

% to truncate the expansion (see (13) in [1]). N = 16 is used as default.

%

% Example:

% x = linspace(-10,10,1001); [X,Y] = meshgrid(x,x);

% W = faddeeva(complex(X,Y));

% figure;

% subplot(121); imagesc(x,x,real(W)); axis xy square; caxis([-1 1]);

% title('re(faddeeva(z))'); xlabel('re(z)'); ylabel('im(z)');

% subplot(122); imagesc(x,x,imag(W)); axis xy square; caxis([-1 1]);

% title('im(faddeeva(z))'); xlabel('re(z)'); ylabel('im(z)');

%

% Reference:

% [1] 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

if nargin<2, N = []; end

if isempty(N), N = 16; end

w = zeros(size(z)); % initialize output

%%%%%

% for purely imaginary-valued inputs, use erf as is if z is real

idx = real(z)==0; %

w(idx) = exp(-z(idx).^2).*erfc(imag(z(idx)));

if all(idx), return; end

idx = ~idx;

%%%%%

% for complex-valued inputs

% make sure all points are in the upper half-plane (positive imag. values)

idx1 = idx & imag(z)<0;

z(idx1) = conj(z(idx1));

M = 2*N;

M2 = 2*M;

k = (-M+1:1:M-1)'; % M2 = no. of sampling points.

L = sqrt(N/sqrt(2)); % Optimal choice of L.

theta = k*pi/M;

t = L*tan(theta/2); % Variables theta and t.

f = exp(-t.^2).*(L^2+t.^2);

f = [0; f]; % Function to be transformed.

a = real(fft(fftshift(f)))/M2; % Coefficients of transform.

a = flipud(a(2:N+1)); % Reorder coefficients.

Z = (L+1i*z(idx))./(L-1i*z(idx));

p = polyval(a,Z); % Polynomial evaluation.

w(idx) = 2*p./(L-1i*z(idx)).^2 + (1/sqrt(pi))./(L-1i*z(idx)); % Evaluate w(z).

% convert the upper half-plane results to the lower half-plane if necesary

w(idx1) = conj(2*exp(-z(idx1).^2) - w(idx1));