GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
cmplx
double complex cmplx
Definition: Faddeeva.cc:230
cbknu
subroutine cbknu(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
Definition: cbknu.f:2
ckscl
subroutine ckscl(ZR, FNU, N, Y, NZ, RZ, ASCLE, TOL, ELIM)
Definition: ckscl.f:2
cshch
subroutine cshch(Z, CSH, CCH)
Definition: cshch.f:2
cuchk
subroutine cuchk(Y, NZ, ASCLE, TOL)
Definition: cuchk.f:2
real
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
atan
Complex atan(const Complex &x)
Definition: lo-mappers.h:71