GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
subroutine zdiv(AR, AI, BR, BI, CR, CI)
Definition: zdiv.f:2
subroutine zrati(ZR, ZI, FNU, N, CYR, CYI, TOL)
Definition: zrati.f:2