1      SUBROUTINE ZBUNI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NUI, NLAST,
2     * FNUL, TOL, ELIM, ALIM)
3C***BEGIN PROLOGUE  ZBUNI
4C***REFER TO  ZBESI,ZBESK
5C
6C     ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z).GT.
7C     FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM
8C     FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING
9C     ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z)
10C     ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2
11C
12C***ROUTINES CALLED  ZUNI1,ZUNI2,ZABS,D1MACH
13C***END PROLOGUE  ZBUNI
14C     COMPLEX CSCL,CSCR,CY,RZ,ST,S1,S2,Y,Z
15      EXTERNAL ZABS
16      DOUBLE PRECISION ALIM, AX, AY, CSCLR, CSCRR, CYI, CYR, DFNU,
17     * ELIM, FNU, FNUI, FNUL, GNU, RAZ, RZI, RZR, STI, STR, S1I, S1R,
18     * S2I, S2R, TOL, YI, YR, ZI, ZR, ZABS, ASCLE, BRY, C1R, C1I, C1M,
19     * D1MACH
20      INTEGER I, IFLAG, IFORM, K, KODE, N, NL, NLAST, NUI, NW, NZ
21      DIMENSION YR(N), YI(N), CYR(2), CYI(2), BRY(3)
22      NZ = 0
23      AX = DABS(ZR)*1.7321D0
24      AY = DABS(ZI)
25      IFORM = 1
26      IF (AY.GT.AX) IFORM = 2
27      IF (NUI.EQ.0) GO TO 60
28      FNUI = DBLE(FLOAT(NUI))
29      DFNU = FNU + DBLE(FLOAT(N-1))
30      GNU = DFNU + FNUI
31      IF (IFORM.EQ.2) GO TO 10
32C-----------------------------------------------------------------------
33C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
34C     -PI/3.LE.ARG(Z).LE.PI/3
35C-----------------------------------------------------------------------
36      CALL ZUNI1(ZR, ZI, GNU, KODE, 2, CYR, CYI, NW, NLAST, FNUL, TOL,
37     * ELIM, ALIM)
38      GO TO 20
39   10 CONTINUE
40C-----------------------------------------------------------------------
41C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
42C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
43C     AND HPI=PI/2
44C-----------------------------------------------------------------------
45      CALL ZUNI2(ZR, ZI, GNU, KODE, 2, CYR, CYI, NW, NLAST, FNUL, TOL,
46     * ELIM, ALIM)
47   20 CONTINUE
48      IF (NW.LT.0) GO TO 50
49      IF (NW.NE.0) GO TO 90
50      STR = ZABS(CYR(1),CYI(1))
51C----------------------------------------------------------------------
52C     SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
53C----------------------------------------------------------------------
54      BRY(1)=1.0D+3*D1MACH(1)/TOL
55      BRY(2) = 1.0D0/BRY(1)
56      BRY(3) = BRY(2)
57      IFLAG = 2
58      ASCLE = BRY(2)
59      CSCLR = 1.0D0
60      IF (STR.GT.BRY(1)) GO TO 21
61      IFLAG = 1
62      ASCLE = BRY(1)
63      CSCLR = 1.0D0/TOL
64      GO TO 25
65   21 CONTINUE
66      IF (STR.LT.BRY(2)) GO TO 25
67      IFLAG = 3
68      ASCLE=BRY(3)
69      CSCLR = TOL
70   25 CONTINUE
71      CSCRR = 1.0D0/CSCLR
72      S1R = CYR(2)*CSCLR
73      S1I = CYI(2)*CSCLR
74      S2R = CYR(1)*CSCLR
75      S2I = CYI(1)*CSCLR
76      RAZ = 1.0D0/ZABS(ZR,ZI)
77      STR = ZR*RAZ
78      STI = -ZI*RAZ
79      RZR = (STR+STR)*RAZ
80      RZI = (STI+STI)*RAZ
81      DO 30 I=1,NUI
82        STR = S2R
83        STI = S2I
84        S2R = (DFNU+FNUI)*(RZR*STR-RZI*STI) + S1R
85        S2I = (DFNU+FNUI)*(RZR*STI+RZI*STR) + S1I
86        S1R = STR
87        S1I = STI
88        FNUI = FNUI - 1.0D0
89        IF (IFLAG.GE.3) GO TO 30
90        STR = S2R*CSCRR
91        STI = S2I*CSCRR
92        C1R = DABS(STR)
93        C1I = DABS(STI)
94        C1M = DMAX1(C1R,C1I)
95        IF (C1M.LE.ASCLE) GO TO 30
96        IFLAG = IFLAG+1
97        ASCLE = BRY(IFLAG)
98        S1R = S1R*CSCRR
99        S1I = S1I*CSCRR
100        S2R = STR
101        S2I = STI
102        CSCLR = CSCLR*TOL
103        CSCRR = 1.0D0/CSCLR
104        S1R = S1R*CSCLR
105        S1I = S1I*CSCLR
106        S2R = S2R*CSCLR
107        S2I = S2I*CSCLR
108   30 CONTINUE
109      YR(N) = S2R*CSCRR
110      YI(N) = S2I*CSCRR
111      IF (N.EQ.1) RETURN
112      NL = N - 1
113      FNUI = DBLE(FLOAT(NL))
114      K = NL
115      DO 40 I=1,NL
116        STR = S2R
117        STI = S2I
118        S2R = (FNU+FNUI)*(RZR*STR-RZI*STI) + S1R
119        S2I = (FNU+FNUI)*(RZR*STI+RZI*STR) + S1I
120        S1R = STR
121        S1I = STI
122        STR = S2R*CSCRR
123        STI = S2I*CSCRR
124        YR(K) = STR
125        YI(K) = STI
126        FNUI = FNUI - 1.0D0
127        K = K - 1
128        IF (IFLAG.GE.3) GO TO 40
129        C1R = DABS(STR)
130        C1I = DABS(STI)
131        C1M = DMAX1(C1R,C1I)
132        IF (C1M.LE.ASCLE) GO TO 40
133        IFLAG = IFLAG+1
134        ASCLE = BRY(IFLAG)
135        S1R = S1R*CSCRR
136        S1I = S1I*CSCRR
137        S2R = STR
138        S2I = STI
139        CSCLR = CSCLR*TOL
140        CSCRR = 1.0D0/CSCLR
141        S1R = S1R*CSCLR
142        S1I = S1I*CSCLR
143        S2R = S2R*CSCLR
144        S2I = S2I*CSCLR
145   40 CONTINUE
146      RETURN
147   50 CONTINUE
148      NZ = -1
149      IF(NW.EQ.(-2)) NZ=-2
150      RETURN
151   60 CONTINUE
152      IF (IFORM.EQ.2) GO TO 70
153C-----------------------------------------------------------------------
154C     ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
155C     -PI/3.LE.ARG(Z).LE.PI/3
156C-----------------------------------------------------------------------
157      CALL ZUNI1(ZR, ZI, FNU, KODE, N, YR, YI, NW, NLAST, FNUL, TOL,
158     * ELIM, ALIM)
159      GO TO 80
160   70 CONTINUE
161C-----------------------------------------------------------------------
162C     ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
163C     APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
164C     AND HPI=PI/2
165C-----------------------------------------------------------------------
166      CALL ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NW, NLAST, FNUL, TOL,
167     * ELIM, ALIM)
168   80 CONTINUE
169      IF (NW.LT.0) GO TO 50
170      NZ = NW
171      RETURN
172   90 CONTINUE
173      NLAST = N
174      RETURN
175      END
176