GNU Octave  4.2.1 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
subroutine cacon(Z, FNU, KODE, MR, N, Y, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: cacon.f:1
s2
Definition: sub2ind.cc:107
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:977
i e
Definition: data.cc:2724
double complex cmplx
octave_value sin(void) const
Definition: ov.h:1386
subroutine cs1s2(ZR, S1, S2, NZ, ASCLE, ALIM, IUF)
Definition: cs1s2.f:1
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to