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
zbknu.f
Go to the documentation of this file.
1  SUBROUTINE zbknu(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
2  * alim)
3 C***BEGIN PROLOGUE ZBKNU
4 C***REFER TO ZBESI,ZBESK,ZAIRY,ZBESH
5 C
6 C ZBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE.
7 C
8 C***ROUTINES CALLED DGAMLN,I1MACH,D1MACH,ZKSCL,ZSHCH,ZUCHK,XZABS,ZDIV,
9 C XZEXP,XZLOG,ZMLT,XZSQRT
10 C***END PROLOGUE ZBKNU
11 C
12  DOUBLE PRECISION aa, ak, alim, ascle, a1, a2, bb, bk, bry, caz,
13  * cbi, cbr, cc, cchi, cchr, cki, ckr, coefi, coefr, conei, coner,
14  * crscr, csclr, cshi, cshr, csi, csr, csrr, cssr, ctwor,
15  * czeroi, czeror, czi, czr, dnu, dnu2, dpi, elim, etest, fc, fhs,
16  * fi, fk, fks, fmui, fmur, fnu, fpi, fr, g1, g2, hpi, pi, pr, pti,
17  * ptr, p1i, p1r, p2i, p2m, p2r, qi, qr, rak, rcaz, rthpi, rzi,
18  * rzr, r1, s, smui, smur, spi, sti, str, s1i, s1r, s2i, s2r, tm,
19  * tol, tth, t1, t2, yi, yr, zi, zr, dgamln, d1mach, xzabs, elm,
20  * celmr, zdr, zdi, as, alas, helim, cyr, cyi
21  INTEGER i, iflag, inu, k, kflag, kk, kmax, kode, koded, n, nz,
22  * idum, i1mach, j, ic, inub, nw
23  dimension yr(n), yi(n), cc(8), cssr(3), csrr(3), bry(3), cyr(2),
24  * cyi(2)
25 C COMPLEX Z,Y,A,B,RZ,SMU,FU,FMU,F,FLRZ,CZ,S1,S2,CSH,CCH
26 C COMPLEX CK,P,Q,COEF,P1,P2,CBK,PT,CZERO,CONE,CTWO,ST,EZ,CS,DK
27 C
28  DATA kmax / 30 /
29  DATA czeror,czeroi,coner,conei,ctwor,r1/
30  1 0.0d0 , 0.0d0 , 1.0d0 , 0.0d0 , 2.0d0 , 2.0d0 /
31  DATA dpi, rthpi, spi ,hpi, fpi, tth /
32  1 3.14159265358979324d0, 1.25331413731550025d0,
33  2 1.90985931710274403d0, 1.57079632679489662d0,
34  3 1.89769999331517738d0, 6.66666666666666666d-01/
35  DATA cc(1), cc(2), cc(3), cc(4), cc(5), cc(6), cc(7), cc(8)/
36  1 5.77215664901532861d-01, -4.20026350340952355d-02,
37  2 -4.21977345555443367d-02, 7.21894324666309954d-03,
38  3 -2.15241674114950973d-04, -2.01348547807882387d-05,
39  4 1.13302723198169588d-06, 6.11609510448141582d-09/
40 C
41  caz = xzabs(zr,zi)
42  csclr = 1.0d0/tol
43  crscr = tol
44  cssr(1) = csclr
45  cssr(2) = 1.0d0
46  cssr(3) = crscr
47  csrr(1) = crscr
48  csrr(2) = 1.0d0
49  csrr(3) = csclr
50  bry(1) = 1.0d+3*d1mach(1)/tol
51  bry(2) = 1.0d0/bry(1)
52  bry(3) = d1mach(2)
53  nz = 0
54  iflag = 0
55  koded = kode
56  rcaz = 1.0d0/caz
57  str = zr*rcaz
58  sti = -zi*rcaz
59  rzr = (str+str)*rcaz
60  rzi = (sti+sti)*rcaz
61  inu = int(sngl(fnu+0.5d0))
62  dnu = fnu - dble(float(inu))
63  IF (dabs(dnu).EQ.0.5d0) go to 110
64  dnu2 = 0.0d0
65  IF (dabs(dnu).GT.tol) dnu2 = dnu*dnu
66  IF (caz.GT.r1) go to 110
67 C-----------------------------------------------------------------------
68 C SERIES FOR CABS(Z).LE.R1
69 C-----------------------------------------------------------------------
70  fc = 1.0d0
71  CALL xzlog(rzr, rzi, smur, smui, idum)
72  fmur = smur*dnu
73  fmui = smui*dnu
74  CALL zshch(fmur, fmui, cshr, cshi, cchr, cchi)
75  IF (dnu.EQ.0.0d0) go to 10
76  fc = dnu*dpi
77  fc = fc/dsin(fc)
78  smur = cshr/dnu
79  smui = cshi/dnu
80  10 CONTINUE
81  a2 = 1.0d0 + dnu
82 C-----------------------------------------------------------------------
83 C GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
84 C-----------------------------------------------------------------------
85  t2 = dexp(-dgamln(a2,idum))
86  t1 = 1.0d0/(t2*fc)
87  IF (dabs(dnu).GT.0.1d0) go to 40
88 C-----------------------------------------------------------------------
89 C SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
90 C-----------------------------------------------------------------------
91  ak = 1.0d0
92  s = cc(1)
93  DO 20 k=2,8
94  ak = ak*dnu2
95  tm = cc(k)*ak
96  s = s + tm
97  IF (dabs(tm).LT.tol) go to 30
98  20 CONTINUE
99  30 g1 = -s
100  go to 50
101  40 CONTINUE
102  g1 = (t1-t2)/(dnu+dnu)
103  50 CONTINUE
104  g2 = (t1+t2)*0.5d0
105  fr = fc*(cchr*g1+smur*g2)
106  fi = fc*(cchi*g1+smui*g2)
107  CALL xzexp(fmur, fmui, str, sti)
108  pr = 0.5d0*str/t2
109  pi = 0.5d0*sti/t2
110  CALL zdiv(0.5d0, 0.0d0, str, sti, ptr, pti)
111  qr = ptr/t1
112  qi = pti/t1
113  s1r = fr
114  s1i = fi
115  s2r = pr
116  s2i = pi
117  ak = 1.0d0
118  a1 = 1.0d0
119  ckr = coner
120  cki = conei
121  bk = 1.0d0 - dnu2
122  IF (inu.GT.0 .OR. n.GT.1) go to 80
123 C-----------------------------------------------------------------------
124 C GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
125 C-----------------------------------------------------------------------
126  IF (caz.LT.tol) go to 70
127  CALL zmlt(zr, zi, zr, zi, czr, czi)
128  czr = 0.25d0*czr
129  czi = 0.25d0*czi
130  t1 = 0.25d0*caz*caz
131  60 CONTINUE
132  fr = (fr*ak+pr+qr)/bk
133  fi = (fi*ak+pi+qi)/bk
134  str = 1.0d0/(ak-dnu)
135  pr = pr*str
136  pi = pi*str
137  str = 1.0d0/(ak+dnu)
138  qr = qr*str
139  qi = qi*str
140  str = ckr*czr - cki*czi
141  rak = 1.0d0/ak
142  cki = (ckr*czi+cki*czr)*rak
143  ckr = str*rak
144  s1r = ckr*fr - cki*fi + s1r
145  s1i = ckr*fi + cki*fr + s1i
146  a1 = a1*t1*rak
147  bk = bk + ak + ak + 1.0d0
148  ak = ak + 1.0d0
149  IF (a1.GT.tol) go to 60
150  70 CONTINUE
151  yr(1) = s1r
152  yi(1) = s1i
153  IF (koded.EQ.1) RETURN
154  CALL xzexp(zr, zi, str, sti)
155  CALL zmlt(s1r, s1i, str, sti, yr(1), yi(1))
156  RETURN
157 C-----------------------------------------------------------------------
158 C GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
159 C-----------------------------------------------------------------------
160  80 CONTINUE
161  IF (caz.LT.tol) go to 100
162  CALL zmlt(zr, zi, zr, zi, czr, czi)
163  czr = 0.25d0*czr
164  czi = 0.25d0*czi
165  t1 = 0.25d0*caz*caz
166  90 CONTINUE
167  fr = (fr*ak+pr+qr)/bk
168  fi = (fi*ak+pi+qi)/bk
169  str = 1.0d0/(ak-dnu)
170  pr = pr*str
171  pi = pi*str
172  str = 1.0d0/(ak+dnu)
173  qr = qr*str
174  qi = qi*str
175  str = ckr*czr - cki*czi
176  rak = 1.0d0/ak
177  cki = (ckr*czi+cki*czr)*rak
178  ckr = str*rak
179  s1r = ckr*fr - cki*fi + s1r
180  s1i = ckr*fi + cki*fr + s1i
181  str = pr - fr*ak
182  sti = pi - fi*ak
183  s2r = ckr*str - cki*sti + s2r
184  s2i = ckr*sti + cki*str + s2i
185  a1 = a1*t1*rak
186  bk = bk + ak + ak + 1.0d0
187  ak = ak + 1.0d0
188  IF (a1.GT.tol) go to 90
189  100 CONTINUE
190  kflag = 2
191  a1 = fnu + 1.0d0
192  ak = a1*dabs(smur)
193  IF (ak.GT.alim) kflag = 3
194  str = cssr(kflag)
195  p2r = s2r*str
196  p2i = s2i*str
197  CALL zmlt(p2r, p2i, rzr, rzi, s2r, s2i)
198  s1r = s1r*str
199  s1i = s1i*str
200  IF (koded.EQ.1) go to 210
201  CALL xzexp(zr, zi, fr, fi)
202  CALL zmlt(s1r, s1i, fr, fi, s1r, s1i)
203  CALL zmlt(s2r, s2i, fr, fi, s2r, s2i)
204  go to 210
205 C-----------------------------------------------------------------------
206 C IFLAG=0 MEANS NO UNDERFLOW OCCURRED
207 C IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
208 C KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
209 C RECURSION
210 C-----------------------------------------------------------------------
211  110 CONTINUE
212  CALL xzsqrt(zr, zi, str, sti)
213  CALL zdiv(rthpi, czeroi, str, sti, coefr, coefi)
214  kflag = 2
215  IF (koded.EQ.2) go to 120
216  IF (zr.GT.alim) go to 290
217 C BLANK LINE
218  str = dexp(-zr)*cssr(kflag)
219  sti = -str*dsin(zi)
220  str = str*dcos(zi)
221  CALL zmlt(coefr, coefi, str, sti, coefr, coefi)
222  120 CONTINUE
223  IF (dabs(dnu).EQ.0.5d0) go to 300
224 C-----------------------------------------------------------------------
225 C MILLER ALGORITHM FOR CABS(Z).GT.R1
226 C-----------------------------------------------------------------------
227  ak = dcos(dpi*dnu)
228  ak = dabs(ak)
229  IF (ak.EQ.czeror) go to 300
230  fhs = dabs(0.25d0-dnu2)
231  IF (fhs.EQ.czeror) go to 300
232 C-----------------------------------------------------------------------
233 C COMPUTE R2=F(E). IF CABS(Z).GE.R2, USE FORWARD RECURRENCE TO
234 C DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
235 C 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))=
236 C TOL WHERE B IS THE BASE OF THE ARITHMETIC.
237 C-----------------------------------------------------------------------
238  t1 = dble(float(i1mach(14)-1))
239  t1 = t1*d1mach(5)*3.321928094d0
240  t1 = dmax1(t1,12.0d0)
241  t1 = dmin1(t1,60.0d0)
242  t2 = tth*t1 - 6.0d0
243  IF (zr.NE.0.0d0) go to 130
244  t1 = hpi
245  go to 140
246  130 CONTINUE
247  t1 = datan(zi/zr)
248  t1 = dabs(t1)
249  140 CONTINUE
250  IF (t2.GT.caz) go to 170
251 C-----------------------------------------------------------------------
252 C FORWARD RECURRENCE LOOP WHEN CABS(Z).GE.R2
253 C-----------------------------------------------------------------------
254  etest = ak/(dpi*caz*tol)
255  fk = coner
256  IF (etest.LT.coner) go to 180
257  fks = ctwor
258  ckr = caz + caz + ctwor
259  p1r = czeror
260  p2r = coner
261  DO 150 i=1,kmax
262  ak = fhs/fks
263  cbr = ckr/(fk+coner)
264  ptr = p2r
265  p2r = cbr*p2r - p1r*ak
266  p1r = ptr
267  ckr = ckr + ctwor
268  fks = fks + fk + fk + ctwor
269  fhs = fhs + fk + fk
270  fk = fk + coner
271  str = dabs(p2r)*fk
272  IF (etest.LT.str) go to 160
273  150 CONTINUE
274  go to 310
275  160 CONTINUE
276  fk = fk + spi*t1*dsqrt(t2/caz)
277  fhs = dabs(0.25d0-dnu2)
278  go to 180
279  170 CONTINUE
280 C-----------------------------------------------------------------------
281 C COMPUTE BACKWARD INDEX K FOR CABS(Z).LT.R2
282 C-----------------------------------------------------------------------
283  a2 = dsqrt(caz)
284  ak = fpi*ak/(tol*dsqrt(a2))
285  aa = 3.0d0*t1/(1.0d0+caz)
286  bb = 14.7d0*t1/(28.0d0+caz)
287  ak = (dlog(ak)+caz*dcos(aa)/(1.0d0+0.008d0*caz))/dcos(bb)
288  fk = 0.12125d0*ak*ak/caz + 1.5d0
289  180 CONTINUE
290 C-----------------------------------------------------------------------
291 C BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
292 C-----------------------------------------------------------------------
293  k = int(sngl(fk))
294  fk = dble(float(k))
295  fks = fk*fk
296  p1r = czeror
297  p1i = czeroi
298  p2r = tol
299  p2i = czeroi
300  csr = p2r
301  csi = p2i
302  DO 190 i=1,k
303  a1 = fks - fk
304  ak = (fks+fk)/(a1+fhs)
305  rak = 2.0d0/(fk+coner)
306  cbr = (fk+zr)*rak
307  cbi = zi*rak
308  ptr = p2r
309  pti = p2i
310  p2r = (ptr*cbr-pti*cbi-p1r)*ak
311  p2i = (pti*cbr+ptr*cbi-p1i)*ak
312  p1r = ptr
313  p1i = pti
314  csr = csr + p2r
315  csi = csi + p2i
316  fks = a1 - fk + coner
317  fk = fk - coner
318  190 CONTINUE
319 C-----------------------------------------------------------------------
320 C COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER
321 C SCALING
322 C-----------------------------------------------------------------------
323  tm = xzabs(csr,csi)
324  ptr = 1.0d0/tm
325  s1r = p2r*ptr
326  s1i = p2i*ptr
327  csr = csr*ptr
328  csi = -csi*ptr
329  CALL zmlt(coefr, coefi, s1r, s1i, str, sti)
330  CALL zmlt(str, sti, csr, csi, s1r, s1i)
331  IF (inu.GT.0 .OR. n.GT.1) go to 200
332  zdr = zr
333  zdi = zi
334  IF(iflag.EQ.1) go to 270
335  go to 240
336  200 CONTINUE
337 C-----------------------------------------------------------------------
338 C COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING
339 C-----------------------------------------------------------------------
340  tm = xzabs(p2r,p2i)
341  ptr = 1.0d0/tm
342  p1r = p1r*ptr
343  p1i = p1i*ptr
344  p2r = p2r*ptr
345  p2i = -p2i*ptr
346  CALL zmlt(p1r, p1i, p2r, p2i, ptr, pti)
347  str = dnu + 0.5d0 - ptr
348  sti = -pti
349  CALL zdiv(str, sti, zr, zi, str, sti)
350  str = str + 1.0d0
351  CALL zmlt(str, sti, s1r, s1i, s2r, s2i)
352 C-----------------------------------------------------------------------
353 C FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH
354 C SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
355 C-----------------------------------------------------------------------
356  210 CONTINUE
357  str = dnu + 1.0d0
358  ckr = str*rzr
359  cki = str*rzi
360  IF (n.EQ.1) inu = inu - 1
361  IF (inu.GT.0) go to 220
362  IF (n.GT.1) go to 215
363  s1r = s2r
364  s1i = s2i
365  215 CONTINUE
366  zdr = zr
367  zdi = zi
368  IF(iflag.EQ.1) go to 270
369  go to 240
370  220 CONTINUE
371  inub = 1
372  IF(iflag.EQ.1) go to 261
373  225 CONTINUE
374  p1r = csrr(kflag)
375  ascle = bry(kflag)
376  DO 230 i=inub,inu
377  str = s2r
378  sti = s2i
379  s2r = ckr*str - cki*sti + s1r
380  s2i = ckr*sti + cki*str + s1i
381  s1r = str
382  s1i = sti
383  ckr = ckr + rzr
384  cki = cki + rzi
385  IF (kflag.GE.3) go to 230
386  p2r = s2r*p1r
387  p2i = s2i*p1r
388  str = dabs(p2r)
389  sti = dabs(p2i)
390  p2m = dmax1(str,sti)
391  IF (p2m.LE.ascle) go to 230
392  kflag = kflag + 1
393  ascle = bry(kflag)
394  s1r = s1r*p1r
395  s1i = s1i*p1r
396  s2r = p2r
397  s2i = p2i
398  str = cssr(kflag)
399  s1r = s1r*str
400  s1i = s1i*str
401  s2r = s2r*str
402  s2i = s2i*str
403  p1r = csrr(kflag)
404  230 CONTINUE
405  IF (n.NE.1) go to 240
406  s1r = s2r
407  s1i = s2i
408  240 CONTINUE
409  str = csrr(kflag)
410  yr(1) = s1r*str
411  yi(1) = s1i*str
412  IF (n.EQ.1) RETURN
413  yr(2) = s2r*str
414  yi(2) = s2i*str
415  IF (n.EQ.2) RETURN
416  kk = 2
417  250 CONTINUE
418  kk = kk + 1
419  IF (kk.GT.n) RETURN
420  p1r = csrr(kflag)
421  ascle = bry(kflag)
422  DO 260 i=kk,n
423  p2r = s2r
424  p2i = s2i
425  s2r = ckr*p2r - cki*p2i + s1r
426  s2i = cki*p2r + ckr*p2i + s1i
427  s1r = p2r
428  s1i = p2i
429  ckr = ckr + rzr
430  cki = cki + rzi
431  p2r = s2r*p1r
432  p2i = s2i*p1r
433  yr(i) = p2r
434  yi(i) = p2i
435  IF (kflag.GE.3) go to 260
436  str = dabs(p2r)
437  sti = dabs(p2i)
438  p2m = dmax1(str,sti)
439  IF (p2m.LE.ascle) go to 260
440  kflag = kflag + 1
441  ascle = bry(kflag)
442  s1r = s1r*p1r
443  s1i = s1i*p1r
444  s2r = p2r
445  s2i = p2i
446  str = cssr(kflag)
447  s1r = s1r*str
448  s1i = s1i*str
449  s2r = s2r*str
450  s2i = s2i*str
451  p1r = csrr(kflag)
452  260 CONTINUE
453  RETURN
454 C-----------------------------------------------------------------------
455 C IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
456 C-----------------------------------------------------------------------
457  261 CONTINUE
458  helim = 0.5d0*elim
459  elm = dexp(-elim)
460  celmr = elm
461  ascle = bry(1)
462  zdr = zr
463  zdi = zi
464  ic = -1
465  j = 2
466  DO 262 i=1,inu
467  str = s2r
468  sti = s2i
469  s2r = str*ckr-sti*cki+s1r
470  s2i = sti*ckr+str*cki+s1i
471  s1r = str
472  s1i = sti
473  ckr = ckr+rzr
474  cki = cki+rzi
475  as = xzabs(s2r,s2i)
476  alas = dlog(as)
477  p2r = -zdr+alas
478  IF(p2r.LT.(-elim)) go to 263
479  CALL xzlog(s2r,s2i,str,sti,idum)
480  p2r = -zdr+str
481  p2i = -zdi+sti
482  p2m = dexp(p2r)/tol
483  p1r = p2m*dcos(p2i)
484  p1i = p2m*dsin(p2i)
485  CALL zuchk(p1r,p1i,nw,ascle,tol)
486  IF(nw.NE.0) go to 263
487  j = 3 - j
488  cyr(j) = p1r
489  cyi(j) = p1i
490  IF(ic.EQ.(i-1)) go to 264
491  ic = i
492  go to 262
493  263 CONTINUE
494  IF(alas.LT.helim) go to 262
495  zdr = zdr-elim
496  s1r = s1r*celmr
497  s1i = s1i*celmr
498  s2r = s2r*celmr
499  s2i = s2i*celmr
500  262 CONTINUE
501  IF(n.NE.1) go to 270
502  s1r = s2r
503  s1i = s2i
504  go to 270
505  264 CONTINUE
506  kflag = 1
507  inub = i+1
508  s2r = cyr(j)
509  s2i = cyi(j)
510  j = 3 - j
511  s1r = cyr(j)
512  s1i = cyi(j)
513  IF(inub.LE.inu) go to 225
514  IF(n.NE.1) go to 240
515  s1r = s2r
516  s1i = s2i
517  go to 240
518  270 CONTINUE
519  yr(1) = s1r
520  yi(1) = s1i
521  IF(n.EQ.1) go to 280
522  yr(2) = s2r
523  yi(2) = s2i
524  280 CONTINUE
525  ascle = bry(1)
526  CALL zkscl(zdr,zdi,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim)
527  inu = n - nz
528  IF (inu.LE.0) RETURN
529  kk = nz + 1
530  s1r = yr(kk)
531  s1i = yi(kk)
532  yr(kk) = s1r*csrr(1)
533  yi(kk) = s1i*csrr(1)
534  IF (inu.EQ.1) RETURN
535  kk = nz + 2
536  s2r = yr(kk)
537  s2i = yi(kk)
538  yr(kk) = s2r*csrr(1)
539  yi(kk) = s2i*csrr(1)
540  IF (inu.EQ.2) RETURN
541  t2 = fnu + dble(float(kk-1))
542  ckr = t2*rzr
543  cki = t2*rzi
544  kflag = 1
545  go to 250
546  290 CONTINUE
547 C-----------------------------------------------------------------------
548 C SCALE BY DEXP(Z), IFLAG = 1 CASES
549 C-----------------------------------------------------------------------
550  koded = 2
551  iflag = 1
552  kflag = 2
553  go to 120
554 C-----------------------------------------------------------------------
555 C FNU=HALF ODD INTEGER CASE, DNU=-0.5
556 C-----------------------------------------------------------------------
557  300 CONTINUE
558  s1r = coefr
559  s1i = coefi
560  s2r = coefr
561  s2i = coefi
562  go to 210
563 C
564 C
565  310 CONTINUE
566  nz=-2
567  RETURN
568  END
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2