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
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
|
SUBROUTINE ZACON(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL,
* TOL, ELIM, ALIM)
C***BEGIN PROLOGUE ZACON
C***REFER TO ZBESK,ZBESH
C
C ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA
C
C K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
C MP=PI*MR*CMPLX(0.0,1.0)
C
C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
C HALF Z PLANE
C
C***ROUTINES CALLED ZBINU,ZBKNU,ZS1S2,D1MACH,AZABS,ZMLT
C***END PROLOGUE ZACON
C COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST,
C *S1,S2,Y,Z,ZN
DOUBLE PRECISION ALIM, ARG, ASCLE, AS2, AZN, BRY, BSCLE, CKI,
* CKR, CONER, CPN, CSCL, CSCR, CSGNI, CSGNR, CSPNI, CSPNR,
* CSR, CSRR, CSSR, CYI, CYR, C1I, C1M, C1R, C2I, C2R, ELIM, FMR,
* FN, FNU, FNUL, PI, PTI, PTR, RAZN, RL, RZI, RZR, SC1I, SC1R,
* SC2I, SC2R, SGN, SPN, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR,
* YY, ZEROR, ZI, ZNI, ZNR, ZR, D1MACH, AZABS
INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ
DIMENSION YR(N), YI(N), CYR(2), CYI(2), CSSR(3), CSRR(3), BRY(3)
DATA PI / 3.14159265358979324D0 /
DATA ZEROR,CONER / 0.0D0,1.0D0 /
NZ = 0
ZNR = -ZR
ZNI = -ZI
NN = N
CALL ZBINU(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, FNUL, TOL,
* ELIM, ALIM)
IF (NW.LT.0) GO TO 90
C-----------------------------------------------------------------------
C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
C-----------------------------------------------------------------------
NN = MIN0(2,N)
CALL ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM)
IF (NW.NE.0) GO TO 90
S1R = CYR(1)
S1I = CYI(1)
FMR = DBLE(FLOAT(MR))
SGN = -DSIGN(PI,FMR)
CSGNR = ZEROR
CSGNI = SGN
IF (KODE.EQ.1) GO TO 10
YY = -ZNI
CPN = DCOS(YY)
SPN = DSIN(YY)
CALL ZMLT(CSGNR, CSGNI, CPN, SPN, CSGNR, CSGNI)
10 CONTINUE
C-----------------------------------------------------------------------
C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
C WHEN FNU IS LARGE
C-----------------------------------------------------------------------
INU = INT(SNGL(FNU))
ARG = (FNU-DBLE(FLOAT(INU)))*SGN
CPN = DCOS(ARG)
SPN = DSIN(ARG)
CSPNR = CPN
CSPNI = SPN
IF (MOD(INU,2).EQ.0) GO TO 20
CSPNR = -CSPNR
CSPNI = -CSPNI
20 CONTINUE
IUF = 0
C1R = S1R
C1I = S1I
C2R = YR(1)
C2I = YI(1)
ASCLE = 1.0D+3*D1MACH(1)/TOL
IF (KODE.EQ.1) GO TO 30
CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
NZ = NZ + NW
SC1R = C1R
SC1I = C1I
30 CONTINUE
CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI)
CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI)
YR(1) = STR + PTR
YI(1) = STI + PTI
IF (N.EQ.1) RETURN
CSPNR = -CSPNR
CSPNI = -CSPNI
S2R = CYR(2)
S2I = CYI(2)
C1R = S2R
C1I = S2I
C2R = YR(2)
C2I = YI(2)
IF (KODE.EQ.1) GO TO 40
CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
NZ = NZ + NW
SC2R = C1R
SC2I = C1I
40 CONTINUE
CALL ZMLT(CSPNR, CSPNI, C1R, C1I, STR, STI)
CALL ZMLT(CSGNR, CSGNI, C2R, C2I, PTR, PTI)
YR(2) = STR + PTR
YI(2) = STI + PTI
IF (N.EQ.2) RETURN
CSPNR = -CSPNR
CSPNI = -CSPNI
AZN = AZABS(ZNR,ZNI)
RAZN = 1.0D0/AZN
STR = ZNR*RAZN
STI = -ZNI*RAZN
RZR = (STR+STR)*RAZN
RZI = (STI+STI)*RAZN
FN = FNU + 1.0D0
CKR = FN*RZR
CKI = FN*RZI
C-----------------------------------------------------------------------
C SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
C-----------------------------------------------------------------------
CSCL = 1.0D0/TOL
CSCR = TOL
CSSR(1) = CSCL
CSSR(2) = CONER
CSSR(3) = CSCR
CSRR(1) = CSCR
CSRR(2) = CONER
CSRR(3) = CSCL
BRY(1) = ASCLE
BRY(2) = 1.0D0/ASCLE
BRY(3) = D1MACH(2)
AS2 = AZABS(S2R,S2I)
KFLAG = 2
IF (AS2.GT.BRY(1)) GO TO 50
KFLAG = 1
GO TO 60
50 CONTINUE
IF (AS2.LT.BRY(2)) GO TO 60
KFLAG = 3
60 CONTINUE
BSCLE = BRY(KFLAG)
S1R = S1R*CSSR(KFLAG)
S1I = S1I*CSSR(KFLAG)
S2R = S2R*CSSR(KFLAG)
S2I = S2I*CSSR(KFLAG)
CSR = CSRR(KFLAG)
DO 80 I=3,N
STR = S2R
STI = S2I
S2R = CKR*STR - CKI*STI + S1R
S2I = CKR*STI + CKI*STR + S1I
S1R = STR
S1I = STI
C1R = S2R*CSR
C1I = S2I*CSR
STR = C1R
STI = C1I
C2R = YR(I)
C2I = YI(I)
IF (KODE.EQ.1) GO TO 70
IF (IUF.LT.0) GO TO 70
CALL ZS1S2(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF)
NZ = NZ + NW
SC1R = SC2R
SC1I = SC2I
SC2R = C1R
SC2I = C1I
IF (IUF.NE.3) GO TO 70
IUF = -4
S1R = SC1R*CSSR(KFLAG)
S1I = SC1I*CSSR(KFLAG)
S2R = SC2R*CSSR(KFLAG)
S2I = SC2I*CSSR(KFLAG)
STR = SC2R
STI = SC2I
70 CONTINUE
PTR = CSPNR*C1R - CSPNI*C1I
PTI = CSPNR*C1I + CSPNI*C1R
YR(I) = PTR + CSGNR*C2R - CSGNI*C2I
YI(I) = PTI + CSGNR*C2I + CSGNI*C2R
CKR = CKR + RZR
CKI = CKI + RZI
CSPNR = -CSPNR
CSPNI = -CSPNI
IF (KFLAG.GE.3) GO TO 80
PTR = DABS(C1R)
PTI = DABS(C1I)
C1M = DMAX1(PTR,PTI)
IF (C1M.LE.BSCLE) GO TO 80
KFLAG = KFLAG + 1
BSCLE = BRY(KFLAG)
S1R = S1R*CSR
S1I = S1I*CSR
S2R = STR
S2I = STI
S1R = S1R*CSSR(KFLAG)
S1I = S1I*CSSR(KFLAG)
S2R = S2R*CSSR(KFLAG)
S2I = S2I*CSSR(KFLAG)
CSR = CSRR(KFLAG)
80 CONTINUE
RETURN
90 CONTINUE
NZ = -1
IF(NW.EQ.(-2)) NZ=-2
RETURN
END
|
