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
|
SUBROUTINE ZUOIK(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL,
* ELIM, ALIM)
C***BEGIN PROLOGUE ZUOIK
C***REFER TO ZBESI,ZBESK,ZBESH
C
C ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
C EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
C (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
C WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
C EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
C THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
C MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
C EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
C EXP(-ELIM)/TOL
C
C IKFLG=1 MEANS THE I SEQUENCE IS TESTED
C =2 MEANS THE K SEQUENCE IS TESTED
C NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
C =-1 MEANS AN OVERFLOW WOULD OCCUR
C IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
C THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
C IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
C IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
C ANOTHER ROUTINE
C
C***ROUTINES CALLED ZUCHK,ZUNHJ,ZUNIK,D1MACH,AZABS,AZLOG
C***END PROLOGUE ZUOIK
C COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN,
C *ZR
DOUBLE PRECISION AARG, AIC, ALIM, APHI, ARGI, ARGR, ASUMI, ASUMR,
* ASCLE, AX, AY, BSUMI, BSUMR, CWRKI, CWRKR, CZI, CZR, ELIM, FNN,
* FNU, GNN, GNU, PHII, PHIR, RCZ, STR, STI, SUMI, SUMR, TOL, YI,
* YR, ZBI, ZBR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI,
* ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, AZABS
INTEGER I, IDUM, IFORM, IKFLG, INIT, KODE, N, NN, NUF, NW
DIMENSION YR(N), YI(N), CWRKR(16), CWRKI(16)
DATA ZEROR,ZEROI / 0.0D0, 0.0D0 /
DATA AIC / 1.265512123484645396D+00 /
NUF = 0
NN = N
ZRR = ZR
ZRI = ZI
IF (ZR.GE.0.0D0) GO TO 10
ZRR = -ZR
ZRI = -ZI
10 CONTINUE
ZBR = ZRR
ZBI = ZRI
AX = DABS(ZR)*1.7321D0
AY = DABS(ZI)
IFORM = 1
IF (AY.GT.AX) IFORM = 2
GNU = DMAX1(FNU,1.0D0)
IF (IKFLG.EQ.1) GO TO 20
FNN = DBLE(FLOAT(NN))
GNN = FNU + FNN - 1.0D0
GNU = DMAX1(GNN,FNN)
20 CONTINUE
C-----------------------------------------------------------------------
C ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
C REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
C THE SIGN OF THE IMAGINARY PART CORRECT.
C-----------------------------------------------------------------------
IF (IFORM.EQ.2) GO TO 30
INIT = 0
CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII,
* ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
CZR = -ZETA1R + ZETA2R
CZI = -ZETA1I + ZETA2I
GO TO 50
30 CONTINUE
ZNR = ZRI
ZNI = -ZRR
IF (ZI.GT.0.0D0) GO TO 40
ZNR = -ZNR
40 CONTINUE
CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
* ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
CZR = -ZETA1R + ZETA2R
CZI = -ZETA1I + ZETA2I
AARG = AZABS(ARGR,ARGI)
50 CONTINUE
IF (KODE.EQ.1) GO TO 60
CZR = CZR - ZBR
CZI = CZI - ZBI
60 CONTINUE
IF (IKFLG.EQ.1) GO TO 70
CZR = -CZR
CZI = -CZI
70 CONTINUE
APHI = AZABS(PHIR,PHII)
RCZ = CZR
C-----------------------------------------------------------------------
C OVERFLOW TEST
C-----------------------------------------------------------------------
IF (RCZ.GT.ELIM) GO TO 210
IF (RCZ.LT.ALIM) GO TO 80
RCZ = RCZ + DLOG(APHI)
IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*DLOG(AARG) - AIC
IF (RCZ.GT.ELIM) GO TO 210
GO TO 130
80 CONTINUE
C-----------------------------------------------------------------------
C UNDERFLOW TEST
C-----------------------------------------------------------------------
IF (RCZ.LT.(-ELIM)) GO TO 90
IF (RCZ.GT.(-ALIM)) GO TO 130
RCZ = RCZ + DLOG(APHI)
IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*DLOG(AARG) - AIC
IF (RCZ.GT.(-ELIM)) GO TO 110
90 CONTINUE
DO 100 I=1,NN
YR(I) = ZEROR
YI(I) = ZEROI
100 CONTINUE
NUF = NN
RETURN
110 CONTINUE
ASCLE = 1.0D+3*D1MACH(1)/TOL
CALL AZLOG(PHIR, PHII, STR, STI, IDUM)
CZR = CZR + STR
CZI = CZI + STI
IF (IFORM.EQ.1) GO TO 120
CALL AZLOG(ARGR, ARGI, STR, STI, IDUM)
CZR = CZR - 0.25D0*STR - AIC
CZI = CZI - 0.25D0*STI
120 CONTINUE
AX = DEXP(RCZ)/TOL
AY = CZI
CZR = AX*DCOS(AY)
CZI = AX*DSIN(AY)
CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL)
IF (NW.NE.0) GO TO 90
130 CONTINUE
IF (IKFLG.EQ.2) RETURN
IF (N.EQ.1) RETURN
C-----------------------------------------------------------------------
C SET UNDERFLOWS ON I SEQUENCE
C-----------------------------------------------------------------------
140 CONTINUE
GNU = FNU + DBLE(FLOAT(NN-1))
IF (IFORM.EQ.2) GO TO 150
INIT = 0
CALL ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, PHIR, PHII,
* ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
CZR = -ZETA1R + ZETA2R
CZI = -ZETA1I + ZETA2I
GO TO 160
150 CONTINUE
CALL ZUNHJ(ZNR, ZNI, GNU, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
* ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
CZR = -ZETA1R + ZETA2R
CZI = -ZETA1I + ZETA2I
AARG = AZABS(ARGR,ARGI)
160 CONTINUE
IF (KODE.EQ.1) GO TO 170
CZR = CZR - ZBR
CZI = CZI - ZBI
170 CONTINUE
APHI = AZABS(PHIR,PHII)
RCZ = CZR
IF (RCZ.LT.(-ELIM)) GO TO 180
IF (RCZ.GT.(-ALIM)) RETURN
RCZ = RCZ + DLOG(APHI)
IF (IFORM.EQ.2) RCZ = RCZ - 0.25D0*DLOG(AARG) - AIC
IF (RCZ.GT.(-ELIM)) GO TO 190
180 CONTINUE
YR(NN) = ZEROR
YI(NN) = ZEROI
NN = NN - 1
NUF = NUF + 1
IF (NN.EQ.0) RETURN
GO TO 140
190 CONTINUE
ASCLE = 1.0D+3*D1MACH(1)/TOL
CALL AZLOG(PHIR, PHII, STR, STI, IDUM)
CZR = CZR + STR
CZI = CZI + STI
IF (IFORM.EQ.1) GO TO 200
CALL AZLOG(ARGR, ARGI, STR, STI, IDUM)
CZR = CZR - 0.25D0*STR - AIC
CZI = CZI - 0.25D0*STI
200 CONTINUE
AX = DEXP(RCZ)/TOL
AY = CZI
CZR = AX*DCOS(AY)
CZI = AX*DSIN(AY)
CALL ZUCHK(CZR, CZI, NW, ASCLE, TOL)
IF (NW.NE.0) GO TO 180
RETURN
210 CONTINUE
NUF = -1
RETURN
END
|
