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
crati.f
Go to the documentation of this file.
1  SUBROUTINE crati(Z, FNU, N, CY, TOL)
2 C***BEGIN PROLOGUE CRATI
3 C***REFER TO CBESI,CBESK,CBESH
4 C
5 C CRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
6 C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD
7 C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
8 C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
9 C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
10 C BY D. J. SOOKNE.
11 C
12 C***ROUTINES CALLED (NONE)
13 C***END PROLOGUE CRATI
14  COMPLEX cdfnu, cone, cy, czero, pt, p1, p2, rz, t1, z
15  REAL ak, amagz, ap1, ap2, arg, az, dfnu, fdnu, flam, fnu, fnup,
16  * rap1, rho, test, test1, tol
17  INTEGER i, id, idnu, inu, itime, k, kk, magz, n
18  dimension cy(n)
19  DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
20  az = cabs(z)
21  inu = int(fnu)
22  idnu = inu + n - 1
23  fdnu = float(idnu)
24  magz = int(az)
25  amagz = float(magz+1)
26  fnup = amax1(amagz,fdnu)
27  id = idnu - magz - 1
28  itime = 1
29  k = 1
30  rz = (cone+cone)/z
31  t1 = cmplx(fnup,0.0e0)*rz
32  p2 = -t1
33  p1 = cone
34  t1 = t1 + rz
35  IF (id.GT.0) id = 0
36  ap2 = cabs(p2)
37  ap1 = cabs(p1)
38 C-----------------------------------------------------------------------
39 C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNX
40 C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
41 C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
42 C PREMATURELY.
43 C-----------------------------------------------------------------------
44  arg = (ap2+ap2)/(ap1*tol)
45  test1 = sqrt(arg)
46  test = test1
47  rap1 = 1.0e0/ap1
48  p1 = p1*cmplx(rap1,0.0e0)
49  p2 = p2*cmplx(rap1,0.0e0)
50  ap2 = ap2*rap1
51  10 CONTINUE
52  k = k + 1
53  ap1 = ap2
54  pt = p2
55  p2 = p1 - t1*p2
56  p1 = pt
57  t1 = t1 + rz
58  ap2 = cabs(p2)
59  IF (ap1.LE.test) go to 10
60  IF (itime.EQ.2) go to 20
61  ak = cabs(t1)*0.5e0
62  flam = ak + sqrt(ak*ak-1.0e0)
63  rho = amin1(ap2/ap1,flam)
64  test = test1*sqrt(rho/(rho*rho-1.0e0))
65  itime = 2
66  go to 10
67  20 CONTINUE
68  kk = k + 1 - id
69  ak = float(kk)
70  dfnu = fnu + float(n-1)
71  cdfnu = cmplx(dfnu,0.0e0)
72  t1 = cmplx(ak,0.0e0)
73  p1 = cmplx(1.0e0/ap2,0.0e0)
74  p2 = czero
75  DO 30 i=1,kk
76  pt = p1
77  p1 = rz*(cdfnu+t1)*p1 + p2
78  p2 = pt
79  t1 = t1 - cone
80  30 CONTINUE
81  IF (REAL(p1).NE.0.0e0 .OR. aimag(p1).NE.0.0e0) go to 40
82  p1 = cmplx(tol,tol)
83  40 CONTINUE
84  cy(n) = p2/p1
85  IF (n.EQ.1) RETURN
86  k = n - 1
87  ak = float(k)
88  t1 = cmplx(ak,0.0e0)
89  cdfnu = cmplx(fnu,0.0e0)*rz
90  DO 60 i=2,n
91  pt = cdfnu + t1*rz + cy(k+1)
92  IF (REAL(pt).NE.0.0e0 .OR. aimag(pt).NE.0.0e0) go to 50
93  pt = cmplx(tol,tol)
94  50 CONTINUE
95  cy(k) = cone/pt
96  t1 = t1 - cone
97  k = k - 1
98  60 CONTINUE
99  RETURN
100  END