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
zacai.f
Go to the documentation of this file.
1  SUBROUTINE zacai(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL,
2  * elim, alim)
3 C***BEGIN PROLOGUE ZACAI
4 C***REFER TO ZAIRY
5 C
6 C ZACAI 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 FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
13 C ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND
14 C RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT IF ZACON
15 C IS CALLED FROM ZAIRY.
16 C
17 C***ROUTINES CALLED ZASYI,ZBKNU,ZMLRI,ZSERI,ZS1S2,D1MACH,XZABS
18 C***END PROLOGUE ZACAI
19 C COMPLEX CSGN,CSPN,C1,C2,Y,Z,ZN,CY
20  DOUBLE PRECISION alim, arg, ascle, az, csgnr, csgni, cspnr,
21  * cspni, c1r, c1i, c2r, c2i, cyr, cyi, dfnu, elim, fmr, fnu, pi,
22  * rl, sgn, tol, yy, yr, yi, zr, zi, znr, zni, d1mach, xzabs
23  INTEGER inu, iuf, kode, mr, n, nn, nw, nz
24  dimension yr(n), yi(n), cyr(2), cyi(2)
25  DATA pi / 3.14159265358979324d0 /
26  nz = 0
27  znr = -zr
28  zni = -zi
29  az = xzabs(zr,zi)
30  nn = n
31  dfnu = fnu + dble(float(n-1))
32  IF (az.LE.2.0d0) go to 10
33  IF (az*az*0.25d0.GT.dfnu+1.0d0) go to 20
34  10 CONTINUE
35 C-----------------------------------------------------------------------
36 C POWER SERIES FOR THE I FUNCTION
37 C-----------------------------------------------------------------------
38  CALL zseri(znr, zni, fnu, kode, nn, yr, yi, nw, tol, elim, alim)
39  go to 40
40  20 CONTINUE
41  IF (az.LT.rl) go to 30
42 C-----------------------------------------------------------------------
43 C ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
44 C-----------------------------------------------------------------------
45  CALL zasyi(znr, zni, fnu, kode, nn, yr, yi, nw, rl, tol, elim,
46  * alim)
47  IF (nw.LT.0) go to 80
48  go to 40
49  30 CONTINUE
50 C-----------------------------------------------------------------------
51 C MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
52 C-----------------------------------------------------------------------
53  CALL zmlri(znr, zni, fnu, kode, nn, yr, yi, nw, tol)
54  IF(nw.LT.0) go to 80
55  40 CONTINUE
56 C-----------------------------------------------------------------------
57 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
58 C-----------------------------------------------------------------------
59  CALL zbknu(znr, zni, fnu, kode, 1, cyr, cyi, nw, tol, elim, alim)
60  IF (nw.NE.0) go to 80
61  fmr = dble(float(mr))
62  sgn = -dsign(pi,fmr)
63  csgnr = 0.0d0
64  csgni = sgn
65  IF (kode.EQ.1) go to 50
66  yy = -zni
67  csgnr = -csgni*dsin(yy)
68  csgni = csgni*dcos(yy)
69  50 CONTINUE
70 C-----------------------------------------------------------------------
71 C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
72 C WHEN FNU IS LARGE
73 C-----------------------------------------------------------------------
74  inu = int(sngl(fnu))
75  arg = (fnu-dble(float(inu)))*sgn
76  cspnr = dcos(arg)
77  cspni = dsin(arg)
78  IF (mod(inu,2).EQ.0) go to 60
79  cspnr = -cspnr
80  cspni = -cspni
81  60 CONTINUE
82  c1r = cyr(1)
83  c1i = cyi(1)
84  c2r = yr(1)
85  c2i = yi(1)
86  IF (kode.EQ.1) go to 70
87  iuf = 0
88  ascle = 1.0d+3*d1mach(1)/tol
89  CALL zs1s2(znr, zni, c1r, c1i, c2r, c2i, nw, ascle, alim, iuf)
90  nz = nz + nw
91  70 CONTINUE
92  yr(1) = cspnr*c1r - cspni*c1i + csgnr*c2r - csgni*c2i
93  yi(1) = cspnr*c1i + cspni*c1r + csgnr*c2i + csgni*c2r
94  RETURN
95  80 CONTINUE
96  nz = -1
97  IF(nw.EQ.(-2)) nz=-2
98  RETURN
99  END