GNU Octave  4.0.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
cseri.f
Go to the documentation of this file.
1  SUBROUTINE cseri(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CSERI
3 C***REFER TO CBESI,CBESK
4 C
5 C CSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
6 C MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE
7 C REGION CABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
8 C NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
9 C DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE
10 C CONDITION CABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE
11 C COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
12 C
13 C***ROUTINES CALLED CUCHK,GAMLN,R1MACH
14 C***END PROLOGUE CSERI
15  COMPLEX AK1, CK, COEF, CONE, CRSC, CZ, CZERO, HZ, RZ, S1, S2, W,
16  * y, z
17  REAL AA, ACZ, AK, ALIM, ARM, ASCLE, ATOL, AZ, DFNU, ELIM, FNU,
18  * fnup, rak1, rs, rtr1, s, ss, tol, x, gamln, r1mach
19  INTEGER I, IB, IDUM, IFLAG, IL, K, KODE, L, M, N, NN, NW, NZ
20  dimension y(n), w(2)
21  DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
22 C
23  nz = 0
24  az = cabs(z)
25  IF (az.EQ.0.0e0) go to 150
26  x = REAL(z)
27  arm = 1.0e+3*r1mach(1)
28  rtr1 = sqrt(arm)
29  crsc = cmplx(1.0e0,0.0e0)
30  iflag = 0
31  IF (az.LT.arm) go to 140
32  hz = z*cmplx(0.5e0,0.0e0)
33  cz = czero
34  IF (az.GT.rtr1) cz = hz*hz
35  acz = cabs(cz)
36  nn = n
37  ck = clog(hz)
38  10 CONTINUE
39  dfnu = fnu + float(nn-1)
40  fnup = dfnu + 1.0e0
41 C-----------------------------------------------------------------------
42 C UNDERFLOW TEST
43 C-----------------------------------------------------------------------
44  ak1 = ck*cmplx(dfnu,0.0e0)
45  ak = gamln(fnup,idum)
46  ak1 = ak1 - cmplx(ak,0.0e0)
47  IF (kode.EQ.2) ak1 = ak1 - cmplx(x,0.0e0)
48  rak1 = REAL(ak1)
49  IF (rak1.GT.(-elim)) go to 30
50  20 CONTINUE
51  nz = nz + 1
52  y(nn) = czero
53  IF (acz.GT.dfnu) go to 170
54  nn = nn - 1
55  IF (nn.EQ.0) RETURN
56  go to 10
57  30 CONTINUE
58  IF (rak1.GT.(-alim)) go to 40
59  iflag = 1
60  ss = 1.0e0/tol
61  crsc = cmplx(tol,0.0e0)
62  ascle = arm*ss
63  40 CONTINUE
64  ak = aimag(ak1)
65  aa = exp(rak1)
66  IF (iflag.EQ.1) aa = aa*ss
67  coef = cmplx(aa,0.0e0)*cmplx(cos(ak),sin(ak))
68  atol = tol*acz/fnup
69  il = min0(2,nn)
70  DO 80 i=1,il
71  dfnu = fnu + float(nn-i)
72  fnup = dfnu + 1.0e0
73  s1 = cone
74  IF (acz.LT.tol*fnup) go to 60
75  ak1 = cone
76  ak = fnup + 2.0e0
77  s = fnup
78  aa = 2.0e0
79  50 CONTINUE
80  rs = 1.0e0/s
81  ak1 = ak1*cz*cmplx(rs,0.0e0)
82  s1 = s1 + ak1
83  s = s + ak
84  ak = ak + 2.0e0
85  aa = aa*acz*rs
86  IF (aa.GT.atol) go to 50
87  60 CONTINUE
88  m = nn - i + 1
89  s2 = s1*coef
90  w(i) = s2
91  IF (iflag.EQ.0) go to 70
92  CALL cuchk(s2, nw, ascle, tol)
93  IF (nw.NE.0) go to 20
94  70 CONTINUE
95  y(m) = s2*crsc
96  IF (i.NE.il) coef = coef*cmplx(dfnu,0.0e0)/hz
97  80 CONTINUE
98  IF (nn.LE.2) RETURN
99  k = nn - 2
100  ak = float(k)
101  rz = (cone+cone)/z
102  IF (iflag.EQ.1) go to 110
103  ib = 3
104  90 CONTINUE
105  DO 100 i=ib,nn
106  y(k) = cmplx(ak+fnu,0.0e0)*rz*y(k+1) + y(k+2)
107  ak = ak - 1.0e0
108  k = k - 1
109  100 CONTINUE
110  RETURN
111 C-----------------------------------------------------------------------
112 C RECUR BACKWARD WITH SCALED VALUES
113 C-----------------------------------------------------------------------
114  110 CONTINUE
115 C-----------------------------------------------------------------------
116 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
117 C UNDERFLOW LIMIT = ASCLE = R1MACH(1)*CSCL*1.0E+3
118 C-----------------------------------------------------------------------
119  s1 = w(1)
120  s2 = w(2)
121  DO 120 l=3,nn
122  ck = s2
123  s2 = s1 + cmplx(ak+fnu,0.0e0)*rz*s2
124  s1 = ck
125  ck = s2*crsc
126  y(k) = ck
127  ak = ak - 1.0e0
128  k = k - 1
129  IF (cabs(ck).GT.ascle) go to 130
130  120 CONTINUE
131  RETURN
132  130 CONTINUE
133  ib = l + 1
134  IF (ib.GT.nn) RETURN
135  go to 90
136  140 CONTINUE
137  nz = n
138  IF (fnu.EQ.0.0e0) nz = nz - 1
139  150 CONTINUE
140  y(1) = czero
141  IF (fnu.EQ.0.0e0) y(1) = cone
142  IF (n.EQ.1) RETURN
143  DO 160 i=2,n
144  y(i) = czero
145  160 CONTINUE
146  RETURN
147 C-----------------------------------------------------------------------
148 C RETURN WITH NZ.LT.0 IF CABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE
149 C THE CALCULATION IN CBINU WITH N=N-IABS(NZ)
150 C-----------------------------------------------------------------------
151  170 CONTINUE
152  nz = -nz
153  RETURN
154  END
subroutine cseri(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
Definition: cseri.f:1
std::string dimension(void) const
int exp(void) const
Definition: DET.h:64
subroutine cuchk(Y, NZ, ASCLE, TOL)
Definition: cuchk.f:1
double complex cmplx
Definition: Faddeeva.cc:226
octave_value sin(void) const
Definition: ov.h:1198
octave_value cos(void) const
Definition: ov.h:1170
octave_value sqrt(void) const
Definition: ov.h:1200
F77_RET_T const double * x
function gamln(Z, IERR)
Definition: gamln.f:1