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
cacon.f
Go to the documentation of this file.
1  SUBROUTINE cacon(Z, FNU, KODE, MR, N, Y, NZ, RL, FNUL, TOL, ELIM,
2  * alim)
3 C***BEGIN PROLOGUE CACON
4 C***REFER TO CBESK,CBESH
5 C
6 C CACON APPLIES THE ANALYTIC CONTINUATION FORMULA
7 C
8 C K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
9 C MP=PI*MR*CMPLX(0.0,1.0)
10 C
11 C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
12 C HALF Z PLANE
13 C
14 C***ROUTINES CALLED CBINU,CBKNU,CS1S2,R1MACH
15 C***END PROLOGUE CACON
16  COMPLEX ck, cone, cs, cscl, cscr, csgn, cspn, css, csr, c1, c2,
17  * rz, sc1, sc2, st, s1, s2, y, z, zn, cy
18  REAL alim, arg, ascle, as2, bscle, bry, cpn, c1i, c1m, c1r, elim,
19  * fmr, fnu, fnul, pi, rl, sgn, spn, tol, yy, r1mach
20  INTEGER i, inu, iuf, kflag, kode, mr, n, nn, nw, nz
21  dimension y(n), cy(2), css(3), csr(3), bry(3)
22  DATA pi / 3.14159265358979324e0 /
23  DATA cone / (1.0e0,0.0e0) /
24  nz = 0
25  zn = -z
26  nn = n
27  CALL cbinu(zn, fnu, kode, nn, y, nw, rl, fnul, tol, elim, alim)
28  IF (nw.LT.0) go to 80
29 C-----------------------------------------------------------------------
30 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
31 C-----------------------------------------------------------------------
32  nn = min0(2,n)
33  CALL cbknu(zn, fnu, kode, nn, cy, nw, tol, elim, alim)
34  IF (nw.NE.0) go to 80
35  s1 = cy(1)
36  fmr = float(mr)
37  sgn = -sign(pi,fmr)
38  csgn = cmplx(0.0e0,sgn)
39  IF (kode.EQ.1) go to 10
40  yy = -aimag(zn)
41  cpn = cos(yy)
42  spn = sin(yy)
43  csgn = csgn*cmplx(cpn,spn)
44  10 CONTINUE
45 C-----------------------------------------------------------------------
46 C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
47 C WHEN FNU IS LARGE
48 C-----------------------------------------------------------------------
49  inu = int(fnu)
50  arg = (fnu-float(inu))*sgn
51  cpn = cos(arg)
52  spn = sin(arg)
53  cspn = cmplx(cpn,spn)
54  IF (mod(inu,2).EQ.1) cspn = -cspn
55  iuf = 0
56  c1 = s1
57  c2 = y(1)
58  ascle = 1.0e+3*r1mach(1)/tol
59  IF (kode.EQ.1) go to 20
60  CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
61  nz = nz + nw
62  sc1 = c1
63  20 CONTINUE
64  y(1) = cspn*c1 + csgn*c2
65  IF (n.EQ.1) RETURN
66  cspn = -cspn
67  s2 = cy(2)
68  c1 = s2
69  c2 = y(2)
70  IF (kode.EQ.1) go to 30
71  CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
72  nz = nz + nw
73  sc2 = c1
74  30 CONTINUE
75  y(2) = cspn*c1 + csgn*c2
76  IF (n.EQ.2) RETURN
77  cspn = -cspn
78  rz = cmplx(2.0e0,0.0e0)/zn
79  ck = cmplx(fnu+1.0e0,0.0e0)*rz
80 C-----------------------------------------------------------------------
81 C SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
82 C-----------------------------------------------------------------------
83  cscl = cmplx(1.0e0/tol,0.0e0)
84  cscr = cmplx(tol,0.0e0)
85  css(1) = cscl
86  css(2) = cone
87  css(3) = cscr
88  csr(1) = cscr
89  csr(2) = cone
90  csr(3) = cscl
91  bry(1) = ascle
92  bry(2) = 1.0e0/ascle
93  bry(3) = r1mach(2)
94  as2 = cabs(s2)
95  kflag = 2
96  IF (as2.GT.bry(1)) go to 40
97  kflag = 1
98  go to 50
99  40 CONTINUE
100  IF (as2.LT.bry(2)) go to 50
101  kflag = 3
102  50 CONTINUE
103  bscle = bry(kflag)
104  s1 = s1*css(kflag)
105  s2 = s2*css(kflag)
106  cs = csr(kflag)
107  DO 70 i=3,n
108  st = s2
109  s2 = ck*s2 + s1
110  s1 = st
111  c1 = s2*cs
112  st = c1
113  c2 = y(i)
114  IF (kode.EQ.1) go to 60
115  IF (iuf.LT.0) go to 60
116  CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
117  nz = nz + nw
118  sc1 = sc2
119  sc2 = c1
120  IF (iuf.NE.3) go to 60
121  iuf = -4
122  s1 = sc1*css(kflag)
123  s2 = sc2*css(kflag)
124  st = sc2
125  60 CONTINUE
126  y(i) = cspn*c1 + csgn*c2
127  ck = ck + rz
128  cspn = -cspn
129  IF (kflag.GE.3) go to 70
130  c1r = REAL(c1)
131  c1i = aimag(c1)
132  c1r = abs(c1r)
133  c1i = abs(c1i)
134  c1m = amax1(c1r,c1i)
135  IF (c1m.LE.bscle) go to 70
136  kflag = kflag + 1
137  bscle = bry(kflag)
138  s1 = s1*cs
139  s2 = st
140  s1 = s1*css(kflag)
141  s2 = s2*css(kflag)
142  cs = csr(kflag)
143  70 CONTINUE
144  RETURN
145  80 CONTINUE
146  nz = -1
147  IF(nw.EQ.(-2)) nz=-2
148  RETURN
149  END