GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
double complex cmplx
Definition: Faddeeva.cc:230
subroutine cmlri(Z, FNU, KODE, N, Y, NZ, TOL)
Definition: cmlri.f:2
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137