|
3217
|
1 SUBROUTINE ZACON(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL, |
|
|
2 * TOL, ELIM, ALIM) |
|
|
3 C***BEGIN PROLOGUE ZACON |
|
|
4 C***REFER TO ZBESK,ZBESH |
|
|
5 C |
|
|
6 C ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA |
|
|
7 C |
|
|
8 C K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN) |
|
|
9 C MP=PI*MR*CMPLX(0.0,1.0) |
|
|
10 C |
|
|
11 C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT |
|
|
12 C HALF Z PLANE |
|
|
13 C |
|
|
14 C***ROUTINES CALLED ZBINU,ZBKNU,ZS1S2,D1MACH,XZABS,ZMLT |
|
|
15 C***END PROLOGUE ZACON |
|
|
16 C COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST, |
|
|
17 C *S1,S2,Y,Z,ZN |
|
|
18 DOUBLE PRECISION ALIM, ARG, ASCLE, AS2, AZN, BRY, BSCLE, CKI, |
|
|
19 * CKR, CONER, CPN, CSCL, CSCR, CSGNI, CSGNR, CSPNI, CSPNR, |
|
|
20 * CSR, CSRR, CSSR, CYI, CYR, C1I, C1M, C1R, C2I, C2R, ELIM, FMR, |
|
|
21 * FN, FNU, FNUL, PI, PTI, PTR, RAZN, RL, RZI, RZR, SC1I, SC1R, |
|
|
22 * SC2I, SC2R, SGN, SPN, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, |
|
|
23 * YY, ZEROR, ZI, ZNI, ZNR, ZR, D1MACH, XZABS |
|
|
24 INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ |
|
|
25 DIMENSION YR(N), YI(N), CYR(2), CYI(2), CSSR(3), CSRR(3), BRY(3) |
|
|
26 DATA PI / 3.14159265358979324D0 / |
|
|
27 DATA ZEROR,CONER / 0.0D0,1.0D0 / |
|
|
28 NZ = 0 |
|
|
29 ZNR = -ZR |
|
|
30 ZNI = -ZI |
|
|
31 NN = N |
|
|
32 CALL ZBINU(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, FNUL, TOL, |
|
|
33 * ELIM, ALIM) |
|
|
34 IF (NW.LT.0) GO TO 90 |
|
|
35 C----------------------------------------------------------------------- |
|
|
36 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION |
|
|
37 C----------------------------------------------------------------------- |
|
|
38 NN = MIN0(2,N) |
|
|
39 CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM) |
|
|
40 IF (NW.NE.0) GO TO 90 |
|
|
41 S1R = CYR(1) |
|
|
42 S1I = CYI(1) |
|
|
43 FMR = DBLE(FLOAT(MR)) |
|
|
44 SGN = -DSIGN(PI,FMR) |
|
|
45 CSGNR = ZEROR |
|
|
46 CSGNI = SGN |
|
|
47 IF (KODE.EQ.1) GO TO 10 |
|
|
48 YY = -ZNI |
|
|
49 CPN = DCOS(YY) |
|
|
50 SPN = DSIN(YY) |
|
|
51 CALL ZMLT(CSGNR, CSGNI, CPN, SPN, CSGNR, CSGNI) |
|
|
52 10 CONTINUE |
|
|
53 C----------------------------------------------------------------------- |
|
|
54 C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE |
|
|
55 C WHEN FNU IS LARGE |
|
|
56 C----------------------------------------------------------------------- |
|
|
57 INU = INT(SNGL(FNU)) |
|
|
58 ARG = (FNU-DBLE(FLOAT(INU)))*SGN |
|
|
59 CPN = DCOS(ARG) |
|
|
60 SPN = DSIN(ARG) |
|
|
61 CSPNR = CPN |
|
|
62 CSPNI = SPN |
|
|
63 IF (MOD(INU,2).EQ.0) GO TO 20 |
|
|
64 CSPNR = -CSPNR |
|
|
65 CSPNI = -CSPNI |
|
|
66 20 CONTINUE |
|
|
67 IUF = 0 |
|
|
68 C1R = S1R |
|
|
69 C1I = S1I |
|
|
70 C2R = YR(1) |
|
|
71 C2I = YI(1) |
|
|
72 ASCLE = 1.0D+3*D1MACH(1)/TOL |
|
|
73 IF (KODE.EQ.1) GO TO 30 |
|
|
74 CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF) |
|
|
75 NZ = NZ + NW |
|
|
76 SC1R = C1R |
|
|
77 SC1I = C1I |
|
|
78 30 CONTINUE |
|
|
79 CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI) |
|
|
80 CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI) |
|
|
81 YR(1) = STR + PTR |
|
|
82 YI(1) = STI + PTI |
|
|
83 IF (N.EQ.1) RETURN |
|
|
84 CSPNR = -CSPNR |
|
|
85 CSPNI = -CSPNI |
|
|
86 S2R = CYR(2) |
|
|
87 S2I = CYI(2) |
|
|
88 C1R = S2R |
|
|
89 C1I = S2I |
|
|
90 C2R = YR(2) |
|
|
91 C2I = YI(2) |
|
|
92 IF (KODE.EQ.1) GO TO 40 |
|
|
93 CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF) |
|
|
94 NZ = NZ + NW |
|
|
95 SC2R = C1R |
|
|
96 SC2I = C1I |
|
|
97 40 CONTINUE |
|
|
98 CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI) |
|
|
99 CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI) |
|
|
100 YR(2) = STR + PTR |
|
|
101 YI(2) = STI + PTI |
|
|
102 IF (N.EQ.2) RETURN |
|
|
103 CSPNR = -CSPNR |
|
|
104 CSPNI = -CSPNI |
|
|
105 AZN = XZABS(ZNR,ZNI) |
|
|
106 RAZN = 1.0D0/AZN |
|
|
107 STR = ZNR*RAZN |
|
|
108 STI = -ZNI*RAZN |
|
|
109 RZR = (STR+STR)*RAZN |
|
|
110 RZI = (STI+STI)*RAZN |
|
|
111 FN = FNU + 1.0D0 |
|
|
112 CKR = FN*RZR |
|
|
113 CKI = FN*RZI |
|
|
114 C----------------------------------------------------------------------- |
|
|
115 C SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS |
|
|
116 C----------------------------------------------------------------------- |
|
|
117 CSCL = 1.0D0/TOL |
|
|
118 CSCR = TOL |
|
|
119 CSSR(1) = CSCL |
|
|
120 CSSR(2) = CONER |
|
|
121 CSSR(3) = CSCR |
|
|
122 CSRR(1) = CSCR |
|
|
123 CSRR(2) = CONER |
|
|
124 CSRR(3) = CSCL |
|
|
125 BRY(1) = ASCLE |
|
|
126 BRY(2) = 1.0D0/ASCLE |
|
|
127 BRY(3) = D1MACH(2) |
|
|
128 AS2 = XZABS(S2R,S2I) |
|
|
129 KFLAG = 2 |
|
|
130 IF (AS2.GT.BRY(1)) GO TO 50 |
|
|
131 KFLAG = 1 |
|
|
132 GO TO 60 |
|
|
133 50 CONTINUE |
|
|
134 IF (AS2.LT.BRY(2)) GO TO 60 |
|
|
135 KFLAG = 3 |
|
|
136 60 CONTINUE |
|
|
137 BSCLE = BRY(KFLAG) |
|
|
138 S1R = S1R*CSSR(KFLAG) |
|
|
139 S1I = S1I*CSSR(KFLAG) |
|
|
140 S2R = S2R*CSSR(KFLAG) |
|
|
141 S2I = S2I*CSSR(KFLAG) |
|
|
142 CSR = CSRR(KFLAG) |
|
|
143 DO 80 I=3,N |
|
|
144 STR = S2R |
|
|
145 STI = S2I |
|
|
146 S2R = CKR*STR - CKI*STI + S1R |
|
|
147 S2I = CKR*STI + CKI*STR + S1I |
|
|
148 S1R = STR |
|
|
149 S1I = STI |
|
|
150 C1R = S2R*CSR |
|
|
151 C1I = S2I*CSR |
|
|
152 STR = C1R |
|
|
153 STI = C1I |
|
|
154 C2R = YR(I) |
|
|
155 C2I = YI(I) |
|
|
156 IF (KODE.EQ.1) GO TO 70 |
|
|
157 IF (IUF.LT.0) GO TO 70 |
|
|
158 CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF) |
|
|
159 NZ = NZ + NW |
|
|
160 SC1R = SC2R |
|
|
161 SC1I = SC2I |
|
|
162 SC2R = C1R |
|
|
163 SC2I = C1I |
|
|
164 IF (IUF.NE.3) GO TO 70 |
|
|
165 IUF = -4 |
|
|
166 S1R = SC1R*CSSR(KFLAG) |
|
|
167 S1I = SC1I*CSSR(KFLAG) |
|
|
168 S2R = SC2R*CSSR(KFLAG) |
|
|
169 S2I = SC2I*CSSR(KFLAG) |
|
|
170 STR = SC2R |
|
|
171 STI = SC2I |
|
|
172 70 CONTINUE |
|
|
173 PTR = CSPNR*C1R - CSPNI*C1I |
|
|
174 PTI = CSPNR*C1I + CSPNI*C1R |
|
|
175 YR(I) = PTR + CSGNR*C2R - CSGNI*C2I |
|
|
176 YI(I) = PTI + CSGNR*C2I + CSGNI*C2R |
|
|
177 CKR = CKR + RZR |
|
|
178 CKI = CKI + RZI |
|
|
179 CSPNR = -CSPNR |
|
|
180 CSPNI = -CSPNI |
|
|
181 IF (KFLAG.GE.3) GO TO 80 |
|
|
182 PTR = DABS(C1R) |
|
|
183 PTI = DABS(C1I) |
|
|
184 C1M = DMAX1(PTR,PTI) |
|
|
185 IF (C1M.LE.BSCLE) GO TO 80 |
|
|
186 KFLAG = KFLAG + 1 |
|
|
187 BSCLE = BRY(KFLAG) |
|
|
188 S1R = S1R*CSR |
|
|
189 S1I = S1I*CSR |
|
|
190 S2R = STR |
|
|
191 S2I = STI |
|
|
192 S1R = S1R*CSSR(KFLAG) |
|
|
193 S1I = S1I*CSSR(KFLAG) |
|
|
194 S2R = S2R*CSSR(KFLAG) |
|
|
195 S2I = S2I*CSSR(KFLAG) |
|
|
196 CSR = CSRR(KFLAG) |
|
|
197 80 CONTINUE |
|
|
198 RETURN |
|
|
199 90 CONTINUE |
|
|
200 NZ = -1 |
|
|
201 IF(NW.EQ.(-2)) NZ=-2 |
|
|
202 RETURN |
|
|
203 END |
