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
cbknu.f
Go to the documentation of this file.
1  SUBROUTINE cbknu(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CBKNU
3 C***REFER TO CBESI,CBESK,CAIRY,CBESH
4 C
5 C CBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE
6 C
7 C***ROUTINES CALLED CKSCL,CSHCH,GAMLN,I1MACH,R1MACH,CUCHK
8 C***END PROLOGUE CBKNU
9 C
10  COMPLEX cch, ck, coef, cone, crsc, cs, cscl, csh, csr, css, ctwo,
11  * cz, czero, f, fmu, p, pt, p1, p2, q, rz, smu, st, s1, s2, y, z,
12  * zd, celm, cy
13  REAL aa, ak, alim, ascle, a1, a2, bb, bk, bry, caz, cc, dnu,
14  * dnu2, elim, etest, fc, fhs, fk, fks, fnu, fpi, g1, g2, hpi, pi,
15  * p2i, p2m, p2r, rk, rthpi, r1, s, spi, tm, tol, tth, t1, t2, xx,
16  * yy, gamln, r1mach, helim, elm, xd, yd, alas, as
17  INTEGER i, idum, iflag, inu, k, kflag, kk, kmax, kode, koded, n,
18  * nz, i1mach, nw, j, ic, inub
19  dimension bry(3), cc(8), css(3), csr(3), y(n), cy(2)
20 C
21  DATA kmax / 30 /
22  DATA r1 / 2.0e0 /
23  DATA czero,cone,ctwo /(0.0e0,0.0e0),(1.0e0,0.0e0),(2.0e0,0.0e0)/
24 C
25  DATA pi, rthpi, spi ,hpi, fpi, tth /
26  1 3.14159265358979324e0, 1.25331413731550025e0,
27  2 1.90985931710274403e0, 1.57079632679489662e0,
28  3 1.89769999331517738e0, 6.66666666666666666e-01/
29 C
30  DATA cc(1), cc(2), cc(3), cc(4), cc(5), cc(6), cc(7), cc(8)/
31  1 5.77215664901532861e-01, -4.20026350340952355e-02,
32  2 -4.21977345555443367e-02, 7.21894324666309954e-03,
33  3 -2.15241674114950973e-04, -2.01348547807882387e-05,
34  4 1.13302723198169588e-06, 6.11609510448141582e-09/
35 C
36  xx = REAL(z)
37  yy = aimag(z)
38  caz = cabs(z)
39  cscl = cmplx(1.0e0/tol,0.0e0)
40  crsc = cmplx(tol,0.0e0)
41  css(1) = cscl
42  css(2) = cone
43  css(3) = crsc
44  csr(1) = crsc
45  csr(2) = cone
46  csr(3) = cscl
47  bry(1) = 1.0e+3*r1mach(1)/tol
48  bry(2) = 1.0e0/bry(1)
49  bry(3) = r1mach(2)
50  nz = 0
51  iflag = 0
52  koded = kode
53  rz = ctwo/z
54  inu = int(fnu+0.5e0)
55  dnu = fnu - float(inu)
56  IF (abs(dnu).EQ.0.5e0) go to 110
57  dnu2 = 0.0e0
58  IF (abs(dnu).GT.tol) dnu2 = dnu*dnu
59  IF (caz.GT.r1) go to 110
60 C-----------------------------------------------------------------------
61 C SERIES FOR CABS(Z).LE.R1
62 C-----------------------------------------------------------------------
63  fc = 1.0e0
64  smu = clog(rz)
65  fmu = smu*cmplx(dnu,0.0e0)
66  CALL cshch(fmu, csh, cch)
67  IF (dnu.EQ.0.0e0) go to 10
68  fc = dnu*pi
69  fc = fc/sin(fc)
70  smu = csh*cmplx(1.0e0/dnu,0.0e0)
71  10 CONTINUE
72  a2 = 1.0e0 + dnu
73 C-----------------------------------------------------------------------
74 C GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
75 C-----------------------------------------------------------------------
76  t2 = exp(-gamln(a2,idum))
77  t1 = 1.0e0/(t2*fc)
78  IF (abs(dnu).GT.0.1e0) go to 40
79 C-----------------------------------------------------------------------
80 C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
81 C-----------------------------------------------------------------------
82  ak = 1.0e0
83  s = cc(1)
84  DO 20 k=2,8
85  ak = ak*dnu2
86  tm = cc(k)*ak
87  s = s + tm
88  IF (abs(tm).LT.tol) go to 30
89  20 CONTINUE
90  30 g1 = -s
91  go to 50
92  40 CONTINUE
93  g1 = (t1-t2)/(dnu+dnu)
94  50 CONTINUE
95  g2 = 0.5e0*(t1+t2)*fc
96  g1 = g1*fc
97  f = cmplx(g1,0.0e0)*cch + smu*cmplx(g2,0.0e0)
98  pt = cexp(fmu)
99  p = cmplx(0.5e0/t2,0.0e0)*pt
100  q = cmplx(0.5e0/t1,0.0e0)/pt
101  s1 = f
102  s2 = p
103  ak = 1.0e0
104  a1 = 1.0e0
105  ck = cone
106  bk = 1.0e0 - dnu2
107  IF (inu.GT.0 .OR. n.GT.1) go to 80
108 C-----------------------------------------------------------------------
109 C GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
110 C-----------------------------------------------------------------------
111  IF (caz.LT.tol) go to 70
112  cz = z*z*cmplx(0.25e0,0.0e0)
113  t1 = 0.25e0*caz*caz
114  60 CONTINUE
115  f = (f*cmplx(ak,0.0e0)+p+q)*cmplx(1.0e0/bk,0.0e0)
116  p = p*cmplx(1.0e0/(ak-dnu),0.0e0)
117  q = q*cmplx(1.0e0/(ak+dnu),0.0e0)
118  rk = 1.0e0/ak
119  ck = ck*cz*cmplx(rk,0.0)
120  s1 = s1 + ck*f
121  a1 = a1*t1*rk
122  bk = bk + ak + ak + 1.0e0
123  ak = ak + 1.0e0
124  IF (a1.GT.tol) go to 60
125  70 CONTINUE
126  y(1) = s1
127  IF (koded.EQ.1) RETURN
128  y(1) = s1*cexp(z)
129  RETURN
130 C-----------------------------------------------------------------------
131 C GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
132 C-----------------------------------------------------------------------
133  80 CONTINUE
134  IF (caz.LT.tol) go to 100
135  cz = z*z*cmplx(0.25e0,0.0e0)
136  t1 = 0.25e0*caz*caz
137  90 CONTINUE
138  f = (f*cmplx(ak,0.0e0)+p+q)*cmplx(1.0e0/bk,0.0e0)
139  p = p*cmplx(1.0e0/(ak-dnu),0.0e0)
140  q = q*cmplx(1.0e0/(ak+dnu),0.0e0)
141  rk = 1.0e0/ak
142  ck = ck*cz*cmplx(rk,0.0e0)
143  s1 = s1 + ck*f
144  s2 = s2 + ck*(p-f*cmplx(ak,0.0e0))
145  a1 = a1*t1*rk
146  bk = bk + ak + ak + 1.0e0
147  ak = ak + 1.0e0
148  IF (a1.GT.tol) go to 90
149  100 CONTINUE
150  kflag = 2
151  bk = REAL(smu)
152  a1 = fnu + 1.0e0
153  ak = a1*abs(bk)
154  IF (ak.GT.alim) kflag = 3
155  p2 = s2*css(kflag)
156  s2 = p2*rz
157  s1 = s1*css(kflag)
158  IF (koded.EQ.1) go to 210
159  f = cexp(z)
160  s1 = s1*f
161  s2 = s2*f
162  go to 210
163 C-----------------------------------------------------------------------
164 C IFLAG=0 MEANS NO UNDERFLOW OCCURRED
165 C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
166 C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
167 C RECURSION
168 C-----------------------------------------------------------------------
169  110 CONTINUE
170  coef = cmplx(rthpi,0.0e0)/csqrt(z)
171  kflag = 2
172  IF (koded.EQ.2) go to 120
173  IF (xx.GT.alim) go to 290
174 C BLANK LINE
175  a1 = exp(-xx)*REAL(css(kflag))
176  pt = cmplx(a1,0.0e0)*cmplx(cos(yy),-sin(yy))
177  coef = coef*pt
178  120 CONTINUE
179  IF (abs(dnu).EQ.0.5e0) go to 300
180 C-----------------------------------------------------------------------
181 C MILLER ALGORITHM FOR CABS(Z).GT.R1
182 C-----------------------------------------------------------------------
183  ak = cos(pi*dnu)
184  ak = abs(ak)
185  IF (ak.EQ.0.0e0) go to 300
186  fhs = abs(0.25e0-dnu2)
187  IF (fhs.EQ.0.0e0) go to 300
188 C-----------------------------------------------------------------------
189 C COMPUTE R2=F(E). IF CABS(Z).GE.R2, USE FORWARD RECURRENCE TO
190 C DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
191 C 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(11))=
192 C TOL WHERE B IS THE BASE OF THE ARITHMETIC.
193 C-----------------------------------------------------------------------
194  t1 = float(i1mach(11)-1)*r1mach(5)*3.321928094e0
195  t1 = amax1(t1,12.0e0)
196  t1 = amin1(t1,60.0e0)
197  t2 = tth*t1 - 6.0e0
198  IF (xx.NE.0.0e0) go to 130
199  t1 = hpi
200  go to 140
201  130 CONTINUE
202  t1 = atan(yy/xx)
203  t1 = abs(t1)
204  140 CONTINUE
205  IF (t2.GT.caz) go to 170
206 C-----------------------------------------------------------------------
207 C FORWARD RECURRENCE LOOP WHEN CABS(Z).GE.R2
208 C-----------------------------------------------------------------------
209  etest = ak/(pi*caz*tol)
210  fk = 1.0e0
211  IF (etest.LT.1.0e0) go to 180
212  fks = 2.0e0
213  rk = caz + caz + 2.0e0
214  a1 = 0.0e0
215  a2 = 1.0e0
216  DO 150 i=1,kmax
217  ak = fhs/fks
218  bk = rk/(fk+1.0e0)
219  tm = a2
220  a2 = bk*a2 - ak*a1
221  a1 = tm
222  rk = rk + 2.0e0
223  fks = fks + fk + fk + 2.0e0
224  fhs = fhs + fk + fk
225  fk = fk + 1.0e0
226  tm = abs(a2)*fk
227  IF (etest.LT.tm) go to 160
228  150 CONTINUE
229  go to 310
230  160 CONTINUE
231  fk = fk + spi*t1*sqrt(t2/caz)
232  fhs = abs(0.25e0-dnu2)
233  go to 180
234  170 CONTINUE
235 C-----------------------------------------------------------------------
236 C COMPUTE BACKWARD INDEX K FOR CABS(Z).LT.R2
237 C-----------------------------------------------------------------------
238  a2 = sqrt(caz)
239  ak = fpi*ak/(tol*sqrt(a2))
240  aa = 3.0e0*t1/(1.0e0+caz)
241  bb = 14.7e0*t1/(28.0e0+caz)
242  ak = (alog(ak)+caz*cos(aa)/(1.0e0+0.008e0*caz))/cos(bb)
243  fk = 0.12125e0*ak*ak/caz + 1.5e0
244  180 CONTINUE
245  k = int(fk)
246 C-----------------------------------------------------------------------
247 C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
248 C-----------------------------------------------------------------------
249  fk = float(k)
250  fks = fk*fk
251  p1 = czero
252  p2 = cmplx(tol,0.0e0)
253  cs = p2
254  DO 190 i=1,k
255  a1 = fks - fk
256  a2 = (fks+fk)/(a1+fhs)
257  rk = 2.0e0/(fk+1.0e0)
258  t1 = (fk+xx)*rk
259  t2 = yy*rk
260  pt = p2
261  p2 = (p2*cmplx(t1,t2)-p1)*cmplx(a2,0.0e0)
262  p1 = pt
263  cs = cs + p2
264  fks = a1 - fk + 1.0e0
265  fk = fk - 1.0e0
266  190 CONTINUE
267 C-----------------------------------------------------------------------
268 C COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER
269 C SCALING
270 C-----------------------------------------------------------------------
271  tm = cabs(cs)
272  pt = cmplx(1.0e0/tm,0.0e0)
273  s1 = pt*p2
274  cs = conjg(cs)*pt
275  s1 = coef*s1*cs
276  IF (inu.GT.0 .OR. n.GT.1) go to 200
277  zd = z
278  IF(iflag.EQ.1) go to 270
279  go to 240
280  200 CONTINUE
281 C-----------------------------------------------------------------------
282 C COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING
283 C-----------------------------------------------------------------------
284  tm = cabs(p2)
285  pt = cmplx(1.0e0/tm,0.0e0)
286  p1 = pt*p1
287  p2 = conjg(p2)*pt
288  pt = p1*p2
289  s2 = s1*(cone+(cmplx(dnu+0.5e0,0.0e0)-pt)/z)
290 C-----------------------------------------------------------------------
291 C FORWARD RECURSION ON THE THREE TERM RECURSION RELATION WITH
292 C SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
293 C-----------------------------------------------------------------------
294  210 CONTINUE
295  ck = cmplx(dnu+1.0e0,0.0e0)*rz
296  IF (n.EQ.1) inu = inu - 1
297  IF (inu.GT.0) go to 220
298  IF (n.EQ.1) s1=s2
299  zd = z
300  IF(iflag.EQ.1) go to 270
301  go to 240
302  220 CONTINUE
303  inub = 1
304  IF (iflag.EQ.1) go to 261
305  225 CONTINUE
306  p1 = csr(kflag)
307  ascle = bry(kflag)
308  DO 230 i=inub,inu
309  st = s2
310  s2 = ck*s2 + s1
311  s1 = st
312  ck = ck + rz
313  IF (kflag.GE.3) go to 230
314  p2 = s2*p1
315  p2r = REAL(p2)
316  p2i = aimag(p2)
317  p2r = abs(p2r)
318  p2i = abs(p2i)
319  p2m = amax1(p2r,p2i)
320  IF (p2m.LE.ascle) go to 230
321  kflag = kflag + 1
322  ascle = bry(kflag)
323  s1 = s1*p1
324  s2 = p2
325  s1 = s1*css(kflag)
326  s2 = s2*css(kflag)
327  p1 = csr(kflag)
328  230 CONTINUE
329  IF (n.EQ.1) s1 = s2
330  240 CONTINUE
331  y(1) = s1*csr(kflag)
332  IF (n.EQ.1) RETURN
333  y(2) = s2*csr(kflag)
334  IF (n.EQ.2) RETURN
335  kk = 2
336  250 CONTINUE
337  kk = kk + 1
338  IF (kk.GT.n) RETURN
339  p1 = csr(kflag)
340  ascle = bry(kflag)
341  DO 260 i=kk,n
342  p2 = s2
343  s2 = ck*s2 + s1
344  s1 = p2
345  ck = ck + rz
346  p2 = s2*p1
347  y(i) = p2
348  IF (kflag.GE.3) go to 260
349  p2r = REAL(p2)
350  p2i = aimag(p2)
351  p2r = abs(p2r)
352  p2i = abs(p2i)
353  p2m = amax1(p2r,p2i)
354  IF (p2m.LE.ascle) go to 260
355  kflag = kflag + 1
356  ascle = bry(kflag)
357  s1 = s1*p1
358  s2 = p2
359  s1 = s1*css(kflag)
360  s2 = s2*css(kflag)
361  p1 = csr(kflag)
362  260 CONTINUE
363  RETURN
364 C-----------------------------------------------------------------------
365 C IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
366 C-----------------------------------------------------------------------
367  261 CONTINUE
368  helim = 0.5e0*elim
369  elm = exp(-elim)
370  celm = cmplx(elm,0.0)
371  ascle = bry(1)
372  zd = z
373  xd = xx
374  yd = yy
375  ic = -1
376  j = 2
377  DO 262 i=1,inu
378  st = s2
379  s2 = ck*s2+s1
380  s1 = st
381  ck = ck+rz
382  as = cabs(s2)
383  alas = alog(as)
384  p2r = -xd+alas
385  IF(p2r.LT.(-elim)) go to 263
386  p2 = -zd+clog(s2)
387  p2r = REAL(p2)
388  p2i = aimag(p2)
389  p2m = exp(p2r)/tol
390  p1 = cmplx(p2m,0.0e0)*cmplx(cos(p2i),sin(p2i))
391  CALL cuchk(p1,nw,ascle,tol)
392  IF(nw.NE.0) go to 263
393  j=3-j
394  cy(j) = p1
395  IF(ic.EQ.(i-1)) go to 264
396  ic = i
397  go to 262
398  263 CONTINUE
399  IF(alas.LT.helim) go to 262
400  xd = xd-elim
401  s1 = s1*celm
402  s2 = s2*celm
403  zd = cmplx(xd,yd)
404  262 CONTINUE
405  IF(n.EQ.1) s1 = s2
406  go to 270
407  264 CONTINUE
408  kflag = 1
409  inub = i+1
410  s2 = cy(j)
411  j = 3 - j
412  s1 = cy(j)
413  IF(inub.LE.inu) go to 225
414  IF(n.EQ.1) s1 = s2
415  go to 240
416  270 CONTINUE
417  y(1) = s1
418  IF (n.EQ.1) go to 280
419  y(2) = s2
420  280 CONTINUE
421  ascle = bry(1)
422  CALL ckscl(zd, fnu, n, y, nz, rz, ascle, tol, elim)
423  inu = n - nz
424  IF (inu.LE.0) RETURN
425  kk = nz + 1
426  s1 = y(kk)
427  y(kk) = s1*csr(1)
428  IF (inu.EQ.1) RETURN
429  kk = nz + 2
430  s2 = y(kk)
431  y(kk) = s2*csr(1)
432  IF (inu.EQ.2) RETURN
433  t2 = fnu + float(kk-1)
434  ck = cmplx(t2,0.0e0)*rz
435  kflag = 1
436  go to 250
437  290 CONTINUE
438 C-----------------------------------------------------------------------
439 C SCALE BY EXP(Z), IFLAG = 1 CASES
440 C-----------------------------------------------------------------------
441  koded = 2
442  iflag = 1
443  kflag = 2
444  go to 120
445 C-----------------------------------------------------------------------
446 C FNU=HALF ODD INTEGER CASE, DNU=-0.5
447 C-----------------------------------------------------------------------
448  300 CONTINUE
449  s1 = coef
450  s2 = coef
451  go to 210
452  310 CONTINUE
453  nz=-2
454  RETURN
455  END
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2