1
SUBROUTINE ZAIRY(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
2
C***BEGIN PROLOGUE ZAIRY
3
C***DATE WRITTEN 830501 (YYMMDD)
4
C***REVISION DATE 890801 (YYMMDD)
6
C***KEYWORDS AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD
7
C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
8
C***PURPOSE TO COMPUTE AIRY FUNCTIONS AI(Z) AND DAI(Z) FOR COMPLEX Z
11
C ***A DOUBLE PRECISION ROUTINE***
12
C ON KODE=1, ZAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR
13
C ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON
14
C KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)*
15
C DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN
16
C -PI/3.LT.ARG(Z).LT.PI/3 AND THE EXPONENTIAL GROWTH IN
17
C PI/3.LT.ABS(ARG(Z)).LT.PI WHERE ZTA=(2/3)*Z*CSQRT(Z).
19
C WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTIC IN
20
C THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED
21
C FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS.
22
C DEFINTIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF
23
C MATHEMATICAL FUNCTIONS (REF. 1).
25
C INPUT ZR,ZI ARE DOUBLE PRECISION
26
C ZR,ZI - Z=CMPLX(ZR,ZI)
27
C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1
28
C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
31
C AI=DAI(Z)/DZ ON ID=1
33
C AI=CEXP(ZTA)*AI(Z) ON ID=0 OR
34
C AI=CEXP(ZTA)*DAI(Z)/DZ ON ID=1 WHERE
35
C ZTA=(2/3)*Z*CSQRT(Z)
37
C OUTPUT AIR,AII ARE DOUBLE PRECISION
38
C AIR,AII- COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND
40
C NZ - UNDERFLOW INDICATOR
41
C NZ= 0 , NORMAL RETURN
42
C NZ= 1 , AI=CMPLX(0.0D0,0.0D0) DUE TO UNDERFLOW IN
43
C -PI/3.LT.ARG(Z).LT.PI/3 ON KODE=1
45
C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
46
C IERR=1, INPUT ERROR - NO COMPUTATION
47
C IERR=2, OVERFLOW - NO COMPUTATION, REAL(ZTA)
49
C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED
50
C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION
51
C PRODUCE LESS THAN HALF OF MACHINE ACCURACY
52
C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION
53
C COMPLETE LOSS OF ACCURACY BY ARGUMENT
55
C IERR=5, ERROR - NO COMPUTATION,
56
C ALGORITHM TERMINATION CONDITION NOT MET
60
C AI AND DAI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE K BESSEL
63
C AI(Z)=C*SQRT(Z)*K(1/3,ZTA) , DAI(Z)=-C*Z*K(2/3,ZTA)
64
C C=1.0/(PI*SQRT(3.0))
67
C WITH THE POWER SERIES FOR CABS(Z).LE.1.0.
69
C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
70
C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES
71
C OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF
72
C THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR),
73
C THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR
74
C FLAG IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
75
C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
76
C ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN
77
C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT
78
C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE
79
C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA
80
C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2,
81
C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE
82
C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE
83
C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT-
84
C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG-
85
C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN
86
C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN
87
C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES,
88
C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE
89
C PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER
92
C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
93
C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
94
C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
95
C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
96
C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
97
C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
98
C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
99
C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
100
C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
101
C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
102
C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
103
C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
104
C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
105
C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
106
C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
107
C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
108
C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
109
C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
112
C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
113
C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
116
C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
117
C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
119
C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
120
C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
123
C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
124
C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
125
C MATH. SOFTWARE, 1986
127
C***ROUTINES CALLED ZACAI,ZBKNU,AZEXP,AZSQRT,I1MACH,D1MACH
128
C***END PROLOGUE ZAIRY
129
C COMPLEX AI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3
130
DOUBLE PRECISION AA, AD, AII, AIR, AK, ALIM, ATRM, AZ, AZ3, BK,
131
* CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2, DIG,
132
* DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR,
133
* S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI,
134
* ZEROR, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, AZABS, ALAZ, BB
135
INTEGER ID, IERR, IFLAG, K, KODE, K1, K2, MR, NN, NZ, I1MACH
136
DIMENSION CYR(1), CYI(1)
137
DATA TTH, C1, C2, COEF /6.66666666666666667D-01,
138
* 3.55028053887817240D-01,2.58819403792806799D-01,
139
* 1.83776298473930683D-01/
140
DATA ZEROR, ZEROI, CONER, CONEI /0.0D0,0.0D0,1.0D0,0.0D0/
141
C***FIRST EXECUTABLE STATEMENT ZAIRY
144
IF (ID.LT.0 .OR. ID.GT.1) IERR=1
145
IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
146
IF (IERR.NE.0) RETURN
148
TOL = DMAX1(D1MACH(4),1.0D-18)
149
FID = DBLE(FLOAT(ID))
150
IF (AZ.GT.1.0D0) GO TO 70
151
C-----------------------------------------------------------------------
152
C POWER SERIES FOR CABS(Z).LE.1.
153
C-----------------------------------------------------------------------
158
IF (AZ.LT.TOL) GO TO 170
160
IF (AA.LT.TOL/AZ) GO TO 40
168
Z3R = STR*ZR - STI*ZI
169
Z3I = STR*ZI + STI*ZR
172
BK = 3.0D0 - FID - FID
174
DK = 3.0D0 + FID + FID
178
AK = 24.0D0 + 9.0D0*FID
179
BK = 30.0D0 - 9.0D0*FID
181
STR = (TRM1R*Z3R-TRM1I*Z3I)/D1
182
TRM1I = (TRM1R*Z3I+TRM1I*Z3R)/D1
186
STR = (TRM2R*Z3R-TRM2I*Z3I)/D2
187
TRM2I = (TRM2R*Z3I+TRM2I*Z3R)/D2
195
IF (ATRM.LT.TOL*AD) GO TO 40
200
IF (ID.EQ.1) GO TO 50
201
AIR = S1R*C1 - C2*(ZR*S2R-ZI*S2I)
202
AII = S1I*C1 - C2*(ZR*S2I+ZI*S2R)
203
IF (KODE.EQ.1) RETURN
204
CALL AZSQRT(ZR, ZI, STR, STI)
205
ZTAR = TTH*(ZR*STR-ZI*STI)
206
ZTAI = TTH*(ZR*STI+ZI*STR)
207
CALL AZEXP(ZTAR, ZTAI, STR, STI)
208
PTR = AIR*STR - AII*STI
209
AII = AIR*STI + AII*STR
215
IF (AZ.LE.TOL) GO TO 60
216
STR = ZR*S1R - ZI*S1I
217
STI = ZR*S1I + ZI*S1R
219
AIR = AIR + CC*(STR*ZR-STI*ZI)
220
AII = AII + CC*(STR*ZI+STI*ZR)
222
IF (KODE.EQ.1) RETURN
223
CALL AZSQRT(ZR, ZI, STR, STI)
224
ZTAR = TTH*(ZR*STR-ZI*STI)
225
ZTAI = TTH*(ZR*STI+ZI*STR)
226
CALL AZEXP(ZTAR, ZTAI, STR, STI)
227
PTR = STR*AIR - STI*AII
228
AII = STR*AII + STI*AIR
231
C-----------------------------------------------------------------------
232
C CASE FOR CABS(Z).GT.1.0
233
C-----------------------------------------------------------------------
235
FNU = (1.0D0+FID)/3.0D0
236
C-----------------------------------------------------------------------
237
C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
238
C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0D-18.
239
C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
240
C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
241
C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
242
C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
243
C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
244
C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
245
C-----------------------------------------------------------------------
249
K = MIN0(IABS(K1),IABS(K2))
250
ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
252
AA = R1M5*DBLE(FLOAT(K1))
253
DIG = DMIN1(AA,18.0D0)
255
ALIM = ELIM + DMAX1(-AA,-41.45D0)
256
RL = 1.2D0*DIG + 3.0D0
258
C--------------------------------------------------------------------------
259
C TEST FOR PROPER RANGE
260
C-----------------------------------------------------------------------
262
BB=DBLE(FLOAT(I1MACH(9)))*0.5D0
265
IF (AZ.GT.AA) GO TO 260
268
CALL AZSQRT(ZR, ZI, CSQR, CSQI)
269
ZTAR = TTH*(ZR*CSQR-ZI*CSQI)
270
ZTAI = TTH*(ZR*CSQI+ZI*CSQR)
271
C-----------------------------------------------------------------------
272
C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
273
C-----------------------------------------------------------------------
277
IF (ZR.GE.0.0D0) GO TO 80
283
IF (ZI.NE.0.0D0) GO TO 90
284
IF (ZR.GT.0.0D0) GO TO 90
289
IF (AA.GE.0.0D0 .AND. ZR.GT.0.0D0) GO TO 110
290
IF (KODE.EQ.2) GO TO 100
291
C-----------------------------------------------------------------------
293
C-----------------------------------------------------------------------
294
IF (AA.GT.(-ALIM)) GO TO 100
295
AA = -AA + 0.25D0*ALAZ
298
IF (AA.GT.ELIM) GO TO 270
300
C-----------------------------------------------------------------------
301
C CBKNU AND CACON RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
302
C-----------------------------------------------------------------------
304
IF (ZI.LT.0.0D0) MR = -1
305
CALL ZACAI(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, NN, RL, TOL,
307
IF (NN.LT.0) GO TO 280
311
IF (KODE.EQ.2) GO TO 120
312
C-----------------------------------------------------------------------
314
C-----------------------------------------------------------------------
315
IF (AA.LT.ALIM) GO TO 120
316
AA = -AA - 0.25D0*ALAZ
319
IF (AA.LT.(-ELIM)) GO TO 210
321
CALL ZBKNU(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM,
326
IF (IFLAG.NE.0) GO TO 150
327
IF (ID.EQ.1) GO TO 140
328
AIR = CSQR*S1R - CSQI*S1I
329
AII = CSQR*S1I + CSQI*S1R
332
AIR = -(ZR*S1R-ZI*S1I)
333
AII = -(ZR*S1I+ZI*S1R)
338
IF (ID.EQ.1) GO TO 160
339
STR = S1R*CSQR - S1I*CSQI
340
S1I = S1R*CSQI + S1I*CSQR
346
STR = -(S1R*ZR-S1I*ZI)
347
S1I = -(S1R*ZI+S1I*ZR)
353
AA = 1.0D+3*D1MACH(1)
356
IF (ID.EQ.1) GO TO 190
357
IF (AZ.LE.AA) GO TO 180
368
IF (AZ.LE.AA) GO TO 200
369
S1R = 0.5D0*(ZR*ZR-ZI*ZI)
385
IF(NN.EQ.(-1)) GO TO 270
added
removed
Lib/special/amos/zairy.f
