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