GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
zunk2.f
Go to the documentation of this file.
1 SUBROUTINE zunk2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
2 * ALIM)
3C***BEGIN PROLOGUE ZUNK2
4C***REFER TO ZBESK
5C
6C ZUNK2 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 EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
9C WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
10C -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
11C HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
12C ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
13C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
14C
15C***ROUTINES CALLED ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,XZABS
16C***END PROLOGUE ZUNK2
17C COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC,
18C *CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ,
19C *S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR
20 DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGDI,
21 * argdr, argi, argr, asc, ascle, asumdi, asumdr, asumi, asumr,
22 * bry, bsumdi, bsumdr, bsumi, bsumr, car, cipi, cipr, cki, ckr,
23 * coner, crsc, cr1i, cr1r, cr2i, cr2r, cscl, csgni, csi,
24 * cspni, cspnr, csr, csrr, cssr, cyi, cyr, c1i, c1r, c2i, c2m,
25 * c2r, daii, dair, elim, fmr, fn, fnf, fnu, hpi, phidi, phidr,
26 * phii, phir, pi, pti, ptr, rast, razr, rs1, rzi, rzr, sar, sgn,
27 * sti, str, s1i, s1r, s2i, s2r, tol, yi, yr, yy, zbi, zbr, zeroi,
28 * zeror, zeta1i, zeta1r, zeta2i, zeta2r, zet1di, zet1dr, zet2di,
29 * zet2dr, zi, zni, znr, zr, zri, zrr, d1mach, xzabs
30 INTEGER I, IB, IFLAG, IFN, IL, IN, INU, IUF, K, KDFLG, KFLAG, KK,
31 * kode, mr, n, nai, ndai, nw, nz, idum, j, ipard, ic
32 dimension bry(3), yr(n), yi(n), asumr(2), asumi(2), bsumr(2),
33 * bsumi(2), phir(2), phii(2), argr(2), argi(2), zeta1r(2),
34 * zeta1i(2), zeta2r(2), zeta2i(2), cyr(2), cyi(2), cipr(4),
35 * cipi(4), cssr(3), csrr(3)
36 DATA zeror,zeroi,coner,cr1r,cr1i,cr2r,cr2i /
37 1 0.0d0, 0.0d0, 1.0d0,
38 1 1.0d0,1.73205080756887729d0 , -0.5d0,-8.66025403784438647d-01 /
39 DATA hpi, pi, aic /
40 1 1.57079632679489662d+00, 3.14159265358979324d+00,
41 1 1.26551212348464539d+00/
42 DATA cipr(1),cipi(1),cipr(2),cipi(2),cipr(3),cipi(3),cipr(4),
43 * cipi(4) /
44 1 1.0d0,0.0d0 , 0.0d0,-1.0d0 , -1.0d0,0.0d0 , 0.0d0,1.0d0 /
45C
46 kdflg = 1
47 nz = 0
48C-----------------------------------------------------------------------
49C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
50C THE UNDERFLOW LIMIT
51C-----------------------------------------------------------------------
52 cscl = 1.0d0/tol
53 crsc = tol
54 cssr(1) = cscl
55 cssr(2) = coner
56 cssr(3) = crsc
57 csrr(1) = crsc
58 csrr(2) = coner
59 csrr(3) = cscl
60 bry(1) = 1.0d+3*d1mach(1)/tol
61 bry(2) = 1.0d0/bry(1)
62 bry(3) = d1mach(2)
63 zrr = zr
64 zri = zi
65 IF (zr.GE.0.0d0) GO TO 10
66 zrr = -zr
67 zri = -zi
68 10 CONTINUE
69 yy = zri
70 znr = zri
71 zni = -zrr
72 zbr = zrr
73 zbi = zri
74 inu = int(sngl(fnu))
75 fnf = fnu - dble(float(inu))
76 ang = -hpi*fnf
77 car = dcos(ang)
78 sar = dsin(ang)
79 c2r = hpi*sar
80 c2i = -hpi*car
81 kk = mod(inu,4) + 1
82 str = c2r*cipr(kk) - c2i*cipi(kk)
83 sti = c2r*cipi(kk) + c2i*cipr(kk)
84 csr = cr1r*str - cr1i*sti
85 csi = cr1r*sti + cr1i*str
86 IF (yy.GT.0.0d0) GO TO 20
87 znr = -znr
88 zbi = -zbi
89 20 CONTINUE
90C-----------------------------------------------------------------------
91C K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
92C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
93C CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
94C-----------------------------------------------------------------------
95 j = 2
96 DO 80 i=1,n
97C-----------------------------------------------------------------------
98C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
99C-----------------------------------------------------------------------
100 j = 3 - j
101 fn = fnu + dble(float(i-1))
102 CALL zunhj(znr, zni, fn, 0, tol, phir(j), phii(j), argr(j),
103 * argi(j), zeta1r(j), zeta1i(j), zeta2r(j), zeta2i(j), asumr(j),
104 * asumi(j), bsumr(j), bsumi(j))
105 IF (kode.EQ.1) GO TO 30
106 str = zbr + zeta2r(j)
107 sti = zbi + zeta2i(j)
108 rast = fn/xzabs(str,sti)
109 str = str*rast*rast
110 sti = -sti*rast*rast
111 s1r = zeta1r(j) - str
112 s1i = zeta1i(j) - sti
113 GO TO 40
114 30 CONTINUE
115 s1r = zeta1r(j) - zeta2r(j)
116 s1i = zeta1i(j) - zeta2i(j)
117 40 CONTINUE
118C-----------------------------------------------------------------------
119C TEST FOR UNDERFLOW AND OVERFLOW
120C-----------------------------------------------------------------------
121 rs1 = s1r
122 IF (dabs(rs1).GT.elim) GO TO 70
123 IF (kdflg.EQ.1) kflag = 2
124 IF (dabs(rs1).LT.alim) GO TO 50
125C-----------------------------------------------------------------------
126C REFINE TEST AND SCALE
127C-----------------------------------------------------------------------
128 aphi = xzabs(phir(j),phii(j))
129 aarg = xzabs(argr(j),argi(j))
130 rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
131 IF (dabs(rs1).GT.elim) GO TO 70
132 IF (kdflg.EQ.1) kflag = 1
133 IF (rs1.LT.0.0d0) GO TO 50
134 IF (kdflg.EQ.1) kflag = 3
135 50 CONTINUE
136C-----------------------------------------------------------------------
137C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
138C EXPONENT EXTREMES
139C-----------------------------------------------------------------------
140 c2r = argr(j)*cr2r - argi(j)*cr2i
141 c2i = argr(j)*cr2i + argi(j)*cr2r
142 CALL zairy(c2r, c2i, 0, 2, air, aii, nai, idum)
143 CALL zairy(c2r, c2i, 1, 2, dair, daii, ndai, idum)
144 str = dair*bsumr(j) - daii*bsumi(j)
145 sti = dair*bsumi(j) + daii*bsumr(j)
146 ptr = str*cr2r - sti*cr2i
147 pti = str*cr2i + sti*cr2r
148 str = ptr + (air*asumr(j)-aii*asumi(j))
149 sti = pti + (air*asumi(j)+aii*asumr(j))
150 ptr = str*phir(j) - sti*phii(j)
151 pti = str*phii(j) + sti*phir(j)
152 s2r = ptr*csr - pti*csi
153 s2i = ptr*csi + pti*csr
154 str = dexp(s1r)*cssr(kflag)
155 s1r = str*dcos(s1i)
156 s1i = str*dsin(s1i)
157 str = s2r*s1r - s2i*s1i
158 s2i = s1r*s2i + s2r*s1i
159 s2r = str
160 IF (kflag.NE.1) GO TO 60
161 CALL zuchk(s2r, s2i, nw, bry(1), tol)
162 IF (nw.NE.0) GO TO 70
163 60 CONTINUE
164 IF (yy.LE.0.0d0) s2i = -s2i
165 cyr(kdflg) = s2r
166 cyi(kdflg) = s2i
167 yr(i) = s2r*csrr(kflag)
168 yi(i) = s2i*csrr(kflag)
169 str = csi
170 csi = -csr
171 csr = str
172 IF (kdflg.EQ.2) GO TO 85
173 kdflg = 2
174 GO TO 80
175 70 CONTINUE
176 IF (rs1.GT.0.0d0) GO TO 320
177C-----------------------------------------------------------------------
178C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
179C-----------------------------------------------------------------------
180 IF (zr.LT.0.0d0) GO TO 320
181 kdflg = 1
182 yr(i)=zeror
183 yi(i)=zeroi
184 nz=nz+1
185 str = csi
186 csi =-csr
187 csr = str
188 IF (i.EQ.1) GO TO 80
189 IF ((yr(i-1).EQ.zeror).AND.(yi(i-1).EQ.zeroi)) GO TO 80
190 yr(i-1)=zeror
191 yi(i-1)=zeroi
192 nz=nz+1
193 80 CONTINUE
194 i = n
195 85 CONTINUE
196 razr = 1.0d0/xzabs(zrr,zri)
197 str = zrr*razr
198 sti = -zri*razr
199 rzr = (str+str)*razr
200 rzi = (sti+sti)*razr
201 ckr = fn*rzr
202 cki = fn*rzi
203 ib = i + 1
204 IF (n.LT.ib) GO TO 180
205C-----------------------------------------------------------------------
206C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
207C ON UNDERFLOW.
208C-----------------------------------------------------------------------
209 fn = fnu + dble(float(n-1))
210 ipard = 1
211 IF (mr.NE.0) ipard = 0
212 CALL zunhj(znr, zni, fn, ipard, tol, phidr, phidi, argdr, argdi,
213 * zet1dr, zet1di, zet2dr, zet2di, asumdr, asumdi, bsumdr, bsumdi)
214 IF (kode.EQ.1) GO TO 90
215 str = zbr + zet2dr
216 sti = zbi + zet2di
217 rast = fn/xzabs(str,sti)
218 str = str*rast*rast
219 sti = -sti*rast*rast
220 s1r = zet1dr - str
221 s1i = zet1di - sti
222 GO TO 100
223 90 CONTINUE
224 s1r = zet1dr - zet2dr
225 s1i = zet1di - zet2di
226 100 CONTINUE
227 rs1 = s1r
228 IF (dabs(rs1).GT.elim) GO TO 105
229 IF (dabs(rs1).LT.alim) GO TO 120
230C----------------------------------------------------------------------------
231C REFINE ESTIMATE AND TEST
232C-------------------------------------------------------------------------
233 aphi = xzabs(phidr,phidi)
234 rs1 = rs1+dlog(aphi)
235 IF (dabs(rs1).LT.elim) GO TO 120
236 105 CONTINUE
237 IF (rs1.GT.0.0d0) GO TO 320
238C-----------------------------------------------------------------------
239C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
240C-----------------------------------------------------------------------
241 IF (zr.LT.0.0d0) GO TO 320
242 nz = n
243 DO 106 i=1,n
244 yr(i) = zeror
245 yi(i) = zeroi
246 106 CONTINUE
247 RETURN
248 120 CONTINUE
249 s1r = cyr(1)
250 s1i = cyi(1)
251 s2r = cyr(2)
252 s2i = cyi(2)
253 c1r = csrr(kflag)
254 ascle = bry(kflag)
255 DO 130 i=ib,n
256 c2r = s2r
257 c2i = s2i
258 s2r = ckr*c2r - cki*c2i + s1r
259 s2i = ckr*c2i + cki*c2r + s1i
260 s1r = c2r
261 s1i = c2i
262 ckr = ckr + rzr
263 cki = cki + rzi
264 c2r = s2r*c1r
265 c2i = s2i*c1r
266 yr(i) = c2r
267 yi(i) = c2i
268 IF (kflag.GE.3) GO TO 130
269 str = dabs(c2r)
270 sti = dabs(c2i)
271 c2m = dmax1(str,sti)
272 IF (c2m.LE.ascle) GO TO 130
273 kflag = kflag + 1
274 ascle = bry(kflag)
275 s1r = s1r*c1r
276 s1i = s1i*c1r
277 s2r = c2r
278 s2i = c2i
279 s1r = s1r*cssr(kflag)
280 s1i = s1i*cssr(kflag)
281 s2r = s2r*cssr(kflag)
282 s2i = s2i*cssr(kflag)
283 c1r = csrr(kflag)
284 130 CONTINUE
285 180 CONTINUE
286 IF (mr.EQ.0) RETURN
287C-----------------------------------------------------------------------
288C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
289C-----------------------------------------------------------------------
290 nz = 0
291 fmr = dble(float(mr))
292 sgn = -dsign(pi,fmr)
293C-----------------------------------------------------------------------
294C CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
295C-----------------------------------------------------------------------
296 csgni = sgn
297 IF (yy.LE.0.0d0) csgni = -csgni
298 ifn = inu + n - 1
299 ang = fnf*sgn
300 cspnr = dcos(ang)
301 cspni = dsin(ang)
302 IF (mod(ifn,2).EQ.0) GO TO 190
303 cspnr = -cspnr
304 cspni = -cspni
305 190 CONTINUE
306C-----------------------------------------------------------------------
307C CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
308C COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
309C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
310C CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
311C-----------------------------------------------------------------------
312 csr = sar*csgni
313 csi = car*csgni
314 in = mod(ifn,4) + 1
315 c2r = cipr(in)
316 c2i = cipi(in)
317 str = csr*c2r + csi*c2i
318 csi = -csr*c2i + csi*c2r
319 csr = str
320 asc = bry(1)
321 iuf = 0
322 kk = n
323 kdflg = 1
324 ib = ib - 1
325 ic = ib - 1
326 DO 290 k=1,n
327 fn = fnu + dble(float(kk-1))
328C-----------------------------------------------------------------------
329C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
330C FUNCTION ABOVE
331C-----------------------------------------------------------------------
332 IF (n.GT.2) GO TO 175
333 172 CONTINUE
334 phidr = phir(j)
335 phidi = phii(j)
336 argdr = argr(j)
337 argdi = argi(j)
338 zet1dr = zeta1r(j)
339 zet1di = zeta1i(j)
340 zet2dr = zeta2r(j)
341 zet2di = zeta2i(j)
342 asumdr = asumr(j)
343 asumdi = asumi(j)
344 bsumdr = bsumr(j)
345 bsumdi = bsumi(j)
346 j = 3 - j
347 GO TO 210
348 175 CONTINUE
349 IF ((kk.EQ.n).AND.(ib.LT.n)) GO TO 210
350 IF ((kk.EQ.ib).OR.(kk.EQ.ic)) GO TO 172
351 CALL zunhj(znr, zni, fn, 0, tol, phidr, phidi, argdr,
352 * argdi, zet1dr, zet1di, zet2dr, zet2di, asumdr,
353 * asumdi, bsumdr, bsumdi)
354 210 CONTINUE
355 IF (kode.EQ.1) GO TO 220
356 str = zbr + zet2dr
357 sti = zbi + zet2di
358 rast = fn/xzabs(str,sti)
359 str = str*rast*rast
360 sti = -sti*rast*rast
361 s1r = -zet1dr + str
362 s1i = -zet1di + sti
363 GO TO 230
364 220 CONTINUE
365 s1r = -zet1dr + zet2dr
366 s1i = -zet1di + zet2di
367 230 CONTINUE
368C-----------------------------------------------------------------------
369C TEST FOR UNDERFLOW AND OVERFLOW
370C-----------------------------------------------------------------------
371 rs1 = s1r
372 IF (dabs(rs1).GT.elim) GO TO 280
373 IF (kdflg.EQ.1) iflag = 2
374 IF (dabs(rs1).LT.alim) GO TO 240
375C-----------------------------------------------------------------------
376C REFINE TEST AND SCALE
377C-----------------------------------------------------------------------
378 aphi = xzabs(phidr,phidi)
379 aarg = xzabs(argdr,argdi)
380 rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
381 IF (dabs(rs1).GT.elim) GO TO 280
382 IF (kdflg.EQ.1) iflag = 1
383 IF (rs1.LT.0.0d0) GO TO 240
384 IF (kdflg.EQ.1) iflag = 3
385 240 CONTINUE
386 CALL zairy(argdr, argdi, 0, 2, air, aii, nai, idum)
387 CALL zairy(argdr, argdi, 1, 2, dair, daii, ndai, idum)
388 str = dair*bsumdr - daii*bsumdi
389 sti = dair*bsumdi + daii*bsumdr
390 str = str + (air*asumdr-aii*asumdi)
391 sti = sti + (air*asumdi+aii*asumdr)
392 ptr = str*phidr - sti*phidi
393 pti = str*phidi + sti*phidr
394 s2r = ptr*csr - pti*csi
395 s2i = ptr*csi + pti*csr
396 str = dexp(s1r)*cssr(iflag)
397 s1r = str*dcos(s1i)
398 s1i = str*dsin(s1i)
399 str = s2r*s1r - s2i*s1i
400 s2i = s2r*s1i + s2i*s1r
401 s2r = str
402 IF (iflag.NE.1) GO TO 250
403 CALL zuchk(s2r, s2i, nw, bry(1), tol)
404 IF (nw.EQ.0) GO TO 250
405 s2r = zeror
406 s2i = zeroi
407 250 CONTINUE
408 IF (yy.LE.0.0d0) s2i = -s2i
409 cyr(kdflg) = s2r
410 cyi(kdflg) = s2i
411 c2r = s2r
412 c2i = s2i
413 s2r = s2r*csrr(iflag)
414 s2i = s2i*csrr(iflag)
415C-----------------------------------------------------------------------
416C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
417C-----------------------------------------------------------------------
418 s1r = yr(kk)
419 s1i = yi(kk)
420 IF (kode.EQ.1) GO TO 270
421 CALL zs1s2(zrr, zri, s1r, s1i, s2r, s2i, nw, asc, alim, iuf)
422 nz = nz + nw
423 270 CONTINUE
424 yr(kk) = s1r*cspnr - s1i*cspni + s2r
425 yi(kk) = s1r*cspni + s1i*cspnr + s2i
426 kk = kk - 1
427 cspnr = -cspnr
428 cspni = -cspni
429 str = csi
430 csi = -csr
431 csr = str
432 IF (c2r.NE.0.0d0 .OR. c2i.NE.0.0d0) GO TO 255
433 kdflg = 1
434 GO TO 290
435 255 CONTINUE
436 IF (kdflg.EQ.2) GO TO 295
437 kdflg = 2
438 GO TO 290
439 280 CONTINUE
440 IF (rs1.GT.0.0d0) GO TO 320
441 s2r = zeror
442 s2i = zeroi
443 GO TO 250
444 290 CONTINUE
445 k = n
446 295 CONTINUE
447 il = n - k
448 IF (il.EQ.0) RETURN
449C-----------------------------------------------------------------------
450C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
451C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
452C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
453C-----------------------------------------------------------------------
454 s1r = cyr(1)
455 s1i = cyi(1)
456 s2r = cyr(2)
457 s2i = cyi(2)
458 csr = csrr(iflag)
459 ascle = bry(iflag)
460 fn = dble(float(inu+il))
461 DO 310 i=1,il
462 c2r = s2r
463 c2i = s2i
464 s2r = s1r + (fn+fnf)*(rzr*c2r-rzi*c2i)
465 s2i = s1i + (fn+fnf)*(rzr*c2i+rzi*c2r)
466 s1r = c2r
467 s1i = c2i
468 fn = fn - 1.0d0
469 c2r = s2r*csr
470 c2i = s2i*csr
471 ckr = c2r
472 cki = c2i
473 c1r = yr(kk)
474 c1i = yi(kk)
475 IF (kode.EQ.1) GO TO 300
476 CALL zs1s2(zrr, zri, c1r, c1i, c2r, c2i, nw, asc, alim, iuf)
477 nz = nz + nw
478 300 CONTINUE
479 yr(kk) = c1r*cspnr - c1i*cspni + c2r
480 yi(kk) = c1r*cspni + c1i*cspnr + c2i
481 kk = kk - 1
482 cspnr = -cspnr
483 cspni = -cspni
484 IF (iflag.GE.3) GO TO 310
485 c2r = dabs(ckr)
486 c2i = dabs(cki)
487 c2m = dmax1(c2r,c2i)
488 IF (c2m.LE.ascle) GO TO 310
489 iflag = iflag + 1
490 ascle = bry(iflag)
491 s1r = s1r*csr
492 s1i = s1i*csr
493 s2r = ckr
494 s2i = cki
495 s1r = s1r*cssr(iflag)
496 s1i = s1i*cssr(iflag)
497 s2r = s2r*cssr(iflag)
498 s2i = s2i*cssr(iflag)
499 csr = csrr(iflag)
500 310 CONTINUE
501 RETURN
502 320 CONTINUE
503 nz = -1
504 RETURN
505 END
double precision function d1mach(i)
Definition d1mach.f:23
T mod(T x, T y)
Definition lo-mappers.h:294
double precision function xzabs(zr, zi)
Definition xzabs.f:2
subroutine zairy(zr, zi, id, kode, air, aii, nz, ierr)
Definition zairy.f:2
subroutine zs1s2(zrr, zri, s1r, s1i, s2r, s2i, nz, ascle, alim, iuf)
Definition zs1s2.f:3
subroutine zuchk(yr, yi, nz, ascle, tol)
Definition zuchk.f:2
subroutine zunhj(zr, zi, fnu, ipmtr, tol, phir, phii, argr, argi, zeta1r, zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
Definition zunhj.f:3
subroutine zunk2(zr, zi, fnu, kode, mr, n, yr, yi, nz, tol, elim, alim)
Definition zunk2.f:3