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
zrati.f
Go to the documentation of this file.
1  SUBROUTINE zrati(ZR, ZI, FNU, N, CYR, CYI, TOL)
2 C***BEGIN PROLOGUE ZRATI
3 C***REFER TO ZBESI,ZBESK,ZBESH
4 C
5 C ZRATI 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 XZABS,ZDIV
13 C***END PROLOGUE ZRATI
14 C COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU
15  DOUBLE PRECISION ak, amagz, ap1, ap2, arg, az, cdfnui, cdfnur,
16  * conei, coner, cyi, cyr, czeroi, czeror, dfnu, fdnu, flam, fnu,
17  * fnup, pti, ptr, p1i, p1r, p2i, p2r, rak, rap1, rho, rt2, rzi,
18  * rzr, test, test1, tol, tti, ttr, t1i, t1r, zi, zr, xzabs
19  INTEGER i, id, idnu, inu, itime, k, kk, magz, n
20  dimension cyr(n), cyi(n)
21  DATA czeror,czeroi,coner,conei,rt2/
22  1 0.0d0, 0.0d0, 1.0d0, 0.0d0, 1.41421356237309505d0 /
23  az = xzabs(zr,zi)
24  inu = int(sngl(fnu))
25  idnu = inu + n - 1
26  magz = int(sngl(az))
27  amagz = dble(float(magz+1))
28  fdnu = dble(float(idnu))
29  fnup = dmax1(amagz,fdnu)
30  id = idnu - magz - 1
31  itime = 1
32  k = 1
33  ptr = 1.0d0/az
34  rzr = ptr*(zr+zr)*ptr
35  rzi = -ptr*(zi+zi)*ptr
36  t1r = rzr*fnup
37  t1i = rzi*fnup
38  p2r = -t1r
39  p2i = -t1i
40  p1r = coner
41  p1i = conei
42  t1r = t1r + rzr
43  t1i = t1i + rzi
44  IF (id.GT.0) id = 0
45  ap2 = xzabs(p2r,p2i)
46  ap1 = xzabs(p1r,p1i)
47 C-----------------------------------------------------------------------
48 C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU
49 C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
50 C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
51 C PREMATURELY.
52 C-----------------------------------------------------------------------
53  arg = (ap2+ap2)/(ap1*tol)
54  test1 = dsqrt(arg)
55  test = test1
56  rap1 = 1.0d0/ap1
57  p1r = p1r*rap1
58  p1i = p1i*rap1
59  p2r = p2r*rap1
60  p2i = p2i*rap1
61  ap2 = ap2*rap1
62  10 CONTINUE
63  k = k + 1
64  ap1 = ap2
65  ptr = p2r
66  pti = p2i
67  p2r = p1r - (t1r*ptr-t1i*pti)
68  p2i = p1i - (t1r*pti+t1i*ptr)
69  p1r = ptr
70  p1i = pti
71  t1r = t1r + rzr
72  t1i = t1i + rzi
73  ap2 = xzabs(p2r,p2i)
74  IF (ap1.LE.test) go to 10
75  IF (itime.EQ.2) go to 20
76  ak = xzabs(t1r,t1i)*0.5d0
77  flam = ak + dsqrt(ak*ak-1.0d0)
78  rho = dmin1(ap2/ap1,flam)
79  test = test1*dsqrt(rho/(rho*rho-1.0d0))
80  itime = 2
81  go to 10
82  20 CONTINUE
83  kk = k + 1 - id
84  ak = dble(float(kk))
85  t1r = ak
86  t1i = czeroi
87  dfnu = fnu + dble(float(n-1))
88  p1r = 1.0d0/ap2
89  p1i = czeroi
90  p2r = czeror
91  p2i = czeroi
92  DO 30 i=1,kk
93  ptr = p1r
94  pti = p1i
95  rap1 = dfnu + t1r
96  ttr = rzr*rap1
97  tti = rzi*rap1
98  p1r = (ptr*ttr-pti*tti) + p2r
99  p1i = (ptr*tti+pti*ttr) + p2i
100  p2r = ptr
101  p2i = pti
102  t1r = t1r - coner
103  30 CONTINUE
104  IF (p1r.NE.czeror .OR. p1i.NE.czeroi) go to 40
105  p1r = tol
106  p1i = tol
107  40 CONTINUE
108  CALL zdiv(p2r, p2i, p1r, p1i, cyr(n), cyi(n))
109  IF (n.EQ.1) RETURN
110  k = n - 1
111  ak = dble(float(k))
112  t1r = ak
113  t1i = czeroi
114  cdfnur = fnu*rzr
115  cdfnui = fnu*rzi
116  DO 60 i=2,n
117  ptr = cdfnur + (t1r*rzr-t1i*rzi) + cyr(k+1)
118  pti = cdfnui + (t1r*rzi+t1i*rzr) + cyi(k+1)
119  ak = xzabs(ptr,pti)
120  IF (ak.NE.czeror) go to 50
121  ptr = tol
122  pti = tol
123  ak = tol*rt2
124  50 CONTINUE
125  rak = coner/ak
126  cyr(k) = rak*ptr*rak
127  cyi(k) = -rak*pti*rak
128  t1r = t1r - coner
129  k = k - 1
130  60 CONTINUE
131  RETURN
132  END