GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
double complex cmplx
Definition: Faddeeva.cc:217
subroutine cseri(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
Definition: cseri.f:2
subroutine cuchk(Y, NZ, ASCLE, TOL)
Definition: cuchk.f:2
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137