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
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)
3 C***BEGIN PROLOGUE ZUNK1
4 C***REFER TO ZBESK
5 C
6 C ZUNK1 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 EXPANSION.
9 C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
10 C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
11 C
12 C***ROUTINES CALLED ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,XZABS
13 C***END PROLOGUE ZUNK1
14 C COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
15 C *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 /
30 C
31  kdflg = 1
32  nz = 0
33 C-----------------------------------------------------------------------
34 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
35 C THE UNDERFLOW LIMIT
36 C-----------------------------------------------------------------------
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
56 C-----------------------------------------------------------------------
57 C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
58 C-----------------------------------------------------------------------
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
79 C-----------------------------------------------------------------------
80 C TEST FOR UNDERFLOW AND OVERFLOW
81 C-----------------------------------------------------------------------
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
85 C-----------------------------------------------------------------------
86 C REFINE TEST AND SCALE
87 C-----------------------------------------------------------------------
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
95 C-----------------------------------------------------------------------
96 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
97 C EXPONENT EXTREMES
98 C-----------------------------------------------------------------------
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
120 C-----------------------------------------------------------------------
121 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
122 C-----------------------------------------------------------------------
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
145 C-----------------------------------------------------------------------
146 C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
147 C ON UNDERFLOW.
148 C-----------------------------------------------------------------------
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
172 C----------------------------------------------------------------------------
173 C REFINE ESTIMATE AND TEST
174 C-------------------------------------------------------------------------
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
180 C-----------------------------------------------------------------------
181 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
182 C-----------------------------------------------------------------------
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
190 C---------------------------------------------------------------------------
191 C FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
192 C----------------------------------------------------------------------------
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
232 C-----------------------------------------------------------------------
233 C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
234 C-----------------------------------------------------------------------
235  nz = 0
236  fmr = dble(float(mr))
237  sgn = -dsign(pi,fmr)
238 C-----------------------------------------------------------------------
239 C CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
240 C-----------------------------------------------------------------------
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))
260 C-----------------------------------------------------------------------
261 C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
262 C FUNCTION ABOVE
263 C-----------------------------------------------------------------------
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
300 C-----------------------------------------------------------------------
301 C TEST FOR UNDERFLOW AND OVERFLOW
302 C-----------------------------------------------------------------------
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
307 C-----------------------------------------------------------------------
308 C REFINE TEST AND SCALE
309 C-----------------------------------------------------------------------
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)
339 C-----------------------------------------------------------------------
340 C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
341 C-----------------------------------------------------------------------
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
370 C-----------------------------------------------------------------------
371 C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
372 C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
373 C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
374 C-----------------------------------------------------------------------
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
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2