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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
|
SUBROUTINE ZASYI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM,
* ALIM)
C***BEGIN PROLOGUE ZASYI
C***REFER TO ZBESI,ZBESK
C
C ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
C REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
C
C***ROUTINES CALLED D1MACH,AZABS,ZDIV,AZEXP,ZMLT,AZSQRT
C***END PROLOGUE ZASYI
C COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z
DOUBLE PRECISION AA, AEZ, AK, AK1I, AK1R, ALIM, ARG, ARM, ATOL,
* AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI,
* CZR, DFNU, DKI, DKR, DNU2, ELIM, EZI, EZR, FDN, FNU, PI, P1I,
* P1R, RAZ, RL, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I,
* S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, AZABS
INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
DIMENSION YR(N), YI(N)
DATA PI, RTPI /3.14159265358979324D0 , 0.159154943091895336D0 /
DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
C
NZ = 0
AZ = AZABS(ZR,ZI)
ARM = 1.0D+3*D1MACH(1)
RTR1 = DSQRT(ARM)
IL = MIN0(2,N)
DFNU = FNU + DBLE(FLOAT(N-IL))
C-----------------------------------------------------------------------
C OVERFLOW TEST
C-----------------------------------------------------------------------
RAZ = 1.0D0/AZ
STR = ZR*RAZ
STI = -ZI*RAZ
AK1R = RTPI*STR*RAZ
AK1I = RTPI*STI*RAZ
CALL AZSQRT(AK1R, AK1I, AK1R, AK1I)
CZR = ZR
CZI = ZI
IF (KODE.NE.2) GO TO 10
CZR = ZEROR
CZI = ZI
10 CONTINUE
IF (DABS(CZR).GT.ELIM) GO TO 100
DNU2 = DFNU + DFNU
KODED = 1
IF ((DABS(CZR).GT.ALIM) .AND. (N.GT.2)) GO TO 20
KODED = 0
CALL AZEXP(CZR, CZI, STR, STI)
CALL ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I)
20 CONTINUE
FDN = 0.0D0
IF (DNU2.GT.RTR1) FDN = DNU2*DNU2
EZR = ZR*8.0D0
EZI = ZI*8.0D0
C-----------------------------------------------------------------------
C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
C EXPANSION FOR THE IMAGINARY PART.
C-----------------------------------------------------------------------
AEZ = 8.0D0*AZ
S = TOL/AEZ
JL = INT(SNGL(RL+RL)) + 2
P1R = ZEROR
P1I = ZEROI
IF (ZI.EQ.0.0D0) GO TO 30
C-----------------------------------------------------------------------
C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
C SIGNIFICANCE WHEN FNU OR N IS LARGE
C-----------------------------------------------------------------------
INU = INT(SNGL(FNU))
ARG = (FNU-DBLE(FLOAT(INU)))*PI
INU = INU + N - IL
AK = -DSIN(ARG)
BK = DCOS(ARG)
IF (ZI.LT.0.0D0) BK = -BK
P1R = AK
P1I = BK
IF (MOD(INU,2).EQ.0) GO TO 30
P1R = -P1R
P1I = -P1I
30 CONTINUE
DO 70 K=1,IL
SQK = FDN - 1.0D0
ATOL = S*DABS(SQK)
SGN = 1.0D0
CS1R = CONER
CS1I = CONEI
CS2R = CONER
CS2I = CONEI
CKR = CONER
CKI = CONEI
AK = 0.0D0
AA = 1.0D0
BB = AEZ
DKR = EZR
DKI = EZI
DO 40 J=1,JL
CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI)
CKR = STR*SQK
CKI = STI*SQK
CS2R = CS2R + CKR
CS2I = CS2I + CKI
SGN = -SGN
CS1R = CS1R + CKR*SGN
CS1I = CS1I + CKI*SGN
DKR = DKR + EZR
DKI = DKI + EZI
AA = AA*DABS(SQK)/BB
BB = BB + AEZ
AK = AK + 8.0D0
SQK = SQK - AK
IF (AA.LE.ATOL) GO TO 50
40 CONTINUE
GO TO 110
50 CONTINUE
S2R = CS1R
S2I = CS1I
IF (ZR+ZR.GE.ELIM) GO TO 60
TZR = ZR + ZR
TZI = ZI + ZI
CALL AZEXP(-TZR, -TZI, STR, STI)
CALL ZMLT(STR, STI, P1R, P1I, STR, STI)
CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI)
S2R = S2R + STR
S2I = S2I + STI
60 CONTINUE
FDN = FDN + 8.0D0*DFNU + 4.0D0
P1R = -P1R
P1I = -P1I
M = N - IL + K
YR(M) = S2R*AK1R - S2I*AK1I
YI(M) = S2R*AK1I + S2I*AK1R
70 CONTINUE
IF (N.LE.2) RETURN
NN = N
K = NN - 2
AK = DBLE(FLOAT(K))
STR = ZR*RAZ
STI = -ZI*RAZ
RZR = (STR+STR)*RAZ
RZI = (STI+STI)*RAZ
IB = 3
DO 80 I=IB,NN
YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2)
YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2)
AK = AK - 1.0D0
K = K - 1
80 CONTINUE
IF (KODED.EQ.0) RETURN
CALL AZEXP(CZR, CZI, CKR, CKI)
DO 90 I=1,NN
STR = YR(I)*CKR - YI(I)*CKI
YI(I) = YR(I)*CKI + YI(I)*CKR
YR(I) = STR
90 CONTINUE
RETURN
100 CONTINUE
NZ = -1
RETURN
110 CONTINUE
NZ=-2
RETURN
END
|
