C            THE PROBABILITY INTEGRAL FOR COMPLEX ARGUMENT.             W(Z)  01
C                (ERROR FUNCTION FOR COMPLEX ARGUMENT.)                 W(Z)  02
                         COMPLEX FUNCTION W(Z)                          W(Z)  03
C                    W(Z) = EXP(-Z*Z) * ERFC(-I*Z)                      W(Z)  04
C      = EXP(-Z**2)*(1.+2.*J/ROOT PI*(INTEGRAL, 0 TO Z, (EXP(T**2).DT)))W(Z)  05
C      = J/PI*(INTEGRAL, -INFINITY TO +INFINITY, (EXP(-T**2)/(Z-T).DT)) W(Z)  06
C                                                                       W(Z)  07
C         THIS FUNCTION HAS BEEN TABULATED TO SIX DECIMALS BY           W(Z)  08
C                       FADDEYEVA  AND  TERENT@EV.                      W(Z)  09
C     ( USSR ACADEMY OF SCIENCES.  (1961)  PERGAMON.  ED. BY V.A. FOX.) W(Z)  10
C     ALSO        U.S. NATIONAL BUREAU OF STANDARDS.       (1964)       W(Z)  11
C                @HANDBOOK OF MATHEMATICAL FUNCTIONS@                   W(Z)  12
C                                                                       W(Z)  13
C***********************************************************************W(Z)  14
      COMPLEX Z,MZSQU                                                   W(Z)  15
      DOUBLE PRECISION A,B,C,PART1,PART2,PART3,PART4,PART1S,PART3S      W(Z)  16
      DATA C, D / 1.1283791670955126D0, 4.35             /              W(Z)
      DATA CLOWER,C25,CUPPER/"145440672724,"232400000000,"221652211200/ W(Z)  18
      ABS Z = CABS(Z)                                                   W(Z)  19
      IF(ABS Z.LE.C LOWER) GO TO 10                                     W(Z)  20
      IF(ABS Z.GE.C UPPER) GO TO 30                                     W(Z)  21
      MZSQU = -Z*Z                                                      W(Z)  22
      E =  ( ( CABS(Z+D) + CABS(Z-D) )/2. )**2   - D*D                  W(Z)  23
      IF(E.LT.1.) GO TO 20                                              W(Z)  24
C***********************************************************************W(Z)  25
C     CONTINUED FRACTION METHOD.                                        W(Z)  26
      M = 4. + 40./E                                                    W(Z)  27
      W = FLOAT(M)                                                      W(Z)  28
      DO 5 N = 1,M                                                      W(Z)  29
      RN = M - N + 1                                                    W(Z)  30
    5 W = RN*(RN-.5)/(MZSQU + 2.*RN + .5 - W)                           W(Z)  31
      W = (0.,-.564189584)*Z/(.5+MZSQU-W)                               W(Z)  32
      IF(REAL(MZSQU).LE.(-89.4159863)) RETURN                           W(Z)  35
      IF(AIMAG(Z).EQ.0.) W = CMPLX(EXP(REAL(MZSQU)),AIMAG(W))           W(Z)  33
      IF(AIMAG(Z).GE.0.) RETURN                                         W(Z)  34
      IF(ABS(AIMAG(MZSQU)).GE.C25) GO TO 201                            W(Z)  36
      IF(REAL(MZSQU).LT.(+87.336544 )) GO TO 8                          W(Z)  37
      WRITE(6,9200) Z                                                   W(Z)  38
      MZSQU = CMPLX(87.336544 ,AIMAG(MZSQU))                            W(Z)  39
    8 W = W + 2.*CEXP(MZSQU)                                            W(Z)  40
      RETURN                                                            W(Z)  41
C     RELATIVE INACCURACY AT COMPLEX ZEROS OF W IS IGNORED.             W(Z)  42
C     END OF CONTINUED FRACTION METHOD.                                 W(Z)  43
C***********************************************************************W(Z)  44
   10 W = CMPLX(1.,REAL(Z)*SNGL(C))                                     W(Z)  45
      RETURN                                                            W(Z)  46
C***********************************************************************W(Z)  47
C     DOUBLE-PRECISION CONVERGENT SERIES METHOD.                        W(Z)  48
   20 PART 1 = 1.                                                       W(Z)  49
      PART 2 = 0.                                                       W(Z)  50
      PART 3 = 1.                                                       W(Z)  51
      PART 4 = 0.                                                       W(Z)  52
      M = 4. + ABS Z * (7.5 + 2.0 * ABS Z)                              W(Z)  53
      A = REAL(MZSQU)                                                   W(Z)  54
      B = AIMAG(MZSQU)                                                  W(Z)  55
      DO 25 N = 1,M                                                     W(Z)  56
      PART 1S = PART 1                                                  W(Z)  57
      PART 3S = PART 3                                                  W(Z)  58
      RN = M - N + 1                                                    W(Z)  59
      RNPLUS = RN + 0.5                                                 W(Z)  60
      PART 1 = 1. + (A*PART 1 - B*PART 2)/RN                            W(Z)  61
      PART 2 = (B*PART 1S + A*PART 2)/RN                                W(Z)  62
      PART 3 = 1. + (A*PART 3 - B*PART 4) /RNPLUS                       W(Z)  63
   25 PART 4 = (A*PART 4 + B*PART 3S)/RNPLUS                            W(Z)  64
      PART 3 = PART 3 * C                                               W(Z)  65
      PART 4 = PART 4 * C                                               W(Z)  67
      PART 1 = PART 1 - PART 3 * AIMAG(Z) - PART 4 * REAL(Z)            W(Z)  68
      PART 2 = PART 2 + PART 3 * REAL(Z) - PART 4 * AIMAG(Z)            W(Z)  69
      W = CMPLX(SNGL(PART 1), SNGL(PART 2))                             W(Z)  70
      RETURN                                                            W(Z)  71
C     END OF CONVERGENT SERIES.                                         W(Z)  72
C***********************************************************************W(Z)  73
   30 W = (0.,.564189584)/Z                                             W(Z)  74
      IF(ABS(REAL(Z)).GT.(-AIMAG(Z))) RETURN                            W(Z)  75
  201 TYPE 9201,Z                                                       W(Z)  76
      RETURN                                                            W(Z)  77
 9200 FORMAT(31H SINCE W(Z) OVERFLOWS FOR Z = (,1PE17.9,1H,,E17.9,2H),,/W(Z)  78
     1,34H), AN APPROXIMATION HAS BEEN MADE.)                           W(Z)  79
 9201 FORMAT(38H SINCE W(Z) IS INDETERMINATE FOR Z = (,1PE17.9,1H,,E17.9W(Z)  80
     1,2H),,/,35H -W(-Z) HAS BEEN INSERTED FOR W(Z).)                   W(Z)  81
      END                                                               W(Z)  82

