|
3217
|
1 SUBROUTINE ZASYI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, |
|
|
2 * ALIM) |
|
|
3 C***BEGIN PROLOGUE ZASYI |
|
|
4 C***REFER TO ZBESI,ZBESK |
|
|
5 C |
|
|
6 C ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY |
|
|
7 C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE |
|
|
8 C REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN. |
|
|
9 C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1. |
|
|
10 C |
|
|
11 C***ROUTINES CALLED D1MACH,XZABS,ZDIV,XZEXP,ZMLT,XZSQRT |
|
|
12 C***END PROLOGUE ZASYI |
|
|
13 C COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z |
|
|
14 DOUBLE PRECISION AA, AEZ, AK, AK1I, AK1R, ALIM, ARG, ARM, ATOL, |
|
|
15 * AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI, |
|
|
16 * CZR, DFNU, DKI, DKR, DNU2, ELIM, EZI, EZR, FDN, FNU, PI, P1I, |
|
|
17 * P1R, RAZ, RL, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I, |
|
|
18 * S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, XZABS |
|
|
19 INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ |
|
|
20 DIMENSION YR(N), YI(N) |
|
|
21 DATA PI, RTPI /3.14159265358979324D0 , 0.159154943091895336D0 / |
|
|
22 DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 / |
|
|
23 C |
|
|
24 NZ = 0 |
|
|
25 AZ = XZABS(ZR,ZI) |
|
|
26 ARM = 1.0D+3*D1MACH(1) |
|
|
27 RTR1 = DSQRT(ARM) |
|
|
28 IL = MIN0(2,N) |
|
|
29 DFNU = FNU + DBLE(FLOAT(N-IL)) |
|
|
30 C----------------------------------------------------------------------- |
|
|
31 C OVERFLOW TEST |
|
|
32 C----------------------------------------------------------------------- |
|
|
33 RAZ = 1.0D0/AZ |
|
|
34 STR = ZR*RAZ |
|
|
35 STI = -ZI*RAZ |
|
|
36 AK1R = RTPI*STR*RAZ |
|
|
37 AK1I = RTPI*STI*RAZ |
|
|
38 CALL XZSQRT(AK1R, AK1I, AK1R, AK1I) |
|
|
39 CZR = ZR |
|
|
40 CZI = ZI |
|
|
41 IF (KODE.NE.2) GO TO 10 |
|
|
42 CZR = ZEROR |
|
|
43 CZI = ZI |
|
|
44 10 CONTINUE |
|
|
45 IF (DABS(CZR).GT.ELIM) GO TO 100 |
|
|
46 DNU2 = DFNU + DFNU |
|
|
47 KODED = 1 |
|
|
48 IF ((DABS(CZR).GT.ALIM) .AND. (N.GT.2)) GO TO 20 |
|
|
49 KODED = 0 |
|
|
50 CALL XZEXP(CZR, CZI, STR, STI) |
|
|
51 CALL ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I) |
|
|
52 20 CONTINUE |
|
|
53 FDN = 0.0D0 |
|
|
54 IF (DNU2.GT.RTR1) FDN = DNU2*DNU2 |
|
|
55 EZR = ZR*8.0D0 |
|
|
56 EZI = ZI*8.0D0 |
|
|
57 C----------------------------------------------------------------------- |
|
|
58 C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE |
|
|
59 C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE |
|
|
60 C EXPANSION FOR THE IMAGINARY PART. |
|
|
61 C----------------------------------------------------------------------- |
|
|
62 AEZ = 8.0D0*AZ |
|
|
63 S = TOL/AEZ |
|
|
64 JL = INT(SNGL(RL+RL)) + 2 |
|
|
65 P1R = ZEROR |
|
|
66 P1I = ZEROI |
|
|
67 IF (ZI.EQ.0.0D0) GO TO 30 |
|
|
68 C----------------------------------------------------------------------- |
|
|
69 C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF |
|
|
70 C SIGNIFICANCE WHEN FNU OR N IS LARGE |
|
|
71 C----------------------------------------------------------------------- |
|
|
72 INU = INT(SNGL(FNU)) |
|
|
73 ARG = (FNU-DBLE(FLOAT(INU)))*PI |
|
|
74 INU = INU + N - IL |
|
|
75 AK = -DSIN(ARG) |
|
|
76 BK = DCOS(ARG) |
|
|
77 IF (ZI.LT.0.0D0) BK = -BK |
|
|
78 P1R = AK |
|
|
79 P1I = BK |
|
|
80 IF (MOD(INU,2).EQ.0) GO TO 30 |
|
|
81 P1R = -P1R |
|
|
82 P1I = -P1I |
|
|
83 30 CONTINUE |
|
|
84 DO 70 K=1,IL |
|
|
85 SQK = FDN - 1.0D0 |
|
|
86 ATOL = S*DABS(SQK) |
|
|
87 SGN = 1.0D0 |
|
|
88 CS1R = CONER |
|
|
89 CS1I = CONEI |
|
|
90 CS2R = CONER |
|
|
91 CS2I = CONEI |
|
|
92 CKR = CONER |
|
|
93 CKI = CONEI |
|
|
94 AK = 0.0D0 |
|
|
95 AA = 1.0D0 |
|
|
96 BB = AEZ |
|
|
97 DKR = EZR |
|
|
98 DKI = EZI |
|
|
99 DO 40 J=1,JL |
|
|
100 CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI) |
|
|
101 CKR = STR*SQK |
|
|
102 CKI = STI*SQK |
|
|
103 CS2R = CS2R + CKR |
|
|
104 CS2I = CS2I + CKI |
|
|
105 SGN = -SGN |
|
|
106 CS1R = CS1R + CKR*SGN |
|
|
107 CS1I = CS1I + CKI*SGN |
|
|
108 DKR = DKR + EZR |
|
|
109 DKI = DKI + EZI |
|
|
110 AA = AA*DABS(SQK)/BB |
|
|
111 BB = BB + AEZ |
|
|
112 AK = AK + 8.0D0 |
|
|
113 SQK = SQK - AK |
|
|
114 IF (AA.LE.ATOL) GO TO 50 |
|
|
115 40 CONTINUE |
|
|
116 GO TO 110 |
|
|
117 50 CONTINUE |
|
|
118 S2R = CS1R |
|
|
119 S2I = CS1I |
|
|
120 IF (ZR+ZR.GE.ELIM) GO TO 60 |
|
|
121 TZR = ZR + ZR |
|
|
122 TZI = ZI + ZI |
|
|
123 CALL XZEXP(-TZR, -TZI, STR, STI) |
|
|
124 CALL ZMLT(STR, STI, P1R, P1I, STR, STI) |
|
|
125 CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI) |
|
|
126 S2R = S2R + STR |
|
|
127 S2I = S2I + STI |
|
|
128 60 CONTINUE |
|
|
129 FDN = FDN + 8.0D0*DFNU + 4.0D0 |
|
|
130 P1R = -P1R |
|
|
131 P1I = -P1I |
|
|
132 M = N - IL + K |
|
|
133 YR(M) = S2R*AK1R - S2I*AK1I |
|
|
134 YI(M) = S2R*AK1I + S2I*AK1R |
|
|
135 70 CONTINUE |
|
|
136 IF (N.LE.2) RETURN |
|
|
137 NN = N |
|
|
138 K = NN - 2 |
|
|
139 AK = DBLE(FLOAT(K)) |
|
|
140 STR = ZR*RAZ |
|
|
141 STI = -ZI*RAZ |
|
|
142 RZR = (STR+STR)*RAZ |
|
|
143 RZI = (STI+STI)*RAZ |
|
|
144 IB = 3 |
|
|
145 DO 80 I=IB,NN |
|
|
146 YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2) |
|
|
147 YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2) |
|
|
148 AK = AK - 1.0D0 |
|
|
149 K = K - 1 |
|
|
150 80 CONTINUE |
|
|
151 IF (KODED.EQ.0) RETURN |
|
|
152 CALL XZEXP(CZR, CZI, CKR, CKI) |
|
|
153 DO 90 I=1,NN |
|
|
154 STR = YR(I)*CKR - YI(I)*CKI |
|
|
155 YI(I) = YR(I)*CKI + YI(I)*CKR |
|
|
156 YR(I) = STR |
|
|
157 90 CONTINUE |
|
|
158 RETURN |
|
|
159 100 CONTINUE |
|
|
160 NZ = -1 |
|
|
161 RETURN |
|
|
162 110 CONTINUE |
|
|
163 NZ=-2 |
|
|
164 RETURN |
|
|
165 END |
