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
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