GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cacai.f
Go to the documentation of this file.
1  SUBROUTINE cacai(Z, FNU, KODE, MR, N, Y, NZ, RL, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CACAI
3 C***REFER TO CAIRY
4 C
5 C CACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
6 C
7 C K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
8 C MP=PI*MR*CMPLX(0.0,1.0)
9 C
10 C TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
11 C HALF Z PLANE FOR USE WITH CAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
12 C CACAI IS THE SAME AS CACON WITH THE PARTS FOR LARGER ORDERS AND
13 C RECURRENCE REMOVED. A RECURSIVE CALL TO CACON CAN RESULT IF CACON
14 C IS CALLED FROM CAIRY.
15 C
16 C***ROUTINES CALLED CASYI,CBKNU,CMLRI,CSERI,CS1S2,R1MACH
17 C***END PROLOGUE CACAI
18  COMPLEX CSGN, CSPN, C1, C2, Y, Z, ZN, CY
19  REAL ALIM, ARG, ASCLE, AZ, CPN, DFNU, ELIM, FMR, FNU, PI, RL,
20  * SGN, SPN, TOL, YY, R1MACH
21  INTEGER INU, IUF, KODE, MR, N, NN, NW, NZ
22  dimension y(n), cy(2)
23  DATA pi / 3.14159265358979324e0 /
24  nz = 0
25  zn = -z
26  az = cabs(z)
27  nn = n
28  dfnu = fnu + float(n-1)
29  IF (az.LE.2.0e0) GO TO 10
30  IF (az*az*0.25e0.GT.dfnu+1.0e0) GO TO 20
31  10 CONTINUE
32 C-----------------------------------------------------------------------
33 C POWER SERIES FOR THE I FUNCTION
34 C-----------------------------------------------------------------------
35  CALL cseri(zn, fnu, kode, nn, y, nw, tol, elim, alim)
36  GO TO 40
37  20 CONTINUE
38  IF (az.LT.rl) GO TO 30
39 C-----------------------------------------------------------------------
40 C ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
41 C-----------------------------------------------------------------------
42  CALL casyi(zn, fnu, kode, nn, y, nw, rl, tol, elim, alim)
43  IF (nw.LT.0) GO TO 70
44  GO TO 40
45  30 CONTINUE
46 C-----------------------------------------------------------------------
47 C MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
48 C-----------------------------------------------------------------------
49  CALL cmlri(zn, fnu, kode, nn, y, nw, tol)
50  IF(nw.LT.0) GO TO 70
51  40 CONTINUE
52 C-----------------------------------------------------------------------
53 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
54 C-----------------------------------------------------------------------
55  CALL cbknu(zn, fnu, kode, 1, cy, nw, tol, elim, alim)
56  IF (nw.NE.0) GO TO 70
57  fmr = float(mr)
58  sgn = -sign(pi,fmr)
59  csgn = cmplx(0.0e0,sgn)
60  IF (kode.EQ.1) GO TO 50
61  yy = -aimag(zn)
62  cpn = cos(yy)
63  spn = sin(yy)
64  csgn = csgn*cmplx(cpn,spn)
65  50 CONTINUE
66 C-----------------------------------------------------------------------
67 C CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
68 C WHEN FNU IS LARGE
69 C-----------------------------------------------------------------------
70  inu = int(fnu)
71  arg = (fnu-float(inu))*sgn
72  cpn = cos(arg)
73  spn = sin(arg)
74  cspn = cmplx(cpn,spn)
75  IF (mod(inu,2).EQ.1) cspn = -cspn
76  c1 = cy(1)
77  c2 = y(1)
78  IF (kode.EQ.1) GO TO 60
79  iuf = 0
80  ascle = 1.0e+3*r1mach(1)/tol
81  CALL cs1s2(zn, c1, c2, nw, ascle, alim, iuf)
82  nz = nz + nw
83  60 CONTINUE
84  y(1) = cspn*c1 + csgn*c2
85  RETURN
86  70 CONTINUE
87  nz = -1
88  IF(nw.EQ.(-2)) nz=-2
89  RETURN
90  END
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:860
double complex cmplx
Definition: Faddeeva.cc:217
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)