GNU Octave  4.4.1
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
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:860
double precision function d1mach(i)
Definition: d1mach.f:20
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)