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
zmlri.f
Go to the documentation of this file.
1  SUBROUTINE zmlri(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
2 C***BEGIN PROLOGUE ZMLRI
3 C***REFER TO ZBESI,ZBESK
4 C
5 C ZMLRI 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 DGAMLN,D1MACH,XZABS,XZEXP,XZLOG,ZMLT
9 C***END PROLOGUE ZMLRI
10 C COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
11  DOUBLE PRECISION ack, ak, ap, at, az, bk, cki, ckr, cnormi,
12  * cnormr, conei, coner, fkap, fkk, flam, fnf, fnu, pti, ptr, p1i,
13  * p1r, p2i, p2r, raz, rho, rho2, rzi, rzr, scle, sti, str, sumi,
14  * sumr, tfnf, tol, tst, yi, yr, zeroi, zeror, zi, zr, dgamln,
15  * d1mach, xzabs
16  INTEGER i, iaz, idum, ifnu, inu, itime, k, kk, km, kode, m, n, nz
17  dimension yr(n), yi(n)
18  DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
19  scle = d1mach(1)/tol
20  nz=0
21  az = xzabs(zr,zi)
22  iaz = int(sngl(az))
23  ifnu = int(sngl(fnu))
24  inu = ifnu + n - 1
25  at = dble(float(iaz)) + 1.0d0
26  raz = 1.0d0/az
27  str = zr*raz
28  sti = -zi*raz
29  ckr = str*at*raz
30  cki = sti*at*raz
31  rzr = (str+str)*raz
32  rzi = (sti+sti)*raz
33  p1r = zeror
34  p1i = zeroi
35  p2r = coner
36  p2i = conei
37  ack = (at+1.0d0)*raz
38  rho = ack + dsqrt(ack*ack-1.0d0)
39  rho2 = rho*rho
40  tst = (rho2+rho2)/((rho2-1.0d0)*(rho-1.0d0))
41  tst = tst/tol
42 C-----------------------------------------------------------------------
43 C COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
44 C-----------------------------------------------------------------------
45  ak = at
46  DO 10 i=1,80
47  ptr = p2r
48  pti = p2i
49  p2r = p1r - (ckr*ptr-cki*pti)
50  p2i = p1i - (cki*ptr+ckr*pti)
51  p1r = ptr
52  p1i = pti
53  ckr = ckr + rzr
54  cki = cki + rzi
55  ap = xzabs(p2r,p2i)
56  IF (ap.GT.tst*ak*ak) go to 20
57  ak = ak + 1.0d0
58  10 CONTINUE
59  go to 110
60  20 CONTINUE
61  i = i + 1
62  k = 0
63  IF (inu.LT.iaz) go to 40
64 C-----------------------------------------------------------------------
65 C COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
66 C-----------------------------------------------------------------------
67  p1r = zeror
68  p1i = zeroi
69  p2r = coner
70  p2i = conei
71  at = dble(float(inu)) + 1.0d0
72  str = zr*raz
73  sti = -zi*raz
74  ckr = str*at*raz
75  cki = sti*at*raz
76  ack = at*raz
77  tst = dsqrt(ack/tol)
78  itime = 1
79  DO 30 k=1,80
80  ptr = p2r
81  pti = p2i
82  p2r = p1r - (ckr*ptr-cki*pti)
83  p2i = p1i - (ckr*pti+cki*ptr)
84  p1r = ptr
85  p1i = pti
86  ckr = ckr + rzr
87  cki = cki + rzi
88  ap = xzabs(p2r,p2i)
89  IF (ap.LT.tst) go to 30
90  IF (itime.EQ.2) go to 40
91  ack = xzabs(ckr,cki)
92  flam = ack + dsqrt(ack*ack-1.0d0)
93  fkap = ap/xzabs(p1r,p1i)
94  rho = dmin1(flam,fkap)
95  tst = tst*dsqrt(rho/(rho*rho-1.0d0))
96  itime = 2
97  30 CONTINUE
98  go to 110
99  40 CONTINUE
100 C-----------------------------------------------------------------------
101 C BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
102 C-----------------------------------------------------------------------
103  k = k + 1
104  kk = max0(i+iaz,k+inu)
105  fkk = dble(float(kk))
106  p1r = zeror
107  p1i = zeroi
108 C-----------------------------------------------------------------------
109 C SCALE P2 AND SUM BY SCLE
110 C-----------------------------------------------------------------------
111  p2r = scle
112  p2i = zeroi
113  fnf = fnu - dble(float(ifnu))
114  tfnf = fnf + fnf
115  bk = dgamln(fkk+tfnf+1.0d0,idum) - dgamln(fkk+1.0d0,idum) -
116  * dgamln(tfnf+1.0d0,idum)
117  bk = dexp(bk)
118  sumr = zeror
119  sumi = zeroi
120  km = kk - inu
121  DO 50 i=1,km
122  ptr = p2r
123  pti = p2i
124  p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
125  p2i = p1i + (fkk+fnf)*(rzi*ptr+rzr*pti)
126  p1r = ptr
127  p1i = pti
128  ak = 1.0d0 - tfnf/(fkk+tfnf)
129  ack = bk*ak
130  sumr = sumr + (ack+bk)*p1r
131  sumi = sumi + (ack+bk)*p1i
132  bk = ack
133  fkk = fkk - 1.0d0
134  50 CONTINUE
135  yr(n) = p2r
136  yi(n) = p2i
137  IF (n.EQ.1) go to 70
138  DO 60 i=2,n
139  ptr = p2r
140  pti = p2i
141  p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
142  p2i = p1i + (fkk+fnf)*(rzi*ptr+rzr*pti)
143  p1r = ptr
144  p1i = pti
145  ak = 1.0d0 - tfnf/(fkk+tfnf)
146  ack = bk*ak
147  sumr = sumr + (ack+bk)*p1r
148  sumi = sumi + (ack+bk)*p1i
149  bk = ack
150  fkk = fkk - 1.0d0
151  m = n - i + 1
152  yr(m) = p2r
153  yi(m) = p2i
154  60 CONTINUE
155  70 CONTINUE
156  IF (ifnu.LE.0) go to 90
157  DO 80 i=1,ifnu
158  ptr = p2r
159  pti = p2i
160  p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
161  p2i = p1i + (fkk+fnf)*(rzr*pti+rzi*ptr)
162  p1r = ptr
163  p1i = pti
164  ak = 1.0d0 - tfnf/(fkk+tfnf)
165  ack = bk*ak
166  sumr = sumr + (ack+bk)*p1r
167  sumi = sumi + (ack+bk)*p1i
168  bk = ack
169  fkk = fkk - 1.0d0
170  80 CONTINUE
171  90 CONTINUE
172  ptr = zr
173  pti = zi
174  IF (kode.EQ.2) ptr = zeror
175  CALL xzlog(rzr, rzi, str, sti, idum)
176  p1r = -fnf*str + ptr
177  p1i = -fnf*sti + pti
178  ap = dgamln(1.0d0+fnf,idum)
179  ptr = p1r - ap
180  pti = p1i
181 C-----------------------------------------------------------------------
182 C THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
183 C IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
184 C-----------------------------------------------------------------------
185  p2r = p2r + sumr
186  p2i = p2i + sumi
187  ap = xzabs(p2r,p2i)
188  p1r = 1.0d0/ap
189  CALL xzexp(ptr, pti, str, sti)
190  ckr = str*p1r
191  cki = sti*p1r
192  ptr = p2r*p1r
193  pti = -p2i*p1r
194  CALL zmlt(ckr, cki, ptr, pti, cnormr, cnormi)
195  DO 100 i=1,n
196  str = yr(i)*cnormr - yi(i)*cnormi
197  yi(i) = yr(i)*cnormi + yi(i)*cnormr
198  yr(i) = str
199  100 CONTINUE
200  RETURN
201  110 CONTINUE
202  nz=-2
203  RETURN
204  END