GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
std::string dimension(void) const
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:959
double complex cmplx