GNU Octave  4.2.1
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
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:977
s
Definition: file-io.cc:2682
i e
Definition: data.cc:2724
double complex cmplx
Definition: Faddeeva.cc:217
octave_value sin(void) const
Definition: ov.h:1386
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
real function r1mach(i)
Definition: r1mach.f:1
octave_value cos(void) const
Definition: ov.h:1358
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
subroutine casyi(Z, FNU, KODE, N, Y, NZ, RL, TOL, ELIM, ALIM)
Definition: casyi.f:1
octave_value sqrt(void) const
Definition: ov.h:1388
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x