GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
subroutine zbinu(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: zbinu.f:1
std::string str(char sep= 'x') const
Definition: dim-vector.cc:63
std::string dimension(void) const
subroutine zacon(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: zacon.f:1
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:959
double precision function d1mach(i)
Definition: d1mach.f:1
F77_RET_T const double const double double * d
subroutine zs1s2(ZRR, ZRI, S1R, S1I, S2R, S2I, NZ, ASCLE, ALIM, IUF)
Definition: zs1s2.f:1
subroutine zmlt(AR, AI, BR, BI, CR, CI)
Definition: zmlt.f:1
subroutine zbknu(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM)
Definition: zbknu.f:1