Description of erfcc

0001 C      ALGORITHM 680, COLLECTED ALGORITHMS FROM ACM.
0002 C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
0003 C      VOL. 16, NO. 1, PP. 47.
0004       SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
0005 C
0006 C  GIVEN A COMPLEX NUMBER Z = (XI,YI), THIS SUBROUTINE COMPUTES
0007 C  THE VALUE OF THE FADDEEVA-FUNCTION W(Z) = EXP(-Z**2)*ERFC(-I*Z),
0008 C  WHERE ERFC IS THE COMPLEX COMPLEMENTARY ERROR-FUNCTION AND I
0009 C  MEANS SQRT(-1).
0010 C  THE ACCURACY OF THE ALGORITHM FOR Z IN THE 1ST AND 2ND QUADRANT
0011 C  IS 14 SIGNIFICANT DIGITS; IN THE 3RD AND 4TH IT IS 13 SIGNIFICANT
0012 C  DIGITS OUTSIDE A CIRCULAR REGION WITH RADIUS 0.126 AROUND A ZERO
0013 C  OF THE FUNCTION.
0014 C  ALL REAL VARIABLES IN THE PROGRAM ARE DOUBLE PRECISION.
0015 C
0016 C
0017 C  THE CODE CONTAINS A FEW COMPILER-DEPENDENT PARAMETERS :
0018 C     RMAXREAL = THE MAXIMUM VALUE OF RMAXREAL EQUALS THE ROOT OF
0019 C                RMAX = THE LARGEST NUMBER WHICH CAN STILL BE
0020 C                IMPLEMENTED ON THE COMPUTER IN DOUBLE PRECISION
0021 C                FLOATING-POINT ARITHMETIC
0022 C     RMAXEXP  = LN(RMAX) - LN(2)
0023 C     RMAXGONI = THE LARGEST POSSIBLE ARGUMENT OF A DOUBLE PRECISION
0024 C                GONIOMETRIC FUNCTION (DCOS, DSIN, ...)
0025 C  THE REASON WHY THESE PARAMETERS ARE NEEDED AS THEY ARE DEFINED WILL
0026 C  BE EXPLAINED IN THE CODE BY MEANS OF COMMENTS
0027 C
0028 C
0029 C  PARAMETER LIST
0030 C     XI     = REAL      PART OF Z
0031 C     YI     = IMAGINARY PART OF Z
0032 C     U      = REAL      PART OF W(Z)
0033 C     V      = IMAGINARY PART OF W(Z)
0034 C     FLAG   = AN ERROR FLAG INDICATING WHETHER OVERFLOW WILL
0035 C              OCCUR OR NOT; TYPE LOGICAL;
0036 C              THE VALUES OF THIS VARIABLE HAVE THE FOLLOWING
0037 C              MEANING :
0038 C              FLAG=.FALSE. : NO ERROR CONDITION
0039 C              FLAG=.TRUE.  : OVERFLOW WILL OCCUR, THE ROUTINE
0040 C                             BECOMES INACTIVE
0041 C  XI, YI      ARE THE INPUT-PARAMETERS
0042 C  U, V, FLAG  ARE THE OUTPUT-PARAMETERS
0043 C
0044 C  FURTHERMORE THE PARAMETER FACTOR EQUALS 2/SQRT(PI)
0045 C
0046 C  THE ROUTINE IS NOT UNDERFLOW-PROTECTED BUT ANY VARIABLE CAN BE
0047 C  PUT TO 0 UPON UNDERFLOW;
0048 C
0049 C  REFERENCE - GPM POPPE, CMJ WIJERS; MORE EFFICIENT COMPUTATION OF
0050 C  THE COMPLEX ERROR-FUNCTION, ACM TRANS. MATH. SOFTWARE.
0051 C
0052 *
0053 *
0054 *
0055 *
0056       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
0057 *
0058       LOGICAL A, B, FLAG
0059       PARAMETER (FACTOR   = 1.12837916709551257388D0,
0060      *           RMAXREAL = 0.5D+154,
0061      *           RMAXEXP  = 708.503061461606D0,
0062      *           RMAXGONI = 3.53711887601422D+15)
0063 *
0064       FLAG = .FALSE.
0065 *
0066       XABS = DABS(XI)
0067       YABS = DABS(YI)
0068       X    = XABS/6.3
0069       Y    = YABS/4.4
0070 *
0071 C
0072 C     THE FOLLOWING IF-STATEMENT PROTECTS
0073 C     QRHO = (X**2 + Y**2) AGAINST OVERFLOW
0074 C
0075       IF ((XABS.GT.RMAXREAL).OR.(YABS.GT.RMAXREAL)) GOTO 100
0076 *
0077       QRHO = X**2 + Y**2
0078 *
0079       XABSQ = XABS**2
0080       XQUAD = XABSQ - YABS**2
0081       YQUAD = 2*XABS*YABS
0082 *
0083       A     = QRHO.LT.0.085264D0
0084 *
0085       IF (A) THEN
0086 C
0087 C  IF (QRHO.LT.0.085264D0) THEN THE FADDEEVA-FUNCTION IS EVALUATED
0088 C  USING A POWER-SERIES (ABRAMOWITZ/STEGUN, EQUATION (7.1.5), P.297)
0089 C  N IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED
0090 C  ACCURACY
0091 C
0092         QRHO  = (1-0.85*Y)*DSQRT(QRHO)
0093         N     = IDNINT(6 + 72*QRHO)
0094         J     = 2*N+1
0095         XSUM  = 1.0D0/J
0096         YSUM  = 0.0D0
0097         DO 10 I=N, 1, -1
0098           J    = J - 2
0099           XAUX = (XSUM*XQUAD - YSUM*YQUAD)/I
0100           YSUM = (XSUM*YQUAD + YSUM*XQUAD)/I
0101           XSUM = XAUX + 1.0D0/J
0102  10     CONTINUE
0103         U1   = -FACTOR*(XSUM*YABS + YSUM*XABS) + 1.0D0
0104         V1   =  FACTOR*(XSUM*XABS - YSUM*YABS)
0105         DAUX =  DEXP(-XQUAD)
0106         U2   =  DAUX*DCOS(YQUAD)
0107         V2   = -DAUX*DSIN(YQUAD)
0108 *
0109         U    = U1*U2 - V1*V2
0110         V    = U1*V2 + V1*U2
0111 *
0112       ELSE
0113 C
0114 C  IF (QRHO.GT.1.O) THEN W(Z) IS EVALUATED USING THE LAPLACE
0115 C  CONTINUED FRACTION
0116 C  NU IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED
0117 C  ACCURACY
0118 C
0119 C  IF ((QRHO.GT.0.085264D0).AND.(QRHO.LT.1.0)) THEN W(Z) IS EVALUATED
0120 C  BY A TRUNCATED TAYLOR EXPANSION, WHERE THE LAPLACE CONTINUED FRACTION
0121 C  IS USED TO CALCULATE THE DERIVATIVES OF W(Z)
0122 C  KAPN IS THE MINIMUM NUMBER OF TERMS IN THE TAYLOR EXPANSION NEEDED
0123 C  TO OBTAIN THE REQUIRED ACCURACY
0124 C  NU IS THE MINIMUM NUMBER OF TERMS OF THE CONTINUED FRACTION NEEDED
0125 C  TO CALCULATE THE DERIVATIVES WITH THE REQUIRED ACCURACY
0126 C
0127 *
0128         IF (QRHO.GT.1.0) THEN
0129           H    = 0.0D0
0130           KAPN = 0
0131           QRHO = DSQRT(QRHO)
0132           NU   = IDINT(3 + (1442/(26*QRHO+77)))
0133         ELSE
0134           QRHO = (1-Y)*DSQRT(1-QRHO)
0135           H    = 1.88*QRHO
0136           H2   = 2*H
0137           KAPN = IDNINT(7  + 34*QRHO)
0138           NU   = IDNINT(16 + 26*QRHO)
0139         ENDIF
0140 *
0141         B = (H.GT.0.0)
0142 *
0143         IF (B) QLAMBDA = H2**KAPN
0144 *
0145         RX = 0.0
0146         RY = 0.0
0147         SX = 0.0
0148         SY = 0.0
0149 *
0150         DO 11 N=NU, 0, -1
0151           NP1 = N + 1
0152           TX  = YABS + H + NP1*RX
0153           TY  = XABS - NP1*RY
0154           C   = 0.5/(TX**2 + TY**2)
0155           RX  = C*TX
0156           RY  = C*TY
0157           IF ((B).AND.(N.LE.KAPN)) THEN
0158             TX = QLAMBDA + SX
0159             SX = RX*TX - RY*SY
0160             SY = RY*TX + RX*SY
0161             QLAMBDA = QLAMBDA/H2
0162           ENDIF
0163  11     CONTINUE
0164 *
0165         IF (H.EQ.0.0) THEN
0166           U = FACTOR*RX
0167           V = FACTOR*RY
0168         ELSE
0169           U = FACTOR*SX
0170           V = FACTOR*SY
0171         END IF
0172 *
0173         IF (YABS.EQ.0.0) U = DEXP(-XABS**2)
0174 *
0175       END IF
0176 *
0177 *
0178 C
0179 C  EVALUATION OF W(Z) IN THE OTHER QUADRANTS
0180 C
0181 *
0182       IF (YI.LT.0.0) THEN
0183 *
0184         IF (A) THEN
0185           U2    = 2*U2
0186           V2    = 2*V2
0187         ELSE
0188           XQUAD =  -XQUAD
0189 *
0190 C
0191 C         THE FOLLOWING IF-STATEMENT PROTECTS 2*EXP(-Z**2)
0192 C         AGAINST OVERFLOW
0193 C
0194           IF ((YQUAD.GT.RMAXGONI).OR.
0195      *        (XQUAD.GT.RMAXEXP)) GOTO 100
0196 *
0197           W1 =  2*DEXP(XQUAD)
0198           U2  =  W1*DCOS(YQUAD)
0199           V2  = -W1*DSIN(YQUAD)
0200         END IF
0201 *
0202         U = U2 - U
0203         V = V2 - V
0204         IF (XI.GT.0.0) V = -V
0205       ELSE
0206         IF (XI.LT.0.0) V = -V
0207       END IF
0208 *
0209       RETURN
0210 *
0211   100 FLAG = .TRUE.
0212       RETURN
0213 *
0214       END
erfcc

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
