GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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