From: Bo =

From: Bo =?iso-8859-1?Q?Thid=E9?= <bt@lufsen.plasma.uu.se>
Organization: Dept of Space and Plasma Physics, Uppsala University
Subject: Re: Plasma dispersion function, C or Fortran
References: <87k0jl$6f0$1@jupiter.cs.uml.edu> <87pvai$63g$1@jupiter.cs.uml.edu>
Newsgroups: sci.physics.plasma

Paul Probert wrote:
>
> Hey Mirko,
>   I've got one that I wrote in IDL, which is mostly a rewrite of a
> fortran error function W(z) that I found on the web. Let me know and
> I'll e-mail it to you.
>
> Mirko Vukovic wrote:
> >
> > Can anyone send me a C or Fortran routine for computing
> > the plasma dispersion function, for real and complex
> > arguments. Or point me to a reference that explains
> > how to do it.
> >
> > I developed one based on the expansion given in Swanson's book
> > on plasma waves,
> > but I'm not sure whether it is valid for complex arguments.
> >
> > Thanks,
> >
> > Mirko
> >
> > Sent via Deja.com http://www.deja.com/
> > Before you buy.
>
> --
> _______________________________
> Paul Probert
> Associate Scientist
> The University of Wisconsin-Madison
> 608 263 3752
> fax 608 262 1267
> pprobert@facstaff.wisc.edu

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