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
zseri.f
Go to the documentation of this file.
1  SUBROUTINE zseri(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM,
2  * alim)
3 C***BEGIN PROLOGUE ZSERI
4 C***REFER TO ZBESI,ZBESK
5 C
6 C ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
7 C MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE
8 C REGION CABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
9 C NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
10 C DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE
11 C CONDITION CABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE
12 C COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
13 C
14 C***ROUTINES CALLED DGAMLN,D1MACH,ZUCHK,XZABS,ZDIV,XZLOG,ZMLT
15 C***END PROLOGUE ZSERI
16 C COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z
17  DOUBLE PRECISION aa, acz, ak, ak1i, ak1r, alim, arm, ascle, atol,
18  * az, cki, ckr, coefi, coefr, conei, coner, crscr, czi, czr, dfnu,
19  * elim, fnu, fnup, hzi, hzr, raz, rs, rtr1, rzi, rzr, s, ss, sti,
20  * str, s1i, s1r, s2i, s2r, tol, yi, yr, wi, wr, zeroi, zeror, zi,
21  * zr, dgamln, d1mach, xzabs
22  INTEGER i, ib, idum, iflag, il, k, kode, l, m, n, nn, nz, nw
23  dimension yr(n), yi(n), wr(2), wi(2)
24  DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
25 C
26  nz = 0
27  az = xzabs(zr,zi)
28  IF (az.EQ.0.0d0) go to 160
29  arm = 1.0d+3*d1mach(1)
30  rtr1 = dsqrt(arm)
31  crscr = 1.0d0
32  iflag = 0
33  IF (az.LT.arm) go to 150
34  hzr = 0.5d0*zr
35  hzi = 0.5d0*zi
36  czr = zeror
37  czi = zeroi
38  IF (az.LE.rtr1) go to 10
39  CALL zmlt(hzr, hzi, hzr, hzi, czr, czi)
40  10 CONTINUE
41  acz = xzabs(czr,czi)
42  nn = n
43  CALL xzlog(hzr, hzi, ckr, cki, idum)
44  20 CONTINUE
45  dfnu = fnu + dble(float(nn-1))
46  fnup = dfnu + 1.0d0
47 C-----------------------------------------------------------------------
48 C UNDERFLOW TEST
49 C-----------------------------------------------------------------------
50  ak1r = ckr*dfnu
51  ak1i = cki*dfnu
52  ak = dgamln(fnup,idum)
53  ak1r = ak1r - ak
54  IF (kode.EQ.2) ak1r = ak1r - zr
55  IF (ak1r.GT.(-elim)) go to 40
56  30 CONTINUE
57  nz = nz + 1
58  yr(nn) = zeror
59  yi(nn) = zeroi
60  IF (acz.GT.dfnu) go to 190
61  nn = nn - 1
62  IF (nn.EQ.0) RETURN
63  go to 20
64  40 CONTINUE
65  IF (ak1r.GT.(-alim)) go to 50
66  iflag = 1
67  ss = 1.0d0/tol
68  crscr = tol
69  ascle = arm*ss
70  50 CONTINUE
71  aa = dexp(ak1r)
72  IF (iflag.EQ.1) aa = aa*ss
73  coefr = aa*dcos(ak1i)
74  coefi = aa*dsin(ak1i)
75  atol = tol*acz/fnup
76  il = min0(2,nn)
77  DO 90 i=1,il
78  dfnu = fnu + dble(float(nn-i))
79  fnup = dfnu + 1.0d0
80  s1r = coner
81  s1i = conei
82  IF (acz.LT.tol*fnup) go to 70
83  ak1r = coner
84  ak1i = conei
85  ak = fnup + 2.0d0
86  s = fnup
87  aa = 2.0d0
88  60 CONTINUE
89  rs = 1.0d0/s
90  str = ak1r*czr - ak1i*czi
91  sti = ak1r*czi + ak1i*czr
92  ak1r = str*rs
93  ak1i = sti*rs
94  s1r = s1r + ak1r
95  s1i = s1i + ak1i
96  s = s + ak
97  ak = ak + 2.0d0
98  aa = aa*acz*rs
99  IF (aa.GT.atol) go to 60
100  70 CONTINUE
101  s2r = s1r*coefr - s1i*coefi
102  s2i = s1r*coefi + s1i*coefr
103  wr(i) = s2r
104  wi(i) = s2i
105  IF (iflag.EQ.0) go to 80
106  CALL zuchk(s2r, s2i, nw, ascle, tol)
107  IF (nw.NE.0) go to 30
108  80 CONTINUE
109  m = nn - i + 1
110  yr(m) = s2r*crscr
111  yi(m) = s2i*crscr
112  IF (i.EQ.il) go to 90
113  CALL zdiv(coefr, coefi, hzr, hzi, str, sti)
114  coefr = str*dfnu
115  coefi = sti*dfnu
116  90 CONTINUE
117  IF (nn.LE.2) RETURN
118  k = nn - 2
119  ak = dble(float(k))
120  raz = 1.0d0/az
121  str = zr*raz
122  sti = -zi*raz
123  rzr = (str+str)*raz
124  rzi = (sti+sti)*raz
125  IF (iflag.EQ.1) go to 120
126  ib = 3
127  100 CONTINUE
128  DO 110 i=ib,nn
129  yr(k) = (ak+fnu)*(rzr*yr(k+1)-rzi*yi(k+1)) + yr(k+2)
130  yi(k) = (ak+fnu)*(rzr*yi(k+1)+rzi*yr(k+1)) + yi(k+2)
131  ak = ak - 1.0d0
132  k = k - 1
133  110 CONTINUE
134  RETURN
135 C-----------------------------------------------------------------------
136 C RECUR BACKWARD WITH SCALED VALUES
137 C-----------------------------------------------------------------------
138  120 CONTINUE
139 C-----------------------------------------------------------------------
140 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
141 C UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1.0D+3
142 C-----------------------------------------------------------------------
143  s1r = wr(1)
144  s1i = wi(1)
145  s2r = wr(2)
146  s2i = wi(2)
147  DO 130 l=3,nn
148  ckr = s2r
149  cki = s2i
150  s2r = s1r + (ak+fnu)*(rzr*ckr-rzi*cki)
151  s2i = s1i + (ak+fnu)*(rzr*cki+rzi*ckr)
152  s1r = ckr
153  s1i = cki
154  ckr = s2r*crscr
155  cki = s2i*crscr
156  yr(k) = ckr
157  yi(k) = cki
158  ak = ak - 1.0d0
159  k = k - 1
160  IF (xzabs(ckr,cki).GT.ascle) go to 140
161  130 CONTINUE
162  RETURN
163  140 CONTINUE
164  ib = l + 1
165  IF (ib.GT.nn) RETURN
166  go to 100
167  150 CONTINUE
168  nz = n
169  IF (fnu.EQ.0.0d0) nz = nz - 1
170  160 CONTINUE
171  yr(1) = zeror
172  yi(1) = zeroi
173  IF (fnu.NE.0.0d0) go to 170
174  yr(1) = coner
175  yi(1) = conei
176  170 CONTINUE
177  IF (n.EQ.1) RETURN
178  DO 180 i=2,n
179  yr(i) = zeror
180  yi(i) = zeroi
181  180 CONTINUE
182  RETURN
183 C-----------------------------------------------------------------------
184 C RETURN WITH NZ.LT.0 IF CABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE
185 C THE CALCULATION IN CBINU WITH N=N-IABS(NZ)
186 C-----------------------------------------------------------------------
187  190 CONTINUE
188  nz = -nz
189  RETURN
190  END