GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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)
3 C***BEGIN PROLOGUE ZUNK2
4 C***REFER TO ZBESK
5 C
6 C ZUNK2 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
7 C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
8 C UNIFORM ASYMPTOTIC EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN)
9 C WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR
10 C -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT
11 C HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC-
12 C ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
13 C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
14 C
15 C***ROUTINES CALLED ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,XZABS
16 C***END PROLOGUE ZUNK2
17 C COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC,
18 C *CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ,
19 C *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 /
45 C
46  kdflg = 1
47  nz = 0
48 C-----------------------------------------------------------------------
49 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
50 C THE UNDERFLOW LIMIT
51 C-----------------------------------------------------------------------
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
90 C-----------------------------------------------------------------------
91 C K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST
92 C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
93 C CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS
94 C-----------------------------------------------------------------------
95  j = 2
96  DO 80 i=1,n
97 C-----------------------------------------------------------------------
98 C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
99 C-----------------------------------------------------------------------
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
118 C-----------------------------------------------------------------------
119 C TEST FOR UNDERFLOW AND OVERFLOW
120 C-----------------------------------------------------------------------
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
125 C-----------------------------------------------------------------------
126 C REFINE TEST AND SCALE
127 C-----------------------------------------------------------------------
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
136 C-----------------------------------------------------------------------
137 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
138 C EXPONENT EXTREMES
139 C-----------------------------------------------------------------------
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
177 C-----------------------------------------------------------------------
178 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
179 C-----------------------------------------------------------------------
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
205 C-----------------------------------------------------------------------
206 C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
207 C ON UNDERFLOW.
208 C-----------------------------------------------------------------------
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
230 C----------------------------------------------------------------------------
231 C REFINE ESTIMATE AND TEST
232 C-------------------------------------------------------------------------
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
238 C-----------------------------------------------------------------------
239 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
240 C-----------------------------------------------------------------------
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
287 C-----------------------------------------------------------------------
288 C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
289 C-----------------------------------------------------------------------
290  nz = 0
291  fmr = dble(float(mr))
292  sgn = -dsign(pi,fmr)
293 C-----------------------------------------------------------------------
294 C CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
295 C-----------------------------------------------------------------------
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
306 C-----------------------------------------------------------------------
307 C CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
308 C COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
309 C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
310 C CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
311 C-----------------------------------------------------------------------
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))
328 C-----------------------------------------------------------------------
329 C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
330 C FUNCTION ABOVE
331 C-----------------------------------------------------------------------
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
368 C-----------------------------------------------------------------------
369 C TEST FOR UNDERFLOW AND OVERFLOW
370 C-----------------------------------------------------------------------
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
375 C-----------------------------------------------------------------------
376 C REFINE TEST AND SCALE
377 C-----------------------------------------------------------------------
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)
415 C-----------------------------------------------------------------------
416 C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
417 C-----------------------------------------------------------------------
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
449 C-----------------------------------------------------------------------
450 C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
451 C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
452 C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
453 C-----------------------------------------------------------------------
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
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2