GNU Octave  4.2.1 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
subroutine xzlog(AR, AI, BR, BI, IERR)
Definition: xzlog.f:1
subroutine xzexp(AR, AI, BR, BI)
Definition: xzexp.f:1
double precision function d1mach(i)
Definition: d1mach.f:1
std::string str
Definition: hash.cc:118
double precision function dgamln(Z, IERR)
Definition: dgamln.f:1
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to