% pldf.m: complex plasma dispersion function  
% GUISDAP v.1.60 96-05-27 Copyright Asko Huuskonen and Markku Lehtinen
%
%
% accuracy:  8 numbers in the whole complex plane
%            (main algorithm accurate to 1e-12, accuracy is
%            limited by subroutine cerfexp and numerical accuracy)
%
% method:    approximation described by Salzer in
%            Math. Tables Aids Comput., Vol. 5 (1951) pp. 67-70.
%            The same formula appears in
%            Abramovitch-Stegun p.299 formula 7.1.29.
%            (the exponentials in the formulas have been written
%            in a slightly modified way to get rid of unnecessary
%            overflows)
%                                           Markku Lehtinen 10.3.1979
%
% res=pldf(z)
%
  function res=pldf(z)
%
%  if max(size(z))==1
%    res=0;
%  else
%    res=zeros(z); 
%  end
  [nrow ncol]=size(z);
  res=zeros(nrow, ncol);
  j=sqrt(-1);
  for irow=1:nrow
    for icol=1:ncol
      z1=-j*z(irow,icol);
      x=real(z1);
      y=imag(z1);
      if(abs(2*x*y)<30000)
        cs=cos(2*x*y);
        sn=sin(2*x*y);
      else
        cs=1; sn=0;
      end
      pisqr=sqrt(pi);
      fn=max(1,floor(abs(2*y))-11):floor(abs(2*y))+11;
      term1=exp(-fn.*fn/4-y*y);
      if max(abs(term1))<1e-100, term1=zeros(fn); end
      term2=exp(-(fn/2-y).*(fn/2-y))/2;
      term3=exp(-(fn/2+y).*(fn/2+y))/2;
      factor=ones(size(fn))./(fn.*fn+4*x*x);
      sume=sum(factor.*(term1*cs-term2-term3));
      sumf=sum(factor.*(term1*2*x*sn+fn.*(term2-term3)));
      if (abs(x*y)<1e-4)
        e1=(-exp(-y*y)*x*y*y/2+2*x*sume)/pisqr;
        f1=(exp(-y*y)*y/2+sumf)/pisqr;
      else
         e1=(exp(-y*y)*(cs-1)/4/x+2*x*sume)/pisqr;
         f1=(exp(-y*y)*sn/4/x+sumf)/pisqr;
      end
      if abs(e1)<1e-100, e1=0; end
      z1=2*j*(e1+j*f1)+j*pisqr*exp(-y*y)*(cs+j*sn)*cerfexp(-x);
      res(irow,icol)=z1;
    end 
  end

% cerfexp.m: Approximation for the function  (1-erf(x))*exp(x**2), used in pldf
% GUISDAP v.1.60 96-05-27 Copyright Asko Huuskonen and Markku Lehtinen
%
%
% Approximation for the function  (1-erf(x))*exp(x**2) with
% accuracy 1e-8 for all x.  Method: power series for
% erf(x) for small x, asymptotic continued fractions expansion
% for large x (Abramovitz-Stegun formula 7.1.14)
%                                             Markku Lehtinen 2.7.1979
%
% res=cerfexp(x)
%
  function res=cerfexp(x)
%
  x1=abs(x);
  sqx=x*x;
     if (x1<1.) 
        tail=0.;
        term=sqx*sqx/10.;
        erf=1.-sqx/3.+term;
        for n=3:12 
           term=-term*sqx*(2*n-1.)/(2*n+1.)/n;
           tail=tail+term;
        end 
        erf=(erf+tail)*x*2./sqrt(pi);
        res=(1-erf)*exp(sqx);
      else
        dum=0;
        for n=(55:-1:1) 
           dum=n/2/(dum+x1);
         end 
         res=1/(dum+x1)/sqrt(pi);
         if (x<0), res=2*exp(sqx)-res; end 
      end
