1      SUBROUTINE ZAIRY(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
2C***BEGIN PROLOGUE  ZAIRY
3C***DATE WRITTEN   830501   (YYMMDD)
4C***REVISION DATE  890801   (YYMMDD)
5C***CATEGORY NO.  B5K
6C***KEYWORDS  AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD
7C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
8C***PURPOSE  TO COMPUTE AIRY FUNCTIONS AI(Z) AND DAI(Z) FOR COMPLEX Z
9C***DESCRIPTION
10C
11C                      ***A DOUBLE PRECISION ROUTINE***
12C         ON KODE=1, ZAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR
13C         ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON
14C         KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)*
15C         DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN
16C         -PI/3.LT.ARG(Z).LT.PI/3 AND THE EXPONENTIAL GROWTH IN
17C         PI/3.LT.ABS(ARG(Z)).LT.PI WHERE ZTA=(2/3)*Z*CSQRT(Z).
18C
19C         WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTIC IN
20C         THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED
21C         FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS.
22C         DEFINTIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF
23C         MATHEMATICAL FUNCTIONS (REF. 1).
24C
25C         INPUT      ZR,ZI ARE DOUBLE PRECISION
26C           ZR,ZI  - Z=CMPLX(ZR,ZI)
27C           ID     - ORDER OF DERIVATIVE, ID=0 OR ID=1
28C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
29C                    KODE= 1  RETURNS
30C                             AI=AI(Z)                ON ID=0 OR
31C                             AI=DAI(Z)/DZ            ON ID=1
32C                        = 2  RETURNS
33C                             AI=CEXP(ZTA)*AI(Z)       ON ID=0 OR
34C                             AI=CEXP(ZTA)*DAI(Z)/DZ   ON ID=1 WHERE
35C                             ZTA=(2/3)*Z*CSQRT(Z)
36C
37C         OUTPUT     AIR,AII ARE DOUBLE PRECISION
38C           AIR,AII- COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND
39C                    KODE
40C           NZ     - UNDERFLOW INDICATOR
41C                    NZ= 0   , NORMAL RETURN
42C                    NZ= 1   , AI=CMPLX(0.0D0,0.0D0) DUE TO UNDERFLOW IN
43C                              -PI/3.LT.ARG(Z).LT.PI/3 ON KODE=1
44C           IERR   - ERROR FLAG
45C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
46C                    IERR=1, INPUT ERROR   - NO COMPUTATION
47C                    IERR=2, OVERFLOW      - NO COMPUTATION, REAL(ZTA)
48C                            TOO LARGE ON KODE=1
49C                    IERR=3, CABS(Z) LARGE      - COMPUTATION COMPLETED
50C                            LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION
51C                            PRODUCE LESS THAN HALF OF MACHINE ACCURACY
52C                    IERR=4, CABS(Z) TOO LARGE  - NO COMPUTATION
53C                            COMPLETE LOSS OF ACCURACY BY ARGUMENT
54C                            REDUCTION
55C                    IERR=5, ERROR              - NO COMPUTATION,
56C                            ALGORITHM TERMINATION CONDITION NOT MET
57C
58C***LONG DESCRIPTION
59C
60C         AI AND DAI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE K BESSEL
61C         FUNCTIONS BY
62C
63C            AI(Z)=C*SQRT(Z)*K(1/3,ZTA) , DAI(Z)=-C*Z*K(2/3,ZTA)
64C                           C=1.0/(PI*SQRT(3.0))
65C                            ZTA=(2/3)*Z**(3/2)
66C
67C         WITH THE POWER SERIES FOR CABS(Z).LE.1.0.
68C
69C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
70C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES
71C         OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF
72C         THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR),
73C         THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR
74C         FLAG IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
75C         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
76C         ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN
77C         ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT
78C         FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE
79C         LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA
80C         MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2,
81C         AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE
82C         PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE
83C         PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT-
84C         ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG-
85C         NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN
86C         DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN
87C         EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES,
88C         NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE
89C         PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER
90C         MACHINES.
91C
92C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
93C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
94C         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
95C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
96C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
97C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
98C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
99C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
100C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
101C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
102C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
103C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
104C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
105C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
106C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
107C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
108C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
109C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
110C         OR -PI/2+P.
111C
112C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
113C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
114C                 COMMERCE, 1955.
115C
116C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
117C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
118C
119C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
120C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
121C                 1018, MAY, 1985
122C
123C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
124C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
125C                 MATH. SOFTWARE, 1986
126C
127C***ROUTINES CALLED  ZACAI,ZBKNU,ZEXP,ZSQRT,I1MACH,D1MACH
128C***END PROLOGUE  ZAIRY
129C     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, ZABS, 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/
141C***FIRST EXECUTABLE STATEMENT  ZAIRY
142      IERR = 0
143      NZ=0
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
147      AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0)))
148      TOL = DMAX1(D1MACH(4),1.0D-18)
149      FID = DBLE(FLOAT(ID))
150      IF (AZ.GT.1.0D0) GO TO 70
151C-----------------------------------------------------------------------
152C     POWER SERIES FOR CABS(Z).LE.1.
153C-----------------------------------------------------------------------
154      S1R = CONER
155      S1I = CONEI
156      S2R = CONER
157      S2I = CONEI
158      IF (AZ.LT.TOL) GO TO 170
159      AA = AZ*AZ
160      IF (AA.LT.TOL/AZ) GO TO 40
161      TRM1R = CONER
162      TRM1I = CONEI
163      TRM2R = CONER
164      TRM2I = CONEI
165      ATRM = 1.0D0
166      STR = ZR*ZR - ZI*ZI
167      STI = ZR*ZI + ZI*ZR
168      Z3R = STR*ZR - STI*ZI
169      Z3I = STR*ZI + STI*ZR
170      AZ3 = AZ*AA
171      AK = 2.0D0 + FID
172      BK = 3.0D0 - FID - FID
173      CK = 4.0D0 - FID
174      DK = 3.0D0 + FID + FID
175      D1 = AK*DK
176      D2 = BK*CK
177      AD = DMIN1(D1,D2)
178      AK = 24.0D0 + 9.0D0*FID
179      BK = 30.0D0 - 9.0D0*FID
180      DO 30 K=1,25
181        STR = (TRM1R*Z3R-TRM1I*Z3I)/D1
182        TRM1I = (TRM1R*Z3I+TRM1I*Z3R)/D1
183        TRM1R = STR
184        S1R = S1R + TRM1R
185        S1I = S1I + TRM1I
186        STR = (TRM2R*Z3R-TRM2I*Z3I)/D2
187        TRM2I = (TRM2R*Z3I+TRM2I*Z3R)/D2
188        TRM2R = STR
189        S2R = S2R + TRM2R
190        S2I = S2I + TRM2I
191        ATRM = ATRM*AZ3/AD
192        D1 = D1 + AK
193        D2 = D2 + BK
194        AD = DMIN1(D1,D2)
195        IF (ATRM.LT.TOL*AD) GO TO 40
196        AK = AK + 18.0D0
197        BK = BK + 18.0D0
198   30 CONTINUE
199   40 CONTINUE
200      IF (ID.EQ.1) THEN
201        GO TO 50
202      END IF
203      AIR = S1R*C1 - C2*(ZR*S2R-ZI*S2I)
204      AII = S1I*C1 - C2*(ZR*S2I+ZI*S2R)
205      IF (KODE.EQ.1) RETURN
206      CALL ZSQRT(ZR, ZI, STR, STI)
207      ZTAR = TTH*(ZR*STR-ZI*STI)
208      ZTAI = TTH*(ZR*STI+ZI*STR)
209      CALL ZEXP(ZTAR, ZTAI, STR, STI)
210      PTR = AIR*STR - AII*STI
211      AII = AIR*STI + AII*STR
212      AIR = PTR
213      RETURN
214   50 CONTINUE
215      AIR = -S2R*C2
216      AII = -S2I*C2
217      IF (AZ.LE.TOL) GO TO 60
218      STR = ZR*S1R - ZI*S1I
219      STI = ZR*S1I + ZI*S1R
220      CC = C1/(1.0D0+FID)
221      AIR = AIR + CC*(STR*ZR-STI*ZI)
222      AII = AII + CC*(STR*ZI+STI*ZR)
223   60 CONTINUE
224      IF (KODE.EQ.1) RETURN
225      CALL ZSQRT(ZR, ZI, STR, STI)
226      ZTAR = TTH*(ZR*STR-ZI*STI)
227      ZTAI = TTH*(ZR*STI+ZI*STR)
228      CALL ZEXP(ZTAR, ZTAI, STR, STI)
229      PTR = STR*AIR - STI*AII
230      AII = STR*AII + STI*AIR
231      AIR = PTR
232      RETURN
233C-----------------------------------------------------------------------
234C     CASE FOR CABS(Z).GT.1.0
235C-----------------------------------------------------------------------
236   70 CONTINUE
237      FNU = (1.0D0+FID)/3.0D0
238C-----------------------------------------------------------------------
239C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
240C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0D-18.
241C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
242C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
243C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
244C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
245C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
246C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
247C-----------------------------------------------------------------------
248      K1 = I1MACH(15)
249      K2 = I1MACH(16)
250      R1M5 = D1MACH(5)
251      K = MIN0(IABS(K1),IABS(K2))
252      ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
253      K1 = I1MACH(14) - 1
254      AA = R1M5*DBLE(FLOAT(K1))
255      DIG = DMIN1(AA,18.0D0)
256      AA = AA*2.303D0
257      ALIM = ELIM + DMAX1(-AA,-41.45D0)
258      RL = 1.2D0*DIG + 3.0D0
259      ALAZ = DLOG(AZ)
260C--------------------------------------------------------------------------
261C     TEST FOR PROPER RANGE
262C-----------------------------------------------------------------------
263      AA=0.5D0/TOL
264      BB=DBLE(FLOAT(I1MACH(9)))*0.5D0
265      AA=DMIN1(AA,BB)
266      AA=AA**TTH
267      IF (AZ.GT.AA) GO TO 260
268      AA=DSQRT(AA)
269      IF (AZ.GT.AA) IERR=3
270      CALL ZSQRT(ZR, ZI, CSQR, CSQI)
271      ZTAR = TTH*(ZR*CSQR-ZI*CSQI)
272      ZTAI = TTH*(ZR*CSQI+ZI*CSQR)
273C-----------------------------------------------------------------------
274C     RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
275C-----------------------------------------------------------------------
276      IFLAG = 0
277      SFAC = 1.0D0
278      AK = ZTAI
279      IF (ZR.GE.0.0D0) GO TO 80
280      BK = ZTAR
281      CK = -DABS(BK)
282      ZTAR = CK
283      ZTAI = AK
284   80 CONTINUE
285      IF (ZI.NE.0.0D0) GO TO 90
286      IF (ZR.GT.0.0D0) GO TO 90
287      ZTAR = 0.0D0
288      ZTAI = AK
289   90 CONTINUE
290      AA = ZTAR
291      IF (AA.GE.0.0D0 .AND. ZR.GT.0.0D0) GO TO 110
292      IF (KODE.EQ.2) GO TO 100
293C-----------------------------------------------------------------------
294C     OVERFLOW TEST
295C-----------------------------------------------------------------------
296      IF (AA.GT.(-ALIM)) GO TO 100
297      AA = -AA + 0.25D0*ALAZ
298      IFLAG = 1
299      SFAC = TOL
300      IF (AA.GT.ELIM) GO TO 270
301  100 CONTINUE
302C-----------------------------------------------------------------------
303C     CBKNU AND CACON RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
304C-----------------------------------------------------------------------
305      MR = 1
306      IF (ZI.LT.0.0D0) MR = -1
307      CALL ZACAI(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, NN, RL, TOL,
308     * ELIM, ALIM)
309      IF (NN.LT.0) GO TO 280
310      NZ = NZ + NN
311      GO TO 130
312  110 CONTINUE
313      IF (KODE.EQ.2) GO TO 120
314C-----------------------------------------------------------------------
315C     UNDERFLOW TEST
316C-----------------------------------------------------------------------
317      IF (AA.LT.ALIM) GO TO 120
318      AA = -AA - 0.25D0*ALAZ
319      IFLAG = 2
320      SFAC = 1.0D0/TOL
321      IF (AA.LT.(-ELIM)) GO TO 210
322  120 CONTINUE
323      CALL ZBKNU(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM,
324     * ALIM)
325  130 CONTINUE
326      S1R = CYR(1)*COEF
327      S1I = CYI(1)*COEF
328      IF (IFLAG.NE.0) GO TO 150
329      IF (ID.EQ.1) GO TO 140
330      AIR = CSQR*S1R - CSQI*S1I
331      AII = CSQR*S1I + CSQI*S1R
332      RETURN
333  140 CONTINUE
334      AIR = -(ZR*S1R-ZI*S1I)
335      AII = -(ZR*S1I+ZI*S1R)
336      RETURN
337  150 CONTINUE
338      S1R = S1R*SFAC
339      S1I = S1I*SFAC
340      IF (ID.EQ.1) GO TO 160
341      STR = S1R*CSQR - S1I*CSQI
342      S1I = S1R*CSQI + S1I*CSQR
343      S1R = STR
344      AIR = S1R/SFAC
345      AII = S1I/SFAC
346      RETURN
347  160 CONTINUE
348      STR = -(S1R*ZR-S1I*ZI)
349      S1I = -(S1R*ZI+S1I*ZR)
350      S1R = STR
351      AIR = S1R/SFAC
352      AII = S1I/SFAC
353      RETURN
354  170 CONTINUE
355      AA = 1.0D+3*D1MACH(1)
356      S1R = ZEROR
357      S1I = ZEROI
358      IF (ID.EQ.1) GO TO 190
359      IF (AZ.LE.AA) GO TO 180
360      S1R = C2*ZR
361      S1I = C2*ZI
362  180 CONTINUE
363      AIR = C1 - S1R
364      AII = -S1I
365      RETURN
366  190 CONTINUE
367      AIR = -C2
368      AII = 0.0D0
369      AA = DSQRT(AA)
370      IF (AZ.LE.AA) GO TO 200
371      S1R = 0.5D0*(ZR*ZR-ZI*ZI)
372      S1I = ZR*ZI
373  200 CONTINUE
374      AIR = AIR + C1*S1R
375      AII = AII + C1*S1I
376      RETURN
377  210 CONTINUE
378      NZ = 1
379      AIR = ZEROR
380      AII = ZEROI
381      RETURN
382  270 CONTINUE
383      NZ = 0
384      IERR=2
385      RETURN
386  280 CONTINUE
387      IF(NN.EQ.(-1)) GO TO 270
388      NZ=0
389      IERR=5
390      RETURN
391  260 CONTINUE
392      IERR=4
393      NZ=0
394      RETURN
395      END
396