GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
zunk1.f
Go to the documentation of this file.
1 SUBROUTINE zunk1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
2 * ALIM)
3C***BEGIN PROLOGUE ZUNK1
4C***REFER TO ZBESK
5C
6C ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
7C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
8C UNIFORM ASYMPTOTIC EXPANSION.
9C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
10C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
11C
12C***ROUTINES CALLED ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,XZABS
13C***END PROLOGUE ZUNK1
14C COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
15C *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR
16 DOUBLE PRECISION ALIM, ANG, APHI, ASC, ASCLE, BRY, CKI, CKR,
17 * coner, crsc, cscl, csgni, cspni, cspnr, csr, csrr, cssr,
18 * cwrki, cwrkr, cyi, cyr, c1i, c1r, c2i, c2m, c2r, elim, fmr, fn,
19 * fnf, fnu, phidi, phidr, phii, phir, pi, rast, razr, rs1, rzi,
20 * rzr, sgn, sti, str, sumdi, sumdr, sumi, sumr, s1i, s1r, s2i,
21 * s2r, tol, yi, yr, zeroi, zeror, zeta1i, zeta1r, zeta2i, zeta2r,
22 * zet1di, zet1dr, zet2di, zet2dr, zi, zr, zri, zrr, d1mach, xzabs
23 INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
24 * kk, kode, mr, n, nw, nz, initd, ic, ipard, j
25 dimension bry(3), init(2), yr(n), yi(n), sumr(2), sumi(2),
26 * zeta1r(2), zeta1i(2), zeta2r(2), zeta2i(2), cyr(2), cyi(2),
27 * cwrkr(16,3), cwrki(16,3), cssr(3), csrr(3), phir(2), phii(2)
28 DATA zeror,zeroi,coner / 0.0d0, 0.0d0, 1.0d0 /
29 DATA pi / 3.14159265358979324d0 /
30C
31 kdflg = 1
32 nz = 0
33C-----------------------------------------------------------------------
34C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
35C THE UNDERFLOW LIMIT
36C-----------------------------------------------------------------------
37 cscl = 1.0d0/tol
38 crsc = tol
39 cssr(1) = cscl
40 cssr(2) = coner
41 cssr(3) = crsc
42 csrr(1) = crsc
43 csrr(2) = coner
44 csrr(3) = cscl
45 bry(1) = 1.0d+3*d1mach(1)/tol
46 bry(2) = 1.0d0/bry(1)
47 bry(3) = d1mach(2)
48 zrr = zr
49 zri = zi
50 IF (zr.GE.0.0d0) GO TO 10
51 zrr = -zr
52 zri = -zi
53 10 CONTINUE
54 j = 2
55 DO 70 i=1,n
56C-----------------------------------------------------------------------
57C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
58C-----------------------------------------------------------------------
59 j = 3 - j
60 fn = fnu + dble(float(i-1))
61 init(j) = 0
62 CALL zunik(zrr, zri, fn, 2, 0, tol, init(j), phir(j), phii(j),
63 * zeta1r(j), zeta1i(j), zeta2r(j), zeta2i(j), sumr(j), sumi(j),
64 * cwrkr(1,j), cwrki(1,j))
65 IF (kode.EQ.1) GO TO 20
66 str = zrr + zeta2r(j)
67 sti = zri + zeta2i(j)
68 rast = fn/xzabs(str,sti)
69 str = str*rast*rast
70 sti = -sti*rast*rast
71 s1r = zeta1r(j) - str
72 s1i = zeta1i(j) - sti
73 GO TO 30
74 20 CONTINUE
75 s1r = zeta1r(j) - zeta2r(j)
76 s1i = zeta1i(j) - zeta2i(j)
77 30 CONTINUE
78 rs1 = s1r
79C-----------------------------------------------------------------------
80C TEST FOR UNDERFLOW AND OVERFLOW
81C-----------------------------------------------------------------------
82 IF (dabs(rs1).GT.elim) GO TO 60
83 IF (kdflg.EQ.1) kflag = 2
84 IF (dabs(rs1).LT.alim) GO TO 40
85C-----------------------------------------------------------------------
86C REFINE TEST AND SCALE
87C-----------------------------------------------------------------------
88 aphi = xzabs(phir(j),phii(j))
89 rs1 = rs1 + dlog(aphi)
90 IF (dabs(rs1).GT.elim) GO TO 60
91 IF (kdflg.EQ.1) kflag = 1
92 IF (rs1.LT.0.0d0) GO TO 40
93 IF (kdflg.EQ.1) kflag = 3
94 40 CONTINUE
95C-----------------------------------------------------------------------
96C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
97C EXPONENT EXTREMES
98C-----------------------------------------------------------------------
99 s2r = phir(j)*sumr(j) - phii(j)*sumi(j)
100 s2i = phir(j)*sumi(j) + phii(j)*sumr(j)
101 str = dexp(s1r)*cssr(kflag)
102 s1r = str*dcos(s1i)
103 s1i = str*dsin(s1i)
104 str = s2r*s1r - s2i*s1i
105 s2i = s1r*s2i + s2r*s1i
106 s2r = str
107 IF (kflag.NE.1) GO TO 50
108 CALL zuchk(s2r, s2i, nw, bry(1), tol)
109 IF (nw.NE.0) GO TO 60
110 50 CONTINUE
111 cyr(kdflg) = s2r
112 cyi(kdflg) = s2i
113 yr(i) = s2r*csrr(kflag)
114 yi(i) = s2i*csrr(kflag)
115 IF (kdflg.EQ.2) GO TO 75
116 kdflg = 2
117 GO TO 70
118 60 CONTINUE
119 IF (rs1.GT.0.0d0) GO TO 300
120C-----------------------------------------------------------------------
121C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
122C-----------------------------------------------------------------------
123 IF (zr.LT.0.0d0) GO TO 300
124 kdflg = 1
125 yr(i)=zeror
126 yi(i)=zeroi
127 nz=nz+1
128 IF (i.EQ.1) GO TO 70
129 IF ((yr(i-1).EQ.zeror).AND.(yi(i-1).EQ.zeroi)) GO TO 70
130 yr(i-1)=zeror
131 yi(i-1)=zeroi
132 nz=nz+1
133 70 CONTINUE
134 i = n
135 75 CONTINUE
136 razr = 1.0d0/xzabs(zrr,zri)
137 str = zrr*razr
138 sti = -zri*razr
139 rzr = (str+str)*razr
140 rzi = (sti+sti)*razr
141 ckr = fn*rzr
142 cki = fn*rzi
143 ib = i + 1
144 IF (n.LT.ib) GO TO 160
145C-----------------------------------------------------------------------
146C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
147C ON UNDERFLOW.
148C-----------------------------------------------------------------------
149 fn = fnu + dble(float(n-1))
150 ipard = 1
151 IF (mr.NE.0) ipard = 0
152 initd = 0
153 CALL zunik(zrr, zri, fn, 2, ipard, tol, initd, phidr, phidi,
154 * zet1dr, zet1di, zet2dr, zet2di, sumdr, sumdi, cwrkr(1,3),
155 * cwrki(1,3))
156 IF (kode.EQ.1) GO TO 80
157 str = zrr + zet2dr
158 sti = zri + zet2di
159 rast = fn/xzabs(str,sti)
160 str = str*rast*rast
161 sti = -sti*rast*rast
162 s1r = zet1dr - str
163 s1i = zet1di - sti
164 GO TO 90
165 80 CONTINUE
166 s1r = zet1dr - zet2dr
167 s1i = zet1di - zet2di
168 90 CONTINUE
169 rs1 = s1r
170 IF (dabs(rs1).GT.elim) GO TO 95
171 IF (dabs(rs1).LT.alim) GO TO 100
172C----------------------------------------------------------------------------
173C REFINE ESTIMATE AND TEST
174C-------------------------------------------------------------------------
175 aphi = xzabs(phidr,phidi)
176 rs1 = rs1+dlog(aphi)
177 IF (dabs(rs1).LT.elim) GO TO 100
178 95 CONTINUE
179 IF (dabs(rs1).GT.0.0d0) GO TO 300
180C-----------------------------------------------------------------------
181C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
182C-----------------------------------------------------------------------
183 IF (zr.LT.0.0d0) GO TO 300
184 nz = n
185 DO 96 i=1,n
186 yr(i) = zeror
187 yi(i) = zeroi
188 96 CONTINUE
189 RETURN
190C---------------------------------------------------------------------------
191C FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
192C----------------------------------------------------------------------------
193 100 CONTINUE
194 s1r = cyr(1)
195 s1i = cyi(1)
196 s2r = cyr(2)
197 s2i = cyi(2)
198 c1r = csrr(kflag)
199 ascle = bry(kflag)
200 DO 120 i=ib,n
201 c2r = s2r
202 c2i = s2i
203 s2r = ckr*c2r - cki*c2i + s1r
204 s2i = ckr*c2i + cki*c2r + s1i
205 s1r = c2r
206 s1i = c2i
207 ckr = ckr + rzr
208 cki = cki + rzi
209 c2r = s2r*c1r
210 c2i = s2i*c1r
211 yr(i) = c2r
212 yi(i) = c2i
213 IF (kflag.GE.3) GO TO 120
214 str = dabs(c2r)
215 sti = dabs(c2i)
216 c2m = dmax1(str,sti)
217 IF (c2m.LE.ascle) GO TO 120
218 kflag = kflag + 1
219 ascle = bry(kflag)
220 s1r = s1r*c1r
221 s1i = s1i*c1r
222 s2r = c2r
223 s2i = c2i
224 s1r = s1r*cssr(kflag)
225 s1i = s1i*cssr(kflag)
226 s2r = s2r*cssr(kflag)
227 s2i = s2i*cssr(kflag)
228 c1r = csrr(kflag)
229 120 CONTINUE
230 160 CONTINUE
231 IF (mr.EQ.0) RETURN
232C-----------------------------------------------------------------------
233C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
234C-----------------------------------------------------------------------
235 nz = 0
236 fmr = dble(float(mr))
237 sgn = -dsign(pi,fmr)
238C-----------------------------------------------------------------------
239C CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
240C-----------------------------------------------------------------------
241 csgni = sgn
242 inu = int(sngl(fnu))
243 fnf = fnu - dble(float(inu))
244 ifn = inu + n - 1
245 ang = fnf*sgn
246 cspnr = dcos(ang)
247 cspni = dsin(ang)
248 IF (mod(ifn,2).EQ.0) GO TO 170
249 cspnr = -cspnr
250 cspni = -cspni
251 170 CONTINUE
252 asc = bry(1)
253 iuf = 0
254 kk = n
255 kdflg = 1
256 ib = ib - 1
257 ic = ib - 1
258 DO 270 k=1,n
259 fn = fnu + dble(float(kk-1))
260C-----------------------------------------------------------------------
261C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
262C FUNCTION ABOVE
263C-----------------------------------------------------------------------
264 m=3
265 IF (n.GT.2) GO TO 175
266 172 CONTINUE
267 initd = init(j)
268 phidr = phir(j)
269 phidi = phii(j)
270 zet1dr = zeta1r(j)
271 zet1di = zeta1i(j)
272 zet2dr = zeta2r(j)
273 zet2di = zeta2i(j)
274 sumdr = sumr(j)
275 sumdi = sumi(j)
276 m = j
277 j = 3 - j
278 GO TO 180
279 175 CONTINUE
280 IF ((kk.EQ.n).AND.(ib.LT.n)) GO TO 180
281 IF ((kk.EQ.ib).OR.(kk.EQ.ic)) GO TO 172
282 initd = 0
283 180 CONTINUE
284 CALL zunik(zrr, zri, fn, 1, 0, tol, initd, phidr, phidi,
285 * zet1dr, zet1di, zet2dr, zet2di, sumdr, sumdi,
286 * cwrkr(1,m), cwrki(1,m))
287 IF (kode.EQ.1) GO TO 200
288 str = zrr + zet2dr
289 sti = zri + zet2di
290 rast = fn/xzabs(str,sti)
291 str = str*rast*rast
292 sti = -sti*rast*rast
293 s1r = -zet1dr + str
294 s1i = -zet1di + sti
295 GO TO 210
296 200 CONTINUE
297 s1r = -zet1dr + zet2dr
298 s1i = -zet1di + zet2di
299 210 CONTINUE
300C-----------------------------------------------------------------------
301C TEST FOR UNDERFLOW AND OVERFLOW
302C-----------------------------------------------------------------------
303 rs1 = s1r
304 IF (dabs(rs1).GT.elim) GO TO 260
305 IF (kdflg.EQ.1) iflag = 2
306 IF (dabs(rs1).LT.alim) GO TO 220
307C-----------------------------------------------------------------------
308C REFINE TEST AND SCALE
309C-----------------------------------------------------------------------
310 aphi = xzabs(phidr,phidi)
311 rs1 = rs1 + dlog(aphi)
312 IF (dabs(rs1).GT.elim) GO TO 260
313 IF (kdflg.EQ.1) iflag = 1
314 IF (rs1.LT.0.0d0) GO TO 220
315 IF (kdflg.EQ.1) iflag = 3
316 220 CONTINUE
317 str = phidr*sumdr - phidi*sumdi
318 sti = phidr*sumdi + phidi*sumdr
319 s2r = -csgni*sti
320 s2i = csgni*str
321 str = dexp(s1r)*cssr(iflag)
322 s1r = str*dcos(s1i)
323 s1i = str*dsin(s1i)
324 str = s2r*s1r - s2i*s1i
325 s2i = s2r*s1i + s2i*s1r
326 s2r = str
327 IF (iflag.NE.1) GO TO 230
328 CALL zuchk(s2r, s2i, nw, bry(1), tol)
329 IF (nw.EQ.0) GO TO 230
330 s2r = zeror
331 s2i = zeroi
332 230 CONTINUE
333 cyr(kdflg) = s2r
334 cyi(kdflg) = s2i
335 c2r = s2r
336 c2i = s2i
337 s2r = s2r*csrr(iflag)
338 s2i = s2i*csrr(iflag)
339C-----------------------------------------------------------------------
340C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
341C-----------------------------------------------------------------------
342 s1r = yr(kk)
343 s1i = yi(kk)
344 IF (kode.EQ.1) GO TO 250
345 CALL zs1s2(zrr, zri, s1r, s1i, s2r, s2i, nw, asc, alim, iuf)
346 nz = nz + nw
347 250 CONTINUE
348 yr(kk) = s1r*cspnr - s1i*cspni + s2r
349 yi(kk) = cspnr*s1i + cspni*s1r + s2i
350 kk = kk - 1
351 cspnr = -cspnr
352 cspni = -cspni
353 IF (c2r.NE.0.0d0 .OR. c2i.NE.0.0d0) GO TO 255
354 kdflg = 1
355 GO TO 270
356 255 CONTINUE
357 IF (kdflg.EQ.2) GO TO 275
358 kdflg = 2
359 GO TO 270
360 260 CONTINUE
361 IF (rs1.GT.0.0d0) GO TO 300
362 s2r = zeror
363 s2i = zeroi
364 GO TO 230
365 270 CONTINUE
366 k = n
367 275 CONTINUE
368 il = n - k
369 IF (il.EQ.0) RETURN
370C-----------------------------------------------------------------------
371C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
372C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
373C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
374C-----------------------------------------------------------------------
375 s1r = cyr(1)
376 s1i = cyi(1)
377 s2r = cyr(2)
378 s2i = cyi(2)
379 csr = csrr(iflag)
380 ascle = bry(iflag)
381 fn = dble(float(inu+il))
382 DO 290 i=1,il
383 c2r = s2r
384 c2i = s2i
385 s2r = s1r + (fn+fnf)*(rzr*c2r-rzi*c2i)
386 s2i = s1i + (fn+fnf)*(rzr*c2i+rzi*c2r)
387 s1r = c2r
388 s1i = c2i
389 fn = fn - 1.0d0
390 c2r = s2r*csr
391 c2i = s2i*csr
392 ckr = c2r
393 cki = c2i
394 c1r = yr(kk)
395 c1i = yi(kk)
396 IF (kode.EQ.1) GO TO 280
397 CALL zs1s2(zrr, zri, c1r, c1i, c2r, c2i, nw, asc, alim, iuf)
398 nz = nz + nw
399 280 CONTINUE
400 yr(kk) = c1r*cspnr - c1i*cspni + c2r
401 yi(kk) = c1r*cspni + c1i*cspnr + c2i
402 kk = kk - 1
403 cspnr = -cspnr
404 cspni = -cspni
405 IF (iflag.GE.3) GO TO 290
406 c2r = dabs(ckr)
407 c2i = dabs(cki)
408 c2m = dmax1(c2r,c2i)
409 IF (c2m.LE.ascle) GO TO 290
410 iflag = iflag + 1
411 ascle = bry(iflag)
412 s1r = s1r*csr
413 s1i = s1i*csr
414 s2r = ckr
415 s2i = cki
416 s1r = s1r*cssr(iflag)
417 s1i = s1i*cssr(iflag)
418 s2r = s2r*cssr(iflag)
419 s2i = s2i*cssr(iflag)
420 csr = csrr(iflag)
421 290 CONTINUE
422 RETURN
423 300 CONTINUE
424 nz = -1
425 RETURN
426 END
d1mach
double precision function d1mach(i)
Definition d1mach.f:23
mod
T mod(T x, T y)
Definition lo-mappers.h:294
xzabs
double precision function xzabs(zr, zi)
Definition xzabs.f:2
zs1s2
subroutine zs1s2(zrr, zri, s1r, s1i, s2r, s2i, nz, ascle, alim, iuf)
Definition zs1s2.f:3
zuchk
subroutine zuchk(yr, yi, nz, ascle, tol)
Definition zuchk.f:2
zunik
subroutine zunik(zrr, zri, fnu, ikflg, ipmtr, tol, init, phir, phii, zeta1r, zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
Definition zunik.f:3
zunk1
subroutine zunk1(zr, zi, fnu, kode, mr, n, yr, yi, nz, tol, elim, alim)
Definition zunk1.f:3