1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
SUBROUTINE ZRATI(ZR, ZI, FNU, N, CYR, CYI, TOL)
C***BEGIN PROLOGUE ZRATI
C***REFER TO ZBESI,ZBESK,ZBESH
C
C ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
C RECURRENCE. THE STARTING INDEX IS DETERMINED BY FORWARD
C RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
C MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
C BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
C BY D. J. SOOKNE.
C
C***ROUTINES CALLED AZABS,ZDIV
C***END PROLOGUE ZRATI
C COMPLEX Z,CY(1),CONE,CZERO,P1,P2,T1,RZ,PT,CDFNU
DOUBLE PRECISION AK, AMAGZ, AP1, AP2, ARG, AZ, CDFNUI, CDFNUR,
* CONEI, CONER, CYI, CYR, CZEROI, CZEROR, DFNU, FDNU, FLAM, FNU,
* FNUP, PTI, PTR, P1I, P1R, P2I, P2R, RAK, RAP1, RHO, RT2, RZI,
* RZR, TEST, TEST1, TOL, TTI, TTR, T1I, T1R, ZI, ZR, AZABS
INTEGER I, ID, IDNU, INU, ITIME, K, KK, MAGZ, N
DIMENSION CYR(N), CYI(N)
DATA CZEROR,CZEROI,CONER,CONEI,RT2/
1 0.0D0, 0.0D0, 1.0D0, 0.0D0, 1.41421356237309505D0 /
AZ = AZABS(ZR,ZI)
INU = INT(SNGL(FNU))
IDNU = INU + N - 1
MAGZ = INT(SNGL(AZ))
AMAGZ = DBLE(FLOAT(MAGZ+1))
FDNU = DBLE(FLOAT(IDNU))
FNUP = DMAX1(AMAGZ,FDNU)
ID = IDNU - MAGZ - 1
ITIME = 1
K = 1
PTR = 1.0D0/AZ
RZR = PTR*(ZR+ZR)*PTR
RZI = -PTR*(ZI+ZI)*PTR
T1R = RZR*FNUP
T1I = RZI*FNUP
P2R = -T1R
P2I = -T1I
P1R = CONER
P1I = CONEI
T1R = T1R + RZR
T1I = T1I + RZI
IF (ID.GT.0) ID = 0
AP2 = AZABS(P2R,P2I)
AP1 = AZABS(P1R,P1I)
C-----------------------------------------------------------------------
C THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU
C GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
C P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
C PREMATURELY.
C-----------------------------------------------------------------------
ARG = (AP2+AP2)/(AP1*TOL)
TEST1 = DSQRT(ARG)
TEST = TEST1
RAP1 = 1.0D0/AP1
P1R = P1R*RAP1
P1I = P1I*RAP1
P2R = P2R*RAP1
P2I = P2I*RAP1
AP2 = AP2*RAP1
10 CONTINUE
K = K + 1
AP1 = AP2
PTR = P2R
PTI = P2I
P2R = P1R - (T1R*PTR-T1I*PTI)
P2I = P1I - (T1R*PTI+T1I*PTR)
P1R = PTR
P1I = PTI
T1R = T1R + RZR
T1I = T1I + RZI
AP2 = AZABS(P2R,P2I)
IF (AP1.LE.TEST) GO TO 10
IF (ITIME.EQ.2) GO TO 20
AK = AZABS(T1R,T1I)*0.5D0
FLAM = AK + DSQRT(AK*AK-1.0D0)
RHO = DMIN1(AP2/AP1,FLAM)
TEST = TEST1*DSQRT(RHO/(RHO*RHO-1.0D0))
ITIME = 2
GO TO 10
20 CONTINUE
KK = K + 1 - ID
AK = DBLE(FLOAT(KK))
T1R = AK
T1I = CZEROI
DFNU = FNU + DBLE(FLOAT(N-1))
P1R = 1.0D0/AP2
P1I = CZEROI
P2R = CZEROR
P2I = CZEROI
DO 30 I=1,KK
PTR = P1R
PTI = P1I
RAP1 = DFNU + T1R
TTR = RZR*RAP1
TTI = RZI*RAP1
P1R = (PTR*TTR-PTI*TTI) + P2R
P1I = (PTR*TTI+PTI*TTR) + P2I
P2R = PTR
P2I = PTI
T1R = T1R - CONER
30 CONTINUE
IF (P1R.NE.CZEROR .OR. P1I.NE.CZEROI) GO TO 40
P1R = TOL
P1I = TOL
40 CONTINUE
CALL ZDIV(P2R, P2I, P1R, P1I, CYR(N), CYI(N))
IF (N.EQ.1) RETURN
K = N - 1
AK = DBLE(FLOAT(K))
T1R = AK
T1I = CZEROI
CDFNUR = FNU*RZR
CDFNUI = FNU*RZI
DO 60 I=2,N
PTR = CDFNUR + (T1R*RZR-T1I*RZI) + CYR(K+1)
PTI = CDFNUI + (T1R*RZI+T1I*RZR) + CYI(K+1)
AK = AZABS(PTR,PTI)
IF (AK.NE.CZEROR) GO TO 50
PTR = TOL
PTI = TOL
AK = TOL*RT2
50 CONTINUE
RAK = CONER/AK
CYR(K) = RAK*PTR*RAK
CYI(K) = -RAK*PTI*RAK
T1R = T1R - CONER
K = K - 1
60 CONTINUE
RETURN
END
|
