GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
double precision function d1mach(i)
Definition: d1mach.f:23
double precision function dgamln(Z, IERR)
Definition: dgamln.f:2
double precision function xzabs(ZR, ZI)
Definition: xzabs.f:2
subroutine xzlog(AR, AI, BR, BI, IERR)
Definition: xzlog.f:2
subroutine zdiv(AR, AI, BR, BI, CR, CI)
Definition: zdiv.f:2
subroutine zmlt(AR, AI, BR, BI, CR, CI)
Definition: zmlt.f:2
subroutine zseri(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM)
Definition: zseri.f:3
subroutine zuchk(YR, YI, NZ, ASCLE, TOL)
Definition: zuchk.f:2