Zplasma

 Plasma dispersion function and derivatives

0001 function [Z,Zp,Zpp] = Zplasma(x,N);
0002 %
0003 % Plasma dispersion function and derivatives
0004 %
0005 if imag(x) ~= 0,
0006     error('Not yet implemented for complex arguments')
0007 end
0008 %
0009 if nargin < 2
0010     N = 10000;
0011 end
0012 %
0013 xlim = 5;
0014 gamma = 0.845;
0015 %
0016 n = 0:N;
0017 %
0018 sx = size(x);
0019 ndim = length(sx);
0020 %
0021 signx = sign(x);
0022 x = abs(x);
0023 %
0024 Xx = shiftdim(repmat(x,[ones(1,ndim),N]),ndim);
0025 %
0026 Xrho_S = shiftdim(repmat(x,[ones(1,ndim),N+1]),ndim)/xlim;
0027 %
0028 Xn_S = repmat(n',[1,sx]);
0029 Xnt_S = N./(1 + (1 - gamma)*log((Xrho_S - gamma)/(1 - gamma)));  
0030 %
0031 Xt_S1 = xlim*Xrho_S.*Xn_S/N;
0032 Xt_S2 = xlim*Xn_S./Xnt_S;
0033 Xt_S3 = xlim*(gamma + (Xrho_S - gamma).*((1 - gamma)./(Xrho_S - gamma)).^((N - Xn_S)./(N - Xnt_S)));
0034 %
0035 Xmasklim_S = (Xrho_S <= 1);
0036 Xmasknt_S = (Xn_S <= Xnt_S);
0037 %
0038 Xt_S = zeros([N+1,sx]);
0039 Xt_S(Xmasklim_S) = Xt_S1(Xmasklim_S);
0040 Xt_S(~Xmasklim_S & Xmasknt_S) = Xt_S2(~Xmasklim_S & Xmasknt_S);
0041 Xt_S(~Xmasklim_S & ~Xmasknt_S) = Xt_S3(~Xmasklim_S & ~Xmasknt_S);
0042 %
0043 Xdt = diff(Xt_S);
0044 str = ',:';
0045 strtot = [];
0046 for idim = 1:ndim,
0047     strtot = [strtot,str];
0048 end
0049 eval(['Xt = (Xt_S(1:N',strtot,') + Xt_S(2:N+1',strtot,'))/2;'])
0050 %
0051 Xf = Xt.*exp(-Xt.^2)./sqrt(Xx.^2 - Xt.^2);
0052 %
0053 Z = -2*signx.*shiftdim(sum(Xf.*Xdt,1),1) + i*sqrt(pi)*exp(-x.^2);
0054 %
0055 Xfp = Xt.*exp(-Xt.^2).*(Xx./sqrt(Xx.^2 - Xt.^2) - 1);
0056 %
0057 Zp = -2*(exp(-x.^2) - 2*shiftdim(sum(Xfp.*Xdt,1),1)) - 2*i*sqrt(pi)*signx.*x.*exp(-x.^2);
0058 %
0059 %
0060 Xfpp = Xt.*exp(-Xt.^2).*(2*Xx + (Xt.^2 - 1)./Xx - (2*Xx.^2 - 1)./sqrt(Xx.^2 - Xt.^2));
0061 %
0062 Zpp = 4*signx.*(x.*exp(-x.^2) + shiftdim(sum(Xfpp.*Xdt,1),1)) - 2*i*sqrt(pi)*(1 + x.^2).*exp(-x.^2);
0063 %