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