! posted as Article 2347 of sci.physics.plasma:
! Subject: Re: Plasma dispersion function, C or Fortran
! Date: 3 Mar 2000 13:09:52 -0500
!
! Below, please find my Fortran77 code for the PDF.  It's written
! for HP-UX f77 but should be fairly portable.
! 
! Bo
! 
!
!  ^   Bo Thidé (Professor)               http://www.wavegroup.irfu.se/~bt
! |I|  Wave Group, Institute of Space Physics,  SE-755 91 Uppsala,  Sweden
! |R|  Office: Villavägen 3.  Phone: +46 018-4717269.  Fax: +46 018-554917 
!/|F|\ Mobile Phone: 0705-613670  Home Phone: 018-554184  Ham Call: SM5DFW
!~~U~~ E-Mail: mailto:bt@plasma.uu.se  Mobile mail: mailto:call-bt@irfu.se
!
!==========================================================================

      COMPLEX*16 FUNCTION Pdf(x,y)
*-----------------------------------------------------------------------
*
*  Computes the plasma dispersion function in the whole z-plane
*  by using the Voigt function
*
*    w(z) = exp(-z**2)erfc(-iz)
*
*  related to the plasma dispersion function  Pdf(z) = Z(z), 
*  usually defined as
*
*    Z(z) = 1/SQRT(Pi) * Integral(-Inf,+Inf) exp(-t**2)/(t-z) dt ,
*  
*  in the following way:
*
*    Z(z) = i*SQRT(Pi) * w(z)
*   
*  The algorithm used is based on algorithm 363 for calculating  w(z)
*  in the whole complex plane with an accuracy of at least 10 digits
*  after the decimal point given by:
*
*  W. Gautschi, "Complex Error Function", 
*  Comm. ACM 12, 635 (1969)
*
*  See also:
*
*  W. Gautschi, "Efficient computation of the complex error function",
*  SIAM J. Math. Anal. 7, 187 (1970)
* 
* (c) Copyright 1992 - bt@irfu.se (Bo Thide' Swedish Inst of Space Phys)
*-----------------------------------------------------------------------


      IMPLICIT REAL*8(A-Y)
      IMPLICIT COMPLEX*16(Z)
      INTEGER Capn,Nu,N,Np1
      LOGICAL B
      COMPLEX*16 Ci,CiTwoSqrtPi

      PARAMETER  ( SqrtPi = 1.7724 53850 90551 60272 98167D0)
      PARAMETER  ( CiTwoSqrtPi = (0.D0,3.544907701811032054596334) )
      PARAMETER  ( Ci = (0.D0,1.D0) )

      DATA Lambda /0.D0/

C     IF ((x.EQ.0.D0).AND.(y.EQ.0.D0)) THEN          !Trivial case
C       Pdf = DCMPLX(0.D0,SqrtPi)
C       RETURN
C     END IF
      z = DCMPLX(x,y)


      IF (y.GE.0.) THEN                              !Upper half z-plane
        IF (x.LT.0.D0) THEN                          !Second quadrant
C         WRITE(7,*) 'Second quadrant'
          Pdf=-DCONJG(Pdf(-x,y))                     !NB. Recursive call
C         WRITE(7,*) 'Recursion upper half plane'
          RETURN                                     !Ready
        END IF                                       !Now in 1st quadrant
      ELSE                                           !Lower half z-plane
        IF (x.LT.0.D0) THEN                          !Third quadrant
C         WRITE(7,*) 'Third quadrant'
          Pdf=CiTwoSqrtPi*ZEXP(-z*z) - Pdf(-x,y)     !NB. Recursive call
        ELSE                                         !Fourth quadrant
C         WRITE(7,*) 'Fourth quadrant'
          Pdf=CiTwoSqrtPi*ZEXP(-z*z) +
     &        DCONJG(Pdf(x,-y))                      !NB. Recursive call
        END IF                                       !Now in 3rd quadrant
C         WRITE(7,*) 'Recursion lower half plane'
        RETURN                                       !Ready
      END IF

C...Compute plasma dispersion function Z(z) in first quadrant in the z-plane

C     WRITE(7,*) 'First quadrant'

      IF ((y.LT.4.29D0).AND.(x.LT.5.33D0)) THEN   !Taylor series
C       WRITE(*,*) 'Taylor series',x,y
        s = (1.D0-y/4.29D0)*DSQRT(1.D0-x*x/28.41D0)
        h = 1.6D0*s
        h2 = h+h
        OneOverh2 = 1.D0/h2
        Capn = 6 + IDINT(23*s)
        Nu = 9 + IDINT(21*s)
      ELSE                                        !Gauss-Hermite quadrature
C       WRITE(*,*) 'Gauss-Hermite quadrature',x,y
        h = 0.D0
        Capn = 0
        Nu = 8
      END IF

      IF (h.GT.0.D0) Lambda = H2**Capn
      B = (h.EQ.0.D0).OR.(Lambda.EQ.0.D0)
      R1 = 0.D0
      R2 = 0.D0
      S1 = 0.D0
      S2 = 0.D0


      DO n = Nu, 0, -1
        Np1 = n + 1
        T1 = y + h + Np1*R1
        T2 = x - Np1*R2
        c = .5D0/(T1*T1 + T2*T2)
        R1 = c*T1
        R2 = c*T2
        IF ((n.LE.Capn).AND.(h.GT.0.D0)) THEN
          T1 = Lambda + S1
          S1 = R1*T1 - R2*S2
          S2 = R2*T1 + R1*S2
          Lambda = Lambda*OneOverh2
        END IF
C       WRITE(7,*) 't1,t2,s1,s2,lambda',t1,t2,s1,s2,lambda
      END DO

      IF (B) THEN
        Z_real = -2*R2
      ELSE
        Z_real = -2*S2
      END IF

      IF (y.EQ.0.D0) THEN
        Z_imag = SqrtPi*DEXP(-x*x)
      ELSE
        IF (B) THEN
          Z_imag = 2*R1
        ELSE
          Z_imag = 2*S1
        END IF
      END IF

      Pdf = DCMPLX(Z_real,Z_imag)
      RETURN
      END