GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
cunk1.f
Go to the documentation of this file.
1 SUBROUTINE cunk1(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
2C***BEGIN PROLOGUE CUNK1
3C***REFER TO CBESK
4C
5C CUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
6C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
7C UNIFORM ASYMPTOTIC EXPANSION.
8C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
9C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
10C
11C***ROUTINES CALLED CS1S2,CUCHK,CUNIK,R1MACH
12C***END PROLOGUE CUNK1
13 COMPLEX CFN, CK, CONE, CRSC, CS, CSCL, CSGN, CSPN, CSR, CSS,
14 * CWRK, CY, CZERO, C1, C2, PHI, RZ, SUM, S1, S2, Y, Z,
15 * ZETA1, ZETA2, ZR, PHID, ZETA1D, ZETA2D, SUMD
16 REAL ALIM, ANG, APHI, ASC, ASCLE, BRY, CPN, C2I, C2M, C2R, ELIM,
17 * FMR, FN, FNF, FNU, PI, RS1, SGN, SPN, TOL, X, R1MACH
18 INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
19 * KK, KODE, MR, N, NW, NZ, J, IPARD, INITD, IC
20 dimension bry(3), init(2), y(n), sum(2), phi(2), zeta1(2),
21 * zeta2(2), cy(2), cwrk(16,3), css(3), csr(3)
22 DATA czero, cone / (0.0e0,0.0e0) , (1.0e0,0.0e0) /
23 DATA pi / 3.14159265358979324e0 /
24C
25 kdflg = 1
26 nz = 0
27C-----------------------------------------------------------------------
28C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
29C THE UNDERFLOW LIMIT
30C-----------------------------------------------------------------------
31 cscl = cmplx(1.0e0/tol,0.0e0)
32 crsc = cmplx(tol,0.0e0)
33 css(1) = cscl
34 css(2) = cone
35 css(3) = crsc
36 csr(1) = crsc
37 csr(2) = cone
38 csr(3) = cscl
39 bry(1) = 1.0e+3*r1mach(1)/tol
40 bry(2) = 1.0e0/bry(1)
41 bry(3) = r1mach(2)
42 x = real(z)
43 zr = z
44 IF (x.LT.0.0e0) zr = -z
45 j=2
46 DO 70 i=1,n
47C-----------------------------------------------------------------------
48C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
49C-----------------------------------------------------------------------
50 j = 3 - j
51 fn = fnu + float(i-1)
52 init(j) = 0
53 CALL cunik(zr, fn, 2, 0, tol, init(j), phi(j), zeta1(j),
54 * zeta2(j), sum(j), cwrk(1,j))
55 IF (kode.EQ.1) GO TO 20
56 cfn = cmplx(fn,0.0e0)
57 s1 = zeta1(j) - cfn*(cfn/(zr+zeta2(j)))
58 GO TO 30
59 20 CONTINUE
60 s1 = zeta1(j) - zeta2(j)
61 30 CONTINUE
62C-----------------------------------------------------------------------
63C TEST FOR UNDERFLOW AND OVERFLOW
64C-----------------------------------------------------------------------
65 rs1 = real(s1)
66 IF (abs(rs1).GT.elim) GO TO 60
67 IF (kdflg.EQ.1) kflag = 2
68 IF (abs(rs1).LT.alim) GO TO 40
69C-----------------------------------------------------------------------
70C REFINE TEST AND SCALE
71C-----------------------------------------------------------------------
72 aphi = cabs(phi(j))
73 rs1 = rs1 + alog(aphi)
74 IF (abs(rs1).GT.elim) GO TO 60
75 IF (kdflg.EQ.1) kflag = 1
76 IF (rs1.LT.0.0e0) GO TO 40
77 IF (kdflg.EQ.1) kflag = 3
78 40 CONTINUE
79C-----------------------------------------------------------------------
80C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
81C EXPONENT EXTREMES
82C-----------------------------------------------------------------------
83 s2 = phi(j)*sum(j)
84 c2r = real(s1)
85 c2i = aimag(s1)
86 c2m = exp(c2r)*real(css(kflag))
87 s1 = cmplx(c2m,0.0e0)*cmplx(cos(c2i),sin(c2i))
88 s2 = s2*s1
89 IF (kflag.NE.1) GO TO 50
90 CALL cuchk(s2, nw, bry(1), tol)
91 IF (nw.NE.0) GO TO 60
92 50 CONTINUE
93 cy(kdflg) = s2
94 y(i) = s2*csr(kflag)
95 IF (kdflg.EQ.2) GO TO 75
96 kdflg = 2
97 GO TO 70
98 60 CONTINUE
99 IF (rs1.GT.0.0e0) GO TO 290
100C-----------------------------------------------------------------------
101C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
102C-----------------------------------------------------------------------
103 IF (x.LT.0.0e0) GO TO 290
104 kdflg = 1
105 y(i) = czero
106 nz=nz+1
107 IF (i.EQ.1) GO TO 70
108 IF (y(i-1).EQ.czero) GO TO 70
109 y(i-1) = czero
110 nz=nz+1
111 70 CONTINUE
112 i=n
113 75 CONTINUE
114 rz = cmplx(2.0e0,0.0e0)/zr
115 ck = cmplx(fn,0.0e0)*rz
116 ib = i+1
117 IF (n.LT.ib) GO TO 160
118C-----------------------------------------------------------------------
119C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
120C ON UNDERFLOW
121C-----------------------------------------------------------------------
122 fn = fnu+float(n-1)
123 ipard = 1
124 IF (mr.NE.0) ipard = 0
125 initd = 0
126 CALL cunik(zr,fn,2,ipard,tol,initd,phid,zeta1d,zeta2d,sumd,
127 *cwrk(1,3))
128 IF (kode.EQ.1) GO TO 80
129 cfn=cmplx(fn,0.0e0)
130 s1=zeta1d-cfn*(cfn/(zr+zeta2d))
131 GO TO 90
132 80 CONTINUE
133 s1=zeta1d-zeta2d
134 90 CONTINUE
135 rs1=real(s1)
136 IF (abs(rs1).GT.elim) GO TO 95
137 IF (abs(rs1).LT.alim) GO TO 100
138C-----------------------------------------------------------------------
139C REFINE ESTIMATE AND TEST
140C-----------------------------------------------------------------------
141 aphi=cabs(phid)
142 rs1=rs1+alog(aphi)
143 IF (abs(rs1).LT.elim) GO TO 100
144 95 CONTINUE
145 IF (rs1.GT.0.0e0) GO TO 290
146C-----------------------------------------------------------------------
147C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
148C-----------------------------------------------------------------------
149 IF (x.LT.0.0e0) GO TO 290
150 nz=n
151 DO 96 i=1,n
152 y(i) = czero
153 96 CONTINUE
154 RETURN
155 100 CONTINUE
156C-----------------------------------------------------------------------
157C RECUR FORWARD FOR REMAINDER OF THE SEQUENCE
158C-----------------------------------------------------------------------
159 s1 = cy(1)
160 s2 = cy(2)
161 c1 = csr(kflag)
162 ascle = bry(kflag)
163 DO 120 i=ib,n
164 c2 = s2
165 s2 = ck*s2 + s1
166 s1 = c2
167 ck = ck + rz
168 c2 = s2*c1
169 y(i) = c2
170 IF (kflag.GE.3) GO TO 120
171 c2r = real(c2)
172 c2i = aimag(c2)
173 c2r = abs(c2r)
174 c2i = abs(c2i)
175 c2m = amax1(c2r,c2i)
176 IF (c2m.LE.ascle) GO TO 120
177 kflag = kflag + 1
178 ascle = bry(kflag)
179 s1 = s1*c1
180 s2 = c2
181 s1 = s1*css(kflag)
182 s2 = s2*css(kflag)
183 c1 = csr(kflag)
184 120 CONTINUE
185 160 CONTINUE
186 IF (mr.EQ.0) RETURN
187C-----------------------------------------------------------------------
188C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0
189C-----------------------------------------------------------------------
190 nz = 0
191 fmr = float(mr)
192 sgn = -sign(pi,fmr)
193C-----------------------------------------------------------------------
194C CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
195C-----------------------------------------------------------------------
196 csgn = cmplx(0.0e0,sgn)
197 inu = int(fnu)
198 fnf = fnu - float(inu)
199 ifn = inu + n - 1
200 ang = fnf*sgn
201 cpn = cos(ang)
202 spn = sin(ang)
203 cspn = cmplx(cpn,spn)
204 IF (mod(ifn,2).EQ.1) cspn = -cspn
205 asc = bry(1)
206 kk = n
207 iuf = 0
208 kdflg = 1
209 ib = ib-1
210 ic = ib-1
211 DO 260 k=1,n
212 fn = fnu + float(kk-1)
213C-----------------------------------------------------------------------
214C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
215C FUNCTION ABOVE
216C-----------------------------------------------------------------------
217 m=3
218 IF (n.GT.2) GO TO 175
219 170 CONTINUE
220 initd = init(j)
221 phid = phi(j)
222 zeta1d = zeta1(j)
223 zeta2d = zeta2(j)
224 sumd = sum(j)
225 m = j
226 j = 3 - j
227 GO TO 180
228 175 CONTINUE
229 IF ((kk.EQ.n).AND.(ib.LT.n)) GO TO 180
230 IF ((kk.EQ.ib).OR.(kk.EQ.ic)) GO TO 170
231 initd = 0
232 180 CONTINUE
233 CALL cunik(zr, fn, 1, 0, tol, initd, phid, zeta1d,
234 * zeta2d, sumd, cwrk(1,m))
235 IF (kode.EQ.1) GO TO 190
236 cfn = cmplx(fn,0.0e0)
237 s1 = -zeta1d + cfn*(cfn/(zr+zeta2d))
238 GO TO 200
239 190 CONTINUE
240 s1 = -zeta1d + zeta2d
241 200 CONTINUE
242C-----------------------------------------------------------------------
243C TEST FOR UNDERFLOW AND OVERFLOW
244C-----------------------------------------------------------------------
245 rs1 = real(s1)
246 IF (abs(rs1).GT.elim) GO TO 250
247 IF (kdflg.EQ.1) iflag = 2
248 IF (abs(rs1).LT.alim) GO TO 210
249C-----------------------------------------------------------------------
250C REFINE TEST AND SCALE
251C-----------------------------------------------------------------------
252 aphi = cabs(phid)
253 rs1 = rs1 + alog(aphi)
254 IF (abs(rs1).GT.elim) GO TO 250
255 IF (kdflg.EQ.1) iflag = 1
256 IF (rs1.LT.0.0e0) GO TO 210
257 IF (kdflg.EQ.1) iflag = 3
258 210 CONTINUE
259 s2 = csgn*phid*sumd
260 c2r = real(s1)
261 c2i = aimag(s1)
262 c2m = exp(c2r)*real(css(iflag))
263 s1 = cmplx(c2m,0.0e0)*cmplx(cos(c2i),sin(c2i))
264 s2 = s2*s1
265 IF (iflag.NE.1) GO TO 220
266 CALL cuchk(s2, nw, bry(1), tol)
267 IF (nw.NE.0) s2 = cmplx(0.0e0,0.0e0)
268 220 CONTINUE
269 cy(kdflg) = s2
270 c2 = s2
271 s2 = s2*csr(iflag)
272C-----------------------------------------------------------------------
273C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
274C-----------------------------------------------------------------------
275 s1 = y(kk)
276 IF (kode.EQ.1) GO TO 240
277 CALL cs1s2(zr, s1, s2, nw, asc, alim, iuf)
278 nz = nz + nw
279 240 CONTINUE
280 y(kk) = s1*cspn + s2
281 kk = kk - 1
282 cspn = -cspn
283 IF (c2.NE.czero) GO TO 245
284 kdflg = 1
285 GO TO 260
286 245 CONTINUE
287 IF (kdflg.EQ.2) GO TO 265
288 kdflg = 2
289 GO TO 260
290 250 CONTINUE
291 IF (rs1.GT.0.0e0) GO TO 290
292 s2 = czero
293 GO TO 220
294 260 CONTINUE
295 k = n
296 265 CONTINUE
297 il = n - k
298 IF (il.EQ.0) RETURN
299C-----------------------------------------------------------------------
300C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
301C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
302C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
303C-----------------------------------------------------------------------
304 s1 = cy(1)
305 s2 = cy(2)
306 cs = csr(iflag)
307 ascle = bry(iflag)
308 fn = float(inu+il)
309 DO 280 i=1,il
310 c2 = s2
311 s2 = s1 + cmplx(fn+fnf,0.0e0)*rz*s2
312 s1 = c2
313 fn = fn - 1.0e0
314 c2 = s2*cs
315 ck = c2
316 c1 = y(kk)
317 IF (kode.EQ.1) GO TO 270
318 CALL cs1s2(zr, c1, c2, nw, asc, alim, iuf)
319 nz = nz + nw
320 270 CONTINUE
321 y(kk) = c1*cspn + c2
322 kk = kk - 1
323 cspn = -cspn
324 IF (iflag.GE.3) GO TO 280
325 c2r = real(ck)
326 c2i = aimag(ck)
327 c2r = abs(c2r)
328 c2i = abs(c2i)
329 c2m = amax1(c2r,c2i)
330 IF (c2m.LE.ascle) GO TO 280
331 iflag = iflag + 1
332 ascle = bry(iflag)
333 s1 = s1*cs
334 s2 = ck
335 s1 = s1*css(iflag)
336 s2 = s2*css(iflag)
337 cs = csr(iflag)
338 280 CONTINUE
339 RETURN
340 290 CONTINUE
341 nz = -1
342 RETURN
343 END
cmplx
double complex cmplx
Definition Faddeeva.cc:227
cs1s2
subroutine cs1s2(zr, s1, s2, nz, ascle, alim, iuf)
Definition cs1s2.f:2
cuchk
subroutine cuchk(y, nz, ascle, tol)
Definition cuchk.f:2
cunik
subroutine cunik(zr, fnu, ikflg, ipmtr, tol, init, phi, zeta1, zeta2, sum, cwrk)
Definition cunik.f:3
cunk1
subroutine cunk1(z, fnu, kode, mr, n, y, nz, tol, elim, alim)
Definition cunk1.f:2
real
ColumnVector real(const ComplexColumnVector &a)
Definition dColVector.cc:137
mod
T mod(T x, T y)
Definition lo-mappers.h:294