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
cunk1.f
Go to the documentation of this file.
1  SUBROUTINE cunk1(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CUNK1
3 C***REFER TO CBESK
4 C
5 C CUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
6 C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
7 C UNIFORM ASYMPTOTIC EXPANSION.
8 C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
9 C NZ=-1 MEANS AN OVERFLOW WILL OCCUR
10 C
11 C***ROUTINES CALLED CS1S2,CUCHK,CUNIK,R1MACH
12 C***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 /
24 C
25  kdflg = 1
26  nz = 0
27 C-----------------------------------------------------------------------
28 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
29 C THE UNDERFLOW LIMIT
30 C-----------------------------------------------------------------------
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
47 C-----------------------------------------------------------------------
48 C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
49 C-----------------------------------------------------------------------
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
62 C-----------------------------------------------------------------------
63 C TEST FOR UNDERFLOW AND OVERFLOW
64 C-----------------------------------------------------------------------
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
69 C-----------------------------------------------------------------------
70 C REFINE TEST AND SCALE
71 C-----------------------------------------------------------------------
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
79 C-----------------------------------------------------------------------
80 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
81 C EXPONENT EXTREMES
82 C-----------------------------------------------------------------------
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
100 C-----------------------------------------------------------------------
101 C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
102 C-----------------------------------------------------------------------
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
118 C-----------------------------------------------------------------------
119 C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW, SET SEQUENCE TO ZERO
120 C ON UNDERFLOW
121 C-----------------------------------------------------------------------
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
138 C-----------------------------------------------------------------------
139 C REFINE ESTIMATE AND TEST
140 C-----------------------------------------------------------------------
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
146 C-----------------------------------------------------------------------
147 C FOR X.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
148 C-----------------------------------------------------------------------
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
156 C-----------------------------------------------------------------------
157 C RECUR FORWARD FOR REMAINDER OF THE SEQUENCE
158 C-----------------------------------------------------------------------
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
187 C-----------------------------------------------------------------------
188 C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0
189 C-----------------------------------------------------------------------
190  nz = 0
191  fmr = float(mr)
192  sgn = -sign(pi,fmr)
193 C-----------------------------------------------------------------------
194 C CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
195 C-----------------------------------------------------------------------
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)
213 C-----------------------------------------------------------------------
214 C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
215 C FUNCTION ABOVE
216 C-----------------------------------------------------------------------
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
242 C-----------------------------------------------------------------------
243 C TEST FOR UNDERFLOW AND OVERFLOW
244 C-----------------------------------------------------------------------
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
249 C-----------------------------------------------------------------------
250 C REFINE TEST AND SCALE
251 C-----------------------------------------------------------------------
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)
272 C-----------------------------------------------------------------------
273 C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
274 C-----------------------------------------------------------------------
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
299 C-----------------------------------------------------------------------
300 C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
301 C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
302 C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
303 C-----------------------------------------------------------------------
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
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2