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
zacon.f
Go to the documentation of this file.
1  SUBROUTINE zacon(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL,
2  * tol, elim, alim)
3 C***BEGIN PROLOGUE ZACON
4 C***REFER TO ZBESK,ZBESH
5 C
6 C ZACON 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 ZBINU,ZBKNU,ZS1S2,D1MACH,XZABS,ZMLT
15 C***END PROLOGUE ZACON
16 C COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST,
17 C *S1,S2,Y,Z,ZN
18  DOUBLE PRECISION alim, arg, ascle, as2, azn, bry, bscle, cki,
19  * ckr, coner, cpn, cscl, cscr, csgni, csgnr, cspni, cspnr,
20  * csr, csrr, cssr, cyi, cyr, c1i, c1m, c1r, c2i, c2r, elim, fmr,
21  * fn, fnu, fnul, pi, pti, ptr, razn, rl, rzi, rzr, sc1i, sc1r,
22  * sc2i, sc2r, sgn, spn, sti, str, s1i, s1r, s2i, s2r, tol, yi, yr,
23  * yy, zeror, zi, zni, znr, zr, d1mach, xzabs
24  INTEGER i, inu, iuf, kflag, kode, mr, n, nn, nw, nz
25  dimension yr(n), yi(n), cyr(2), cyi(2), cssr(3), csrr(3), bry(3)
26  DATA pi / 3.14159265358979324d0 /
27  DATA zeror,coner / 0.0d0,1.0d0 /
28  nz = 0
29  znr = -zr
30  zni = -zi
31  nn = n
32  CALL zbinu(znr, zni, fnu, kode, nn, yr, yi, nw, rl, fnul, tol,
33  * elim, alim)
34  IF (nw.LT.0) go to 90
35 C-----------------------------------------------------------------------
36 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
37 C-----------------------------------------------------------------------
38  nn = min0(2,n)
39  CALL zbknu(znr, zni, fnu, kode, nn, cyr, cyi, nw, tol, elim, alim)
40  IF (nw.NE.0) go to 90
41  s1r = cyr(1)
42  s1i = cyi(1)
43  fmr = dble(float(mr))
44  sgn = -dsign(pi,fmr)
45  csgnr = zeror
46  csgni = sgn
47  IF (kode.EQ.1) go to 10
48  yy = -zni
49  cpn = dcos(yy)
50  spn = dsin(yy)
51  CALL zmlt(csgnr, csgni, cpn, spn, csgnr, csgni)
52  10 CONTINUE
53 C-----------------------------------------------------------------------
54 C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
55 C WHEN FNU IS LARGE
56 C-----------------------------------------------------------------------
57  inu = int(sngl(fnu))
58  arg = (fnu-dble(float(inu)))*sgn
59  cpn = dcos(arg)
60  spn = dsin(arg)
61  cspnr = cpn
62  cspni = spn
63  IF (mod(inu,2).EQ.0) go to 20
64  cspnr = -cspnr
65  cspni = -cspni
66  20 CONTINUE
67  iuf = 0
68  c1r = s1r
69  c1i = s1i
70  c2r = yr(1)
71  c2i = yi(1)
72  ascle = 1.0d+3*d1mach(1)/tol
73  IF (kode.EQ.1) go to 30
74  CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
75  nz = nz + nw
76  sc1r = c1r
77  sc1i = c1i
78  30 CONTINUE
79  CALL zmlt(cspnr, cspni, c1r, c1i, str, sti)
80  CALL zmlt(csgnr, csgni, c2r, c2i, ptr, pti)
81  yr(1) = str + ptr
82  yi(1) = sti + pti
83  IF (n.EQ.1) RETURN
84  cspnr = -cspnr
85  cspni = -cspni
86  s2r = cyr(2)
87  s2i = cyi(2)
88  c1r = s2r
89  c1i = s2i
90  c2r = yr(2)
91  c2i = yi(2)
92  IF (kode.EQ.1) go to 40
93  CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
94  nz = nz + nw
95  sc2r = c1r
96  sc2i = c1i
97  40 CONTINUE
98  CALL zmlt(cspnr, cspni, c1r, c1i, str, sti)
99  CALL zmlt(csgnr, csgni, c2r, c2i, ptr, pti)
100  yr(2) = str + ptr
101  yi(2) = sti + pti
102  IF (n.EQ.2) RETURN
103  cspnr = -cspnr
104  cspni = -cspni
105  azn = xzabs(znr,zni)
106  razn = 1.0d0/azn
107  str = znr*razn
108  sti = -zni*razn
109  rzr = (str+str)*razn
110  rzi = (sti+sti)*razn
111  fn = fnu + 1.0d0
112  ckr = fn*rzr
113  cki = fn*rzi
114 C-----------------------------------------------------------------------
115 C SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
116 C-----------------------------------------------------------------------
117  cscl = 1.0d0/tol
118  cscr = tol
119  cssr(1) = cscl
120  cssr(2) = coner
121  cssr(3) = cscr
122  csrr(1) = cscr
123  csrr(2) = coner
124  csrr(3) = cscl
125  bry(1) = ascle
126  bry(2) = 1.0d0/ascle
127  bry(3) = d1mach(2)
128  as2 = xzabs(s2r,s2i)
129  kflag = 2
130  IF (as2.GT.bry(1)) go to 50
131  kflag = 1
132  go to 60
133  50 CONTINUE
134  IF (as2.LT.bry(2)) go to 60
135  kflag = 3
136  60 CONTINUE
137  bscle = bry(kflag)
138  s1r = s1r*cssr(kflag)
139  s1i = s1i*cssr(kflag)
140  s2r = s2r*cssr(kflag)
141  s2i = s2i*cssr(kflag)
142  csr = csrr(kflag)
143  DO 80 i=3,n
144  str = s2r
145  sti = s2i
146  s2r = ckr*str - cki*sti + s1r
147  s2i = ckr*sti + cki*str + s1i
148  s1r = str
149  s1i = sti
150  c1r = s2r*csr
151  c1i = s2i*csr
152  str = c1r
153  sti = c1i
154  c2r = yr(i)
155  c2i = yi(i)
156  IF (kode.EQ.1) go to 70
157  IF (iuf.LT.0) go to 70
158  CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
159  nz = nz + nw
160  sc1r = sc2r
161  sc1i = sc2i
162  sc2r = c1r
163  sc2i = c1i
164  IF (iuf.NE.3) go to 70
165  iuf = -4
166  s1r = sc1r*cssr(kflag)
167  s1i = sc1i*cssr(kflag)
168  s2r = sc2r*cssr(kflag)
169  s2i = sc2i*cssr(kflag)
170  str = sc2r
171  sti = sc2i
172  70 CONTINUE
173  ptr = cspnr*c1r - cspni*c1i
174  pti = cspnr*c1i + cspni*c1r
175  yr(i) = ptr + csgnr*c2r - csgni*c2i
176  yi(i) = pti + csgnr*c2i + csgni*c2r
177  ckr = ckr + rzr
178  cki = cki + rzi
179  cspnr = -cspnr
180  cspni = -cspni
181  IF (kflag.GE.3) go to 80
182  ptr = dabs(c1r)
183  pti = dabs(c1i)
184  c1m = dmax1(ptr,pti)
185  IF (c1m.LE.bscle) go to 80
186  kflag = kflag + 1
187  bscle = bry(kflag)
188  s1r = s1r*csr
189  s1i = s1i*csr
190  s2r = str
191  s2i = sti
192  s1r = s1r*cssr(kflag)
193  s1i = s1i*cssr(kflag)
194  s2r = s2r*cssr(kflag)
195  s2i = s2i*cssr(kflag)
196  csr = csrr(kflag)
197  80 CONTINUE
198  RETURN
199  90 CONTINUE
200  nz = -1
201  IF(nw.EQ.(-2)) nz=-2
202  RETURN
203  END