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