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
cbuni.f
Go to the documentation of this file.
1  SUBROUTINE cbuni(Z, FNU, KODE, N, Y, NZ, NUI, NLAST, FNUL, TOL,
2  * elim, alim)
3 C***BEGIN PROLOGUE CBUNI
4 C***REFER TO CBESI,CBESK
5 C
6 C CBUNI 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 CUNI1,CUNI2,R1MACH
13 C***END PROLOGUE CBUNI
14  COMPLEX cscl, cscr, cy, rz, st, s1, s2, y, z
15  REAL alim, ax, ay, dfnu, elim, fnu, fnui, fnul, gnu, tol, xx, yy,
16  * ascle, bry, str, sti, stm, r1mach
17  INTEGER i, iflag, iform, k, kode, n, nl, nlast, nui, nw, nz
18  dimension y(n), cy(2), bry(3)
19  nz = 0
20  xx = REAL(z)
21  yy = aimag(z)
22  ax = abs(xx)*1.7321e0
23  ay = abs(yy)
24  iform = 1
25  IF (ay.GT.ax) iform = 2
26  IF (nui.EQ.0) go to 60
27  fnui = float(nui)
28  dfnu = fnu + 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 cuni1(z, gnu, kode, 2, cy, nw, nlast, fnul, tol, elim, alim)
36  go to 20
37  10 CONTINUE
38 C-----------------------------------------------------------------------
39 C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
40 C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
41 C AND HPI=PI/2
42 C-----------------------------------------------------------------------
43  CALL cuni2(z, gnu, kode, 2, cy, nw, nlast, fnul, tol, elim, alim)
44  20 CONTINUE
45  IF (nw.LT.0) go to 50
46  IF (nw.NE.0) go to 90
47  ay = cabs(cy(1))
48 C----------------------------------------------------------------------
49 C SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
50 C----------------------------------------------------------------------
51  bry(1) = 1.0e+3*r1mach(1)/tol
52  bry(2) = 1.0e0/bry(1)
53  bry(3) = bry(2)
54  iflag = 2
55  ascle = bry(2)
56  ax = 1.0e0
57  cscl = cmplx(ax,0.0e0)
58  IF (ay.GT.bry(1)) go to 21
59  iflag = 1
60  ascle = bry(1)
61  ax = 1.0e0/tol
62  cscl = cmplx(ax,0.0e0)
63  go to 25
64  21 CONTINUE
65  IF (ay.LT.bry(2)) go to 25
66  iflag = 3
67  ascle = bry(3)
68  ax = tol
69  cscl = cmplx(ax,0.0e0)
70  25 CONTINUE
71  ay = 1.0e0/ax
72  cscr = cmplx(ay,0.0e0)
73  s1 = cy(2)*cscl
74  s2 = cy(1)*cscl
75  rz = cmplx(2.0e0,0.0e0)/z
76  DO 30 i=1,nui
77  st = s2
78  s2 = cmplx(dfnu+fnui,0.0e0)*rz*s2 + s1
79  s1 = st
80  fnui = fnui - 1.0e0
81  IF (iflag.GE.3) go to 30
82  st = s2*cscr
83  str = REAL(st)
84  sti = aimag(st)
85  str = abs(str)
86  sti = abs(sti)
87  stm = amax1(str,sti)
88  IF (stm.LE.ascle) go to 30
89  iflag = iflag+1
90  ascle = bry(iflag)
91  s1 = s1*cscr
92  s2 = st
93  ax = ax*tol
94  ay = 1.0e0/ax
95  cscl = cmplx(ax,0.0e0)
96  cscr = cmplx(ay,0.0e0)
97  s1 = s1*cscl
98  s2 = s2*cscl
99  30 CONTINUE
100  y(n) = s2*cscr
101  IF (n.EQ.1) RETURN
102  nl = n - 1
103  fnui = float(nl)
104  k = nl
105  DO 40 i=1,nl
106  st = s2
107  s2 = cmplx(fnu+fnui,0.0e0)*rz*s2 + s1
108  s1 = st
109  st = s2*cscr
110  y(k) = st
111  fnui = fnui - 1.0e0
112  k = k - 1
113  IF (iflag.GE.3) go to 40
114  str = REAL(st)
115  sti = aimag(st)
116  str = abs(str)
117  sti = abs(sti)
118  stm = amax1(str,sti)
119  IF (stm.LE.ascle) go to 40
120  iflag = iflag+1
121  ascle = bry(iflag)
122  s1 = s1*cscr
123  s2 = st
124  ax = ax*tol
125  ay = 1.0e0/ax
126  cscl = cmplx(ax,0.0e0)
127  cscr = cmplx(ay,0.0e0)
128  s1 = s1*cscl
129  s2 = s2*cscl
130  40 CONTINUE
131  RETURN
132  50 CONTINUE
133  nz = -1
134  IF(nw.EQ.(-2)) nz=-2
135  RETURN
136  60 CONTINUE
137  IF (iform.EQ.2) go to 70
138 C-----------------------------------------------------------------------
139 C ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
140 C -PI/3.LE.ARG(Z).LE.PI/3
141 C-----------------------------------------------------------------------
142  CALL cuni1(z, fnu, kode, n, y, nw, nlast, fnul, tol, elim, alim)
143  go to 80
144  70 CONTINUE
145 C-----------------------------------------------------------------------
146 C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
147 C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
148 C AND HPI=PI/2
149 C-----------------------------------------------------------------------
150  CALL cuni2(z, fnu, kode, n, y, nw, nlast, fnul, tol, elim, alim)
151  80 CONTINUE
152  IF (nw.LT.0) go to 50
153  nz = nw
154  RETURN
155  90 CONTINUE
156  nlast = n
157  RETURN
158  END