C      SUBROUTINE WEIDEMAN(XI, YI, U, V, FLAG)
C  GIVEN A COMPLEX NUMBER Z = (XI,YI), THIS SUBROUTINE COMPUTES
C  THE VALUE OF THE FADDEEVA-FUNCTION W(Z) = EXP(-Z**2)*ERFC(-I*Z),
C  WHERE ERFC IS THE COMPLEX COMPLEMENTARY ERROR-FUNCTION AND I
C  MEANS SQRT(-1).THE ALGORITHM IS BASED ON THE PAPER PUBLISHED ON
C  SIAM JOURNAL ON NUMERICAL ANALYSIS VOL.31 NO.5 PP.1497-1518,OCTOBER
C  1994 BY J.A.C. WEIDEMAN.THE USER CAN CHOOSE BETWEEN N=4, N=8 AND N=16.
      SUBROUTINE WOFZWH2(Z1,Z2,RES,M,FLAG)
      IMPLICIT NONE
      EXTERNAL CPOLYEV2
      DOUBLE PRECISION CFFT(16), L,YI1,YI2,PI
      INTEGER N,M
      PARAMETER (PI = 3.141592653589793)
      LOGICAL FLAG
      DOUBLE COMPLEX Z1,I,Z2,A1,B1,A2,B2,P1,P2,W,RES,ZZ1,ZZ2,
     *     ZR1,ZR2
C     THE N COEFFICIENTS OF FAST FOURIER TRANSFORM
C     COEFFICIENTS TO USE FOR N=4
C      DATA CFFT(1) / -0.036387778131343 /
C      DATA CFFT(2) / -0.009515238439985 /
C      DATA CFFT(3) / 0.263633942774273 /
C      DATA CFFT(4) / 0.715253452379754 /
C     COEFFICIENTS TO USE FOR N=8
C      DATA CFFT(1) / 0.001422984804155 /
C      DATA CFFT(2) / 0.001514258764391 /
C      DATA CFFT(3) / -0.008825132847676 /
C      DATA CFFT(4) / -0.017753679350031 /
C      DATA CFFT(5) / 0.044594093387920 /
C      DATA CFFT(6) / 0.280930761779060 /
C      DATA CFFT(7) / 0.706286816137253 /
C      DATA CFFT(8) / 1.149736350667053 /
C     COEFFICIENTS TO USE FOR N=16
      DATA CFFT(1) / 9.939322535568174e-07 /
      DATA CFFT(2) / 3.981287575088865e-06 /
      DATA CFFT(3) / -5.584233411098927e-06 /
      DATA CFFT(4) / -2.734640462448423e-05 /
      DATA CFFT(5) / 2.170986793223473e-05 /
      DATA CFFT(6) / 2.107105639653287e-04 /
      DATA CFFT(7) / 8.703158428458035e-05 /
      DATA CFFT(8) / -0.001527659740122 /
      DATA CFFT(9) / -0.003881015189023 /
      DATA CFFT(10) / 0.003682567317092 /
      DATA CFFT(11) / 0.051822402431612 /
      DATA CFFT(12) / 0.191241726746695 /
      DATA CFFT(13) / 0.469290900903604 /
      DATA CFFT(14) / 0.886447830205055 /
      DATA CFFT(15) / 1.362240822271959 /
      DATA CFFT(16) / 1.748395886081962 /
      FLAG = .FALSE.
      i = DCMPLX(0,1)
      N = 16
      L = 2**(-0.25)*DSQRT(N*1.0d0)

      ZR1=Z1
      ZR2=Z2
      YI1=DIMAG(ZR1)
      YI2=DIMAG(ZR2)
      IF (YI1.LT.0) ZR1=-ZR1
      IF (YI2.LT.0) ZR2=-ZR2

      ZZ1 = (L+i*ZR1)/(L-i*ZR1)
      ZZ2 = (L+i*ZR2)/(L-i*ZR2)

      A1 = (1/DSQRT(PI))/(L-i*ZR1)
      B1 = 2/((L-i*ZR1)**2)
      A2 = (1/DSQRT(PI))/(L-i*ZR2)
      B2 = 2/((L-i*ZR2)**2)

      CALL CPOLYEV2(N,ZZ1,ZZ2,CFFT,P1,P2)

      ZZ1 = A1+B1*P1
      ZZ2 = A2+B2*P2
      IF (YI1.LT.0) ZZ1=-ZZ1
      IF (YI2.LT.0) ZZ2=-ZZ2
      RES=ZZ1*Z1**M-ZZ2*Z2**M
      RETURN
  100 FLAG = .TRUE.
      RETURN
*
      END
