GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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