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
cmlri.f
Go to the documentation of this file.
1  SUBROUTINE cmlri(Z, FNU, KODE, N, Y, NZ, TOL)
2 C***BEGIN PROLOGUE CMLRI
3 C***REFER TO CBESI,CBESK
4 C
5 C CMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
6 C MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
7 C
8 C***ROUTINES CALLED GAMLN,R1MACH
9 C***END PROLOGUE CMLRI
10  COMPLEX ck, cnorm, cone, ctwo, czero, pt, p1, p2, rz, sum, y, z
11  REAL ack, ak, ap, at, az, bk, fkap, fkk, flam, fnf, fnu, rho,
12  * rho2, scle, tfnf, tol, tst, x, gamln, r1mach
13  INTEGER i, iaz, idum, ifnu, inu, itime, k, kk, km, kode, m, n
14  dimension y(n)
15  DATA czero,cone,ctwo /(0.0e0,0.0e0),(1.0e0,0.0e0),(2.0e0,0.0e0)/
16  scle = 1.0e+3*r1mach(1)/tol
17  nz=0
18  az = cabs(z)
19  x = REAL(z)
20  iaz = int(az)
21  ifnu = int(fnu)
22  inu = ifnu + n - 1
23  at = float(iaz) + 1.0e0
24  ck = cmplx(at,0.0e0)/z
25  rz = ctwo/z
26  p1 = czero
27  p2 = cone
28  ack = (at+1.0e0)/az
29  rho = ack + sqrt(ack*ack-1.0e0)
30  rho2 = rho*rho
31  tst = (rho2+rho2)/((rho2-1.0e0)*(rho-1.0e0))
32  tst = tst/tol
33 C-----------------------------------------------------------------------
34 C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
35 C-----------------------------------------------------------------------
36  ak = at
37  DO 10 i=1,80
38  pt = p2
39  p2 = p1 - ck*p2
40  p1 = pt
41  ck = ck + rz
42  ap = cabs(p2)
43  IF (ap.GT.tst*ak*ak) go to 20
44  ak = ak + 1.0e0
45  10 CONTINUE
46  go to 110
47  20 CONTINUE
48  i = i + 1
49  k = 0
50  IF (inu.LT.iaz) go to 40
51 C-----------------------------------------------------------------------
52 C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
53 C-----------------------------------------------------------------------
54  p1 = czero
55  p2 = cone
56  at = float(inu) + 1.0e0
57  ck = cmplx(at,0.0e0)/z
58  ack = at/az
59  tst = sqrt(ack/tol)
60  itime = 1
61  DO 30 k=1,80
62  pt = p2
63  p2 = p1 - ck*p2
64  p1 = pt
65  ck = ck + rz
66  ap = cabs(p2)
67  IF (ap.LT.tst) go to 30
68  IF (itime.EQ.2) go to 40
69  ack = cabs(ck)
70  flam = ack + sqrt(ack*ack-1.0e0)
71  fkap = ap/cabs(p1)
72  rho = amin1(flam,fkap)
73  tst = tst*sqrt(rho/(rho*rho-1.0e0))
74  itime = 2
75  30 CONTINUE
76  go to 110
77  40 CONTINUE
78 C-----------------------------------------------------------------------
79 C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
80 C-----------------------------------------------------------------------
81  k = k + 1
82  kk = max0(i+iaz,k+inu)
83  fkk = float(kk)
84  p1 = czero
85 C-----------------------------------------------------------------------
86 C SCALE P2 AND SUM BY SCLE
87 C-----------------------------------------------------------------------
88  p2 = cmplx(scle,0.0e0)
89  fnf = fnu - float(ifnu)
90  tfnf = fnf + fnf
91  bk = gamln(fkk+tfnf+1.0e0,idum) - gamln(fkk+1.0e0,idum)
92  * -gamln(tfnf+1.0e0,idum)
93  bk = exp(bk)
94  sum = czero
95  km = kk - inu
96  DO 50 i=1,km
97  pt = p2
98  p2 = p1 + cmplx(fkk+fnf,0.0e0)*rz*p2
99  p1 = pt
100  ak = 1.0e0 - tfnf/(fkk+tfnf)
101  ack = bk*ak
102  sum = sum + cmplx(ack+bk,0.0e0)*p1
103  bk = ack
104  fkk = fkk - 1.0e0
105  50 CONTINUE
106  y(n) = p2
107  IF (n.EQ.1) go to 70
108  DO 60 i=2,n
109  pt = p2
110  p2 = p1 + cmplx(fkk+fnf,0.0e0)*rz*p2
111  p1 = pt
112  ak = 1.0e0 - tfnf/(fkk+tfnf)
113  ack = bk*ak
114  sum = sum + cmplx(ack+bk,0.0e0)*p1
115  bk = ack
116  fkk = fkk - 1.0e0
117  m = n - i + 1
118  y(m) = p2
119  60 CONTINUE
120  70 CONTINUE
121  IF (ifnu.LE.0) go to 90
122  DO 80 i=1,ifnu
123  pt = p2
124  p2 = p1 + cmplx(fkk+fnf,0.0e0)*rz*p2
125  p1 = pt
126  ak = 1.0e0 - tfnf/(fkk+tfnf)
127  ack = bk*ak
128  sum = sum + cmplx(ack+bk,0.0e0)*p1
129  bk = ack
130  fkk = fkk - 1.0e0
131  80 CONTINUE
132  90 CONTINUE
133  pt = z
134  IF (kode.EQ.2) pt = pt - cmplx(x,0.0e0)
135  p1 = -cmplx(fnf,0.0e0)*clog(rz) + pt
136  ap = gamln(1.0e0+fnf,idum)
137  pt = p1 - cmplx(ap,0.0e0)
138 C-----------------------------------------------------------------------
139 C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
140 C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
141 C-----------------------------------------------------------------------
142  p2 = p2 + sum
143  ap = cabs(p2)
144  p1 = cmplx(1.0e0/ap,0.0e0)
145  ck = cexp(pt)*p1
146  pt = conjg(p2)*p1
147  cnorm = ck*pt
148  DO 100 i=1,n
149  y(i) = y(i)*cnorm
150  100 CONTINUE
151  RETURN
152  110 CONTINUE
153  nz=-2
154  RETURN
155  END