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
casyi.f
Go to the documentation of this file.
1  SUBROUTINE casyi(Z, FNU, KODE, N, Y, NZ, RL, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CASYI
3 C***REFER TO CBESI,CBESK
4 C
5 C CASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
6 C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
7 C REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
8 C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
9 C
10 C***ROUTINES CALLED R1MACH
11 C***END PROLOGUE CASYI
12  COMPLEX ak1, ck, cone, cs1, cs2, cz, czero, dk, ez, p1, rz, s2,
13  * y, z
14  REAL aa, acz, aez, ak, alim, arg, arm, atol, az, bb, bk, dfnu,
15  * dnu2, elim, fdn, fnu, pi, rl, rtpi, rtr1, s, sgn, sqk, tol, x,
16  * yy, r1mach
17  INTEGER i, ib, il, inu, j, jl, k, kode, koded, m, n, nn, nz
18  dimension y(n)
19  DATA pi, rtpi /3.14159265358979324e0 , 0.159154943091895336e0 /
20  DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
21 C
22  nz = 0
23  az = cabs(z)
24  x = REAL(z)
25  arm = 1.0e+3*r1mach(1)
26  rtr1 = sqrt(arm)
27  il = min0(2,n)
28  dfnu = fnu + float(n-il)
29 C-----------------------------------------------------------------------
30 C OVERFLOW TEST
31 C-----------------------------------------------------------------------
32  ak1 = cmplx(rtpi,0.0e0)/z
33  ak1 = csqrt(ak1)
34  cz = z
35  IF (kode.EQ.2) cz = z - cmplx(x,0.0e0)
36  acz = REAL(cz)
37  IF (abs(acz).GT.elim) go to 80
38  dnu2 = dfnu + dfnu
39  koded = 1
40  IF ((abs(acz).GT.alim) .AND. (n.GT.2)) go to 10
41  koded = 0
42  ak1 = ak1*cexp(cz)
43  10 CONTINUE
44  fdn = 0.0e0
45  IF (dnu2.GT.rtr1) fdn = dnu2*dnu2
46  ez = z*cmplx(8.0e0,0.0e0)
47 C-----------------------------------------------------------------------
48 C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
49 C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
50 C EXPANSION FOR THE IMAGINARY PART.
51 C-----------------------------------------------------------------------
52  aez = 8.0e0*az
53  s = tol/aez
54  jl = int(rl+rl) + 2
55  yy = aimag(z)
56  p1 = czero
57  IF (yy.EQ.0.0e0) go to 20
58 C-----------------------------------------------------------------------
59 C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
60 C SIGNIFICANCE WHEN FNU OR N IS LARGE
61 C-----------------------------------------------------------------------
62  inu = int(fnu)
63  arg = (fnu-float(inu))*pi
64  inu = inu + n - il
65  ak = -sin(arg)
66  bk = cos(arg)
67  IF (yy.LT.0.0e0) bk = -bk
68  p1 = cmplx(ak,bk)
69  IF (mod(inu,2).EQ.1) p1 = -p1
70  20 CONTINUE
71  DO 50 k=1,il
72  sqk = fdn - 1.0e0
73  atol = s*abs(sqk)
74  sgn = 1.0e0
75  cs1 = cone
76  cs2 = cone
77  ck = cone
78  ak = 0.0e0
79  aa = 1.0e0
80  bb = aez
81  dk = ez
82  DO 30 j=1,jl
83  ck = ck*cmplx(sqk,0.0e0)/dk
84  cs2 = cs2 + ck
85  sgn = -sgn
86  cs1 = cs1 + ck*cmplx(sgn,0.0e0)
87  dk = dk + ez
88  aa = aa*abs(sqk)/bb
89  bb = bb + aez
90  ak = ak + 8.0e0
91  sqk = sqk - ak
92  IF (aa.LE.atol) go to 40
93  30 CONTINUE
94  go to 90
95  40 CONTINUE
96  s2 = cs1
97  IF (x+x.LT.elim) s2 = s2 + p1*cs2*cexp(-z-z)
98  fdn = fdn + 8.0e0*dfnu + 4.0e0
99  p1 = -p1
100  m = n - il + k
101  y(m) = s2*ak1
102  50 CONTINUE
103  IF (n.LE.2) RETURN
104  nn = n
105  k = nn - 2
106  ak = float(k)
107  rz = (cone+cone)/z
108  ib = 3
109  DO 60 i=ib,nn
110  y(k) = cmplx(ak+fnu,0.0e0)*rz*y(k+1) + y(k+2)
111  ak = ak - 1.0e0
112  k = k - 1
113  60 CONTINUE
114  IF (koded.EQ.0) RETURN
115  ck = cexp(cz)
116  DO 70 i=1,nn
117  y(i) = y(i)*ck
118  70 CONTINUE
119  RETURN
120  80 CONTINUE
121  nz = -1
122  RETURN
123  90 CONTINUE
124  nz=-2
125  RETURN
126  END