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