1*DECK ZAIRY
2      SUBROUTINE ZAIRY (ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
3C***BEGIN PROLOGUE  ZAIRY
4C***PURPOSE  Compute the Airy function Ai(z) or its derivative dAi/dz
5C            for complex argument z.  A scaling option is available
6C            to help avoid underflow and overflow.
7C***LIBRARY   SLATEC
8C***CATEGORY  C10D
9C***TYPE      COMPLEX (CAIRY-C, ZAIRY-C)
10C***KEYWORDS  AIRY FUNCTION, BESSEL FUNCTION OF ORDER ONE THIRD,
11C             BESSEL FUNCTION OF ORDER TWO THIRDS
12C***AUTHOR  Amos, D. E., (SNL)
13C***DESCRIPTION
14C
15C                      ***A DOUBLE PRECISION ROUTINE***
16C         On KODE=1, ZAIRY computes the complex Airy function Ai(z)
17C         or its derivative dAi/dz on ID=0 or ID=1 respectively. On
18C         KODE=2, a scaling option exp(zeta)*Ai(z) or exp(zeta)*dAi/dz
19C         is provided to remove the exponential decay in -pi/3<arg(z)
20C         <pi/3 and the exponential growth in pi/3<abs(arg(z))<pi where
21C         zeta=(2/3)*z**(3/2).
22C
23C         While the Airy functions Ai(z) and dAi/dz are analytic in
24C         the whole z-plane, the corresponding scaled functions defined
25C         for KODE=2 have a cut along the negative real axis.
26C
27C         Input
28C           ZR     - DOUBLE PRECISION real part of argument Z
29C           ZI     - DOUBLE PRECISION imag part of argument Z
30C           ID     - Order of derivative, ID=0 or ID=1
31C           KODE   - A parameter to indicate the scaling option
32C                    KODE=1  returns
33C                            AI=Ai(z)  on ID=0
34C                            AI=dAi/dz on ID=1
35C                            at z=Z
36C                        =2  returns
37C                            AI=exp(zeta)*Ai(z)  on ID=0
38C                            AI=exp(zeta)*dAi/dz on ID=1
39C                            at z=Z where zeta=(2/3)*z**(3/2)
40C
41C         Output
42C           AIR    - DOUBLE PRECISION real part of result
43C           AII    - DOUBLE PRECISION imag part of result
44C           NZ     - Underflow indicator
45C                    NZ=0    Normal return
46C                    NZ=1    AI=0 due to underflow in
47C                            -pi/3<arg(Z)<pi/3 on KODE=1
48C           IERR   - Error flag
49C                    IERR=0  Normal return     - COMPUTATION COMPLETED
50C                    IERR=1  Input error       - NO COMPUTATION
51C                    IERR=2  Overflow          - NO COMPUTATION
52C                            (Re(Z) too large with KODE=1)
53C                    IERR=3  Precision warning - COMPUTATION COMPLETED
54C                            (Result has less than half precision)
55C                    IERR=4  Precision error   - NO COMPUTATION
56C                            (Result has no precision)
57C                    IERR=5  Algorithmic error - NO COMPUTATION
58C                            (Termination condition not met)
59C
60C *Long Description:
61C
62C         Ai(z) and dAi/dz are computed from K Bessel functions by
63C
64C                Ai(z) =  c*sqrt(z)*K(1/3,zeta)
65C               dAi/dz = -c*   z   *K(2/3,zeta)
66C                    c =  1/(pi*sqrt(3))
67C                 zeta =  (2/3)*z**(3/2)
68C
69C         when abs(z)>1 and from power series when abs(z)<=1.
70C
71C         In most complex variable computation, one must evaluate ele-
72C         mentary functions.  When the magnitude of Z is large, losses
73C         of significance by argument reduction occur.  Consequently, if
74C         the magnitude of ZETA=(2/3)*Z**(3/2) exceeds U1=SQRT(0.5/UR),
75C         then losses exceeding half precision are likely and an error
76C         flag IERR=3 is triggered where UR=MAX(D1MACH(4),1.0D-18) is
77C         double precision unit roundoff limited to 18 digits precision.
78C         Also, if the magnitude of ZETA is larger than U2=0.5/UR, then
79C         all significance is lost and IERR=4.  In order to use the INT
80C         function, ZETA must be further restricted not to exceed
81C         U3=I1MACH(9)=LARGEST INTEGER.  Thus, the magnitude of ZETA
82C         must be restricted by MIN(U2,U3).  In IEEE arithmetic, U1,U2,
83C         and U3 are approximately 2.0E+3, 4.2E+6, 2.1E+9 in single
84C         precision and 4.7E+7, 2.3E+15, 2.1E+9 in double precision.
85C         This makes U2 limiting is single precision and U3 limiting
86C         in double precision.  This means that the magnitude of Z
87C         cannot exceed approximately 3.4E+4 in single precision and
88C         2.1E+6 in double precision.  This also means that one can
89C         expect to retain, in the worst cases on 32-bit machines,
90C         no digits in single precision and only 6 digits in double
91C         precision.
92C
93C         The approximate relative error in the magnitude of a complex
94C         Bessel function can be expressed as P*10**S where P=MAX(UNIT
95C         ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre-
96C         sents the increase in error due to argument reduction in the
97C         elementary functions.  Here, S=MAX(1,ABS(LOG10(ABS(Z))),
98C         ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF
99C         ABS(Z),ABS(EXPONENT OF FNU)) ).  However, the phase angle may
100C         have only absolute accuracy.  This is most likely to occur
101C         when one component (in magnitude) is larger than the other by
102C         several orders of magnitude.  If one component is 10**K larger
103C         than the other, then one can expect only MAX(ABS(LOG10(P))-K,
104C         0) significant digits; or, stated another way, when K exceeds
105C         the exponent of P, no significant digits remain in the smaller
106C         component.  However, the phase angle retains absolute accuracy
107C         because, in complex arithmetic with precision P, the smaller
108C         component will not (as a rule) decrease below P times the
109C         magnitude of the larger component. In these extreme cases,
110C         the principal phase angle is on the order of +P, -P, PI/2-P,
111C         or -PI/2+P.
112C
113C***REFERENCES  1. M. Abramowitz and I. A. Stegun, Handbook of Mathe-
114C                 matical Functions, National Bureau of Standards
115C                 Applied Mathematics Series 55, U. S. Department
116C                 of Commerce, Tenth Printing (1972) or later.
117C               2. D. E. Amos, Computation of Bessel Functions of
118C                 Complex Argument and Large Order, Report SAND83-0643,
119C                 Sandia National Laboratories, Albuquerque, NM, May
120C                 1983.
121C               3. D. E. Amos, A Subroutine Package for Bessel Functions
122C                 of a Complex Argument and Nonnegative Order, Report
123C                 SAND85-1018, Sandia National Laboratory, Albuquerque,
124C                 NM, May 1985.
125C               4. D. E. Amos, A portable package for Bessel functions
126C                 of a complex argument and nonnegative order, ACM
127C                 Transactions on Mathematical Software, 12 (September
128C                 1986), pp. 265-273.
129C
130C***ROUTINES CALLED  D1MACH, I1MACH, ZABS, ZACAI, ZBKNU, ZEXP, ZSQRT
131C***REVISION HISTORY  (YYMMDD)
132C   830501  DATE WRITTEN
133C   890801  REVISION DATE from Version 3.2
134C   910415  Prologue converted to Version 4.0 format.  (BAB)
135C   920128  Category corrected.  (WRB)
136C   920811  Prologue revised.  (DWL)
137C   930122  Added ZEXP and ZSQRT to EXTERNAL statement.  (RWC)
138C***END PROLOGUE  ZAIRY
139C     COMPLEX AI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3
140      DOUBLE PRECISION AA, AD, AII, AIR, AK, ALIM, ATRM, AZ, AZ3, BK,
141     * CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2, DIG,
142     * DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR,
143     * S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI,
144     * ZEROR, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, ZABS, ALAZ, BB
145      INTEGER ID, IERR, IFLAG, K, KODE, K1, K2, MR, NN, NZ, I1MACH
146      DIMENSION CYR(1), CYI(1)
147      EXTERNAL ZABS, ZEXP, ZSQRT
148      DATA TTH, C1, C2, COEF /6.66666666666666667D-01,
149     * 3.55028053887817240D-01,2.58819403792806799D-01,
150     * 1.83776298473930683D-01/
151      DATA ZEROR, ZEROI, CONER, CONEI /0.0D0,0.0D0,1.0D0,0.0D0/
152C***FIRST EXECUTABLE STATEMENT  ZAIRY
153      IERR = 0
154      NZ=0
155      IF (ID.LT.0 .OR. ID.GT.1) IERR=1
156      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
157      IF (IERR.NE.0) RETURN
158      AZ = ZABS(ZR,ZI)
159      TOL = MAX(D1MACH(4),1.0D-18)
160      FID = ID
161      IF (AZ.GT.1.0D0) GO TO 70
162C-----------------------------------------------------------------------
163C     POWER SERIES FOR ABS(Z).LE.1.
164C-----------------------------------------------------------------------
165      S1R = CONER
166      S1I = CONEI
167      S2R = CONER
168      S2I = CONEI
169      IF (AZ.LT.TOL) GO TO 170
170      AA = AZ*AZ
171      IF (AA.LT.TOL/AZ) GO TO 40
172      TRM1R = CONER
173      TRM1I = CONEI
174      TRM2R = CONER
175      TRM2I = CONEI
176      ATRM = 1.0D0
177      STR = ZR*ZR - ZI*ZI
178      STI = ZR*ZI + ZI*ZR
179      Z3R = STR*ZR - STI*ZI
180      Z3I = STR*ZI + STI*ZR
181      AZ3 = AZ*AA
182      AK = 2.0D0 + FID
183      BK = 3.0D0 - FID - FID
184      CK = 4.0D0 - FID
185      DK = 3.0D0 + FID + FID
186      D1 = AK*DK
187      D2 = BK*CK
188      AD = MIN(D1,D2)
189      AK = 24.0D0 + 9.0D0*FID
190      BK = 30.0D0 - 9.0D0*FID
191      DO 30 K=1,25
192        STR = (TRM1R*Z3R-TRM1I*Z3I)/D1
193        TRM1I = (TRM1R*Z3I+TRM1I*Z3R)/D1
194        TRM1R = STR
195        S1R = S1R + TRM1R
196        S1I = S1I + TRM1I
197        STR = (TRM2R*Z3R-TRM2I*Z3I)/D2
198        TRM2I = (TRM2R*Z3I+TRM2I*Z3R)/D2
199        TRM2R = STR
200        S2R = S2R + TRM2R
201        S2I = S2I + TRM2I
202        ATRM = ATRM*AZ3/AD
203        D1 = D1 + AK
204        D2 = D2 + BK
205        AD = MIN(D1,D2)
206        IF (ATRM.LT.TOL*AD) GO TO 40
207        AK = AK + 18.0D0
208        BK = BK + 18.0D0
209   30 CONTINUE
210   40 CONTINUE
211      IF (ID.EQ.1) GO TO 50
212      AIR = S1R*C1 - C2*(ZR*S2R-ZI*S2I)
213      AII = S1I*C1 - C2*(ZR*S2I+ZI*S2R)
214      IF (KODE.EQ.1) RETURN
215      CALL ZSQRT(ZR, ZI, STR, STI)
216      ZTAR = TTH*(ZR*STR-ZI*STI)
217      ZTAI = TTH*(ZR*STI+ZI*STR)
218      CALL ZEXP(ZTAR, ZTAI, STR, STI)
219      PTR = AIR*STR - AII*STI
220      AII = AIR*STI + AII*STR
221      AIR = PTR
222      RETURN
223   50 CONTINUE
224      AIR = -S2R*C2
225      AII = -S2I*C2
226      IF (AZ.LE.TOL) GO TO 60
227      STR = ZR*S1R - ZI*S1I
228      STI = ZR*S1I + ZI*S1R
229      CC = C1/(1.0D0+FID)
230      AIR = AIR + CC*(STR*ZR-STI*ZI)
231      AII = AII + CC*(STR*ZI+STI*ZR)
232   60 CONTINUE
233      IF (KODE.EQ.1) RETURN
234      CALL ZSQRT(ZR, ZI, STR, STI)
235      ZTAR = TTH*(ZR*STR-ZI*STI)
236      ZTAI = TTH*(ZR*STI+ZI*STR)
237      CALL ZEXP(ZTAR, ZTAI, STR, STI)
238      PTR = STR*AIR - STI*AII
239      AII = STR*AII + STI*AIR
240      AIR = PTR
241      RETURN
242C-----------------------------------------------------------------------
243C     CASE FOR ABS(Z).GT.1.0
244C-----------------------------------------------------------------------
245   70 CONTINUE
246      FNU = (1.0D0+FID)/3.0D0
247C-----------------------------------------------------------------------
248C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
249C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0D-18.
250C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
251C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
252C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
253C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
254C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
255C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
256C-----------------------------------------------------------------------
257      K1 = I1MACH(15)
258      K2 = I1MACH(16)
259      R1M5 = D1MACH(5)
260      K = MIN(ABS(K1),ABS(K2))
261      ELIM = 2.303D0*(K*R1M5-3.0D0)
262      K1 = I1MACH(14) - 1
263      AA = R1M5*K1
264      DIG = MIN(AA,18.0D0)
265      AA = AA*2.303D0
266      ALIM = ELIM + MAX(-AA,-41.45D0)
267      RL = 1.2D0*DIG + 3.0D0
268      ALAZ = LOG(AZ)
269C-----------------------------------------------------------------------
270C     TEST FOR PROPER RANGE
271C-----------------------------------------------------------------------
272      AA=0.5D0/TOL
273      BB=I1MACH(9)*0.5D0
274      AA=MIN(AA,BB)
275      AA=AA**TTH
276      IF (AZ.GT.AA) GO TO 260
277      AA=SQRT(AA)
278      IF (AZ.GT.AA) IERR=3
279      CALL ZSQRT(ZR, ZI, CSQR, CSQI)
280      ZTAR = TTH*(ZR*CSQR-ZI*CSQI)
281      ZTAI = TTH*(ZR*CSQI+ZI*CSQR)
282C-----------------------------------------------------------------------
283C     RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
284C-----------------------------------------------------------------------
285      IFLAG = 0
286      SFAC = 1.0D0
287      AK = ZTAI
288      IF (ZR.GE.0.0D0) GO TO 80
289      BK = ZTAR
290      CK = -ABS(BK)
291      ZTAR = CK
292      ZTAI = AK
293   80 CONTINUE
294      IF (ZI.NE.0.0D0) GO TO 90
295      IF (ZR.GT.0.0D0) GO TO 90
296      ZTAR = 0.0D0
297      ZTAI = AK
298   90 CONTINUE
299      AA = ZTAR
300      IF (AA.GE.0.0D0 .AND. ZR.GT.0.0D0) GO TO 110
301      IF (KODE.EQ.2) GO TO 100
302C-----------------------------------------------------------------------
303C     OVERFLOW TEST
304C-----------------------------------------------------------------------
305      IF (AA.GT.(-ALIM)) GO TO 100
306      AA = -AA + 0.25D0*ALAZ
307      IFLAG = 1
308      SFAC = TOL
309      IF (AA.GT.ELIM) GO TO 270
310  100 CONTINUE
311C-----------------------------------------------------------------------
312C     CBKNU AND CACON RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
313C-----------------------------------------------------------------------
314      MR = 1
315      IF (ZI.LT.0.0D0) MR = -1
316      CALL ZACAI(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, NN, RL, TOL,
317     * ELIM, ALIM)
318      IF (NN.LT.0) GO TO 280
319      NZ = NZ + NN
320      GO TO 130
321  110 CONTINUE
322      IF (KODE.EQ.2) GO TO 120
323C-----------------------------------------------------------------------
324C     UNDERFLOW TEST
325C-----------------------------------------------------------------------
326      IF (AA.LT.ALIM) GO TO 120
327      AA = -AA - 0.25D0*ALAZ
328      IFLAG = 2
329      SFAC = 1.0D0/TOL
330      IF (AA.LT.(-ELIM)) GO TO 210
331  120 CONTINUE
332      CALL ZBKNU(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM,
333     * ALIM)
334  130 CONTINUE
335      S1R = CYR(1)*COEF
336      S1I = CYI(1)*COEF
337      IF (IFLAG.NE.0) GO TO 150
338      IF (ID.EQ.1) GO TO 140
339      AIR = CSQR*S1R - CSQI*S1I
340      AII = CSQR*S1I + CSQI*S1R
341      RETURN
342  140 CONTINUE
343      AIR = -(ZR*S1R-ZI*S1I)
344      AII = -(ZR*S1I+ZI*S1R)
345      RETURN
346  150 CONTINUE
347      S1R = S1R*SFAC
348      S1I = S1I*SFAC
349      IF (ID.EQ.1) GO TO 160
350      STR = S1R*CSQR - S1I*CSQI
351      S1I = S1R*CSQI + S1I*CSQR
352      S1R = STR
353      AIR = S1R/SFAC
354      AII = S1I/SFAC
355      RETURN
356  160 CONTINUE
357      STR = -(S1R*ZR-S1I*ZI)
358      S1I = -(S1R*ZI+S1I*ZR)
359      S1R = STR
360      AIR = S1R/SFAC
361      AII = S1I/SFAC
362      RETURN
363  170 CONTINUE
364      AA = 1.0D+3*D1MACH(1)
365      S1R = ZEROR
366      S1I = ZEROI
367      IF (ID.EQ.1) GO TO 190
368      IF (AZ.LE.AA) GO TO 180
369      S1R = C2*ZR
370      S1I = C2*ZI
371  180 CONTINUE
372      AIR = C1 - S1R
373      AII = -S1I
374      RETURN
375  190 CONTINUE
376      AIR = -C2
377      AII = 0.0D0
378      AA = SQRT(AA)
379      IF (AZ.LE.AA) GO TO 200
380      S1R = 0.5D0*(ZR*ZR-ZI*ZI)
381      S1I = ZR*ZI
382  200 CONTINUE
383      AIR = AIR + C1*S1R
384      AII = AII + C1*S1I
385      RETURN
386  210 CONTINUE
387      NZ = 1
388      AIR = ZEROR
389      AII = ZEROI
390      RETURN
391  270 CONTINUE
392      NZ = 0
393      IERR=2
394      RETURN
395  280 CONTINUE
396      IF(NN.EQ.(-1)) GO TO 270
397      NZ=0
398      IERR=5
399      RETURN
400  260 CONTINUE
401      IERR=4
402      NZ=0
403      RETURN
404      END
405