GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zbuni.f
Go to the documentation of this file.
1  SUBROUTINE zbuni(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NUI, NLAST,
2  * fnul, tol, elim, alim)
3 C***BEGIN PROLOGUE ZBUNI
4 C***REFER TO ZBESI,ZBESK
5 C
6 C ZBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z).GT.
7 C FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM
8 C FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING
9 C ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z)
10 C ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2
11 C
12 C***ROUTINES CALLED ZUNI1,ZUNI2,XZABS,D1MACH
13 C***END PROLOGUE ZBUNI
14 C 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
31 C-----------------------------------------------------------------------
32 C ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
33 C -PI/3.LE.ARG(Z).LE.PI/3
34 C-----------------------------------------------------------------------
35  CALL zuni1(zr, zi, gnu, kode, 2, cyr, cyi, nw, nlast, fnul, tol,
36  * elim, alim)
37  go to 20
38  10 CONTINUE
39 C-----------------------------------------------------------------------
40 C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
41 C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
42 C AND HPI=PI/2
43 C-----------------------------------------------------------------------
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))
50 C----------------------------------------------------------------------
51 C SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
52 C----------------------------------------------------------------------
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
152 C-----------------------------------------------------------------------
153 C ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
154 C -PI/3.LE.ARG(Z).LE.PI/3
155 C-----------------------------------------------------------------------
156  CALL zuni1(zr, zi, fnu, kode, n, yr, yi, nw, nlast, fnul, tol,
157  * elim, alim)
158  go to 80
159  70 CONTINUE
160 C-----------------------------------------------------------------------
161 C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
162 C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
163 C AND HPI=PI/2
164 C-----------------------------------------------------------------------
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