1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
|
C ALGORITHM 680, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 16, NO. 1, PP. 47.
SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
C
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).
C THE ACCURACY OF THE ALGORITHM FOR Z IN THE 1ST AND 2ND QUADRANT
C IS 14 SIGNIFICANT DIGITS; IN THE 3RD AND 4TH IT IS 13 SIGNIFICANT
C DIGITS OUTSIDE A CIRCULAR REGION WITH RADIUS 0.126 AROUND A ZERO
C OF THE FUNCTION.
C ALL REAL VARIABLES IN THE PROGRAM ARE DOUBLE PRECISION.
C
C
C THE CODE CONTAINS A FEW COMPILER-DEPENDENT PARAMETERS :
C RMAXREAL = THE MAXIMUM VALUE OF RMAXREAL EQUALS THE ROOT OF
C RMAX = THE LARGEST NUMBER WHICH CAN STILL BE
C IMPLEMENTED ON THE COMPUTER IN DOUBLE PRECISION
C FLOATING-POINT ARITHMETIC
C RMAXEXP = LN(RMAX) - LN(2)
C RMAXGONI = THE LARGEST POSSIBLE ARGUMENT OF A DOUBLE PRECISION
C GONIOMETRIC FUNCTION (DCOS, DSIN, ...)
C THE REASON WHY THESE PARAMETERS ARE NEEDED AS THEY ARE DEFINED WILL
C BE EXPLAINED IN THE CODE BY MEANS OF COMMENTS
C
C
C PARAMETER LIST
C XI = REAL PART OF Z
C YI = IMAGINARY PART OF Z
C U = REAL PART OF W(Z)
C V = IMAGINARY PART OF W(Z)
C FLAG = AN ERROR FLAG INDICATING WHETHER OVERFLOW WILL
C OCCUR OR NOT; TYPE LOGICAL;
C THE VALUES OF THIS VARIABLE HAVE THE FOLLOWING
C MEANING :
C FLAG=.FALSE. : NO ERROR CONDITION
C FLAG=.TRUE. : OVERFLOW WILL OCCUR, THE ROUTINE
C BECOMES INACTIVE
C XI, YI ARE THE INPUT-PARAMETERS
C U, V, FLAG ARE THE OUTPUT-PARAMETERS
C
C FURTHERMORE THE PARAMETER FACTOR EQUALS 2/SQRT(PI)
C
C THE ROUTINE IS NOT UNDERFLOW-PROTECTED BUT ANY VARIABLE CAN BE
C PUT TO 0 UPON UNDERFLOW;
C
C REFERENCE - GPM POPPE, CMJ WIJERS; MORE EFFICIENT COMPUTATION OF
C THE COMPLEX ERROR-FUNCTION, ACM TRANS. MATH. SOFTWARE.
C
*
*
*
*
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
*
LOGICAL A, B, FLAG
PARAMETER (FACTOR = 1.12837916709551257388D0,
* RMAXREAL = 0.5D+154,
* RMAXEXP = 708.503061461606D0,
* RMAXGONI = 3.53711887601422D+15)
*
FLAG = .FALSE.
*
XABS = DABS(XI)
YABS = DABS(YI)
X = XABS/6.3
Y = YABS/4.4
*
C
C THE FOLLOWING IF-STATEMENT PROTECTS
C QRHO = (X**2 + Y**2) AGAINST OVERFLOW
C
IF ((XABS.GT.RMAXREAL).OR.(YABS.GT.RMAXREAL)) GOTO 100
*
QRHO = X**2 + Y**2
*
XABSQ = XABS**2
XQUAD = XABSQ - YABS**2
YQUAD = 2*XABS*YABS
*
A = QRHO.LT.0.085264D0
*
IF (A) THEN
C
C IF (QRHO.LT.0.085264D0) THEN THE FADDEEVA-FUNCTION IS EVALUATED
C USING A POWER-SERIES (ABRAMOWITZ/STEGUN, EQUATION (7.1.5), P.297)
C N IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED
C ACCURACY
C
QRHO = (1-0.85*Y)*DSQRT(QRHO)
N = IDNINT(6 + 72*QRHO)
J = 2*N+1
XSUM = 1.0/J
YSUM = 0.0D0
DO 10 I=N, 1, -1
J = J - 2
XAUX = (XSUM*XQUAD - YSUM*YQUAD)/I
YSUM = (XSUM*YQUAD + YSUM*XQUAD)/I
XSUM = XAUX + 1.0/J
10 CONTINUE
U1 = -FACTOR*(XSUM*YABS + YSUM*XABS) + 1.0
V1 = FACTOR*(XSUM*XABS - YSUM*YABS)
DAUX = DEXP(-XQUAD)
U2 = DAUX*DCOS(YQUAD)
V2 = -DAUX*DSIN(YQUAD)
*
U = U1*U2 - V1*V2
V = U1*V2 + V1*U2
*
ELSE
C
C IF (QRHO.GT.1.O) THEN W(Z) IS EVALUATED USING THE LAPLACE
C CONTINUED FRACTION
C NU IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED
C ACCURACY
C
C IF ((QRHO.GT.0.085264D0).AND.(QRHO.LT.1.0)) THEN W(Z) IS EVALUATED
C BY A TRUNCATED TAYLOR EXPANSION, WHERE THE LAPLACE CONTINUED FRACTION
C IS USED TO CALCULATE THE DERIVATIVES OF W(Z)
C KAPN IS THE MINIMUM NUMBER OF TERMS IN THE TAYLOR EXPANSION NEEDED
C TO OBTAIN THE REQUIRED ACCURACY
C NU IS THE MINIMUM NUMBER OF TERMS OF THE CONTINUED FRACTION NEEDED
C TO CALCULATE THE DERIVATIVES WITH THE REQUIRED ACCURACY
C
*
IF (QRHO.GT.1.0) THEN
H = 0.0D0
KAPN = 0
QRHO = DSQRT(QRHO)
NU = IDINT(3 + (1442/(26*QRHO+77)))
ELSE
QRHO = (1-Y)*DSQRT(1-QRHO)
H = 1.88*QRHO
H2 = 2*H
KAPN = IDNINT(7 + 34*QRHO)
NU = IDNINT(16 + 26*QRHO)
ENDIF
*
B = (H.GT.0.0)
*
IF (B) QLAMBDA = H2**KAPN
*
RX = 0.0
RY = 0.0
SX = 0.0
SY = 0.0
*
DO 11 N=NU, 0, -1
NP1 = N + 1
TX = YABS + H + NP1*RX
TY = XABS - NP1*RY
C = 0.5/(TX**2 + TY**2)
RX = C*TX
RY = C*TY
IF ((B).AND.(N.LE.KAPN)) THEN
TX = QLAMBDA + SX
SX = RX*TX - RY*SY
SY = RY*TX + RX*SY
QLAMBDA = QLAMBDA/H2
ENDIF
11 CONTINUE
*
IF (H.EQ.0.0) THEN
U = FACTOR*RX
V = FACTOR*RY
ELSE
U = FACTOR*SX
V = FACTOR*SY
END IF
*
IF (YABS.EQ.0.0) U = DEXP(-XABS**2)
*
END IF
*
*
C
C EVALUATION OF W(Z) IN THE OTHER QUADRANTS
C
*
IF (YI.LT.0.0) THEN
*
IF (A) THEN
U2 = 2*U2
V2 = 2*V2
ELSE
XQUAD = -XQUAD
*
C
C THE FOLLOWING IF-STATEMENT PROTECTS 2*EXP(-Z**2)
C AGAINST OVERFLOW
C
IF ((YQUAD.GT.RMAXGONI).OR.
* (XQUAD.GT.RMAXEXP)) GOTO 100
*
W1 = 2*DEXP(XQUAD)
U2 = W1*DCOS(YQUAD)
V2 = -W1*DSIN(YQUAD)
END IF
*
U = U2 - U
V = V2 - V
IF (XI.GT.0.0) V = -V
ELSE
IF (XI.LT.0.0) V = -V
END IF
*
RETURN
*
100 FLAG = .TRUE.
RETURN
*
END
|
