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
zuni1.f
Go to the documentation of this file.
1  SUBROUTINE zuni1(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
2  * tol, elim, alim)
3 C***BEGIN PROLOGUE ZUNI1
4 C***REFER TO ZBESI,ZBESK
5 C
6 C ZUNI1 COMPUTES I(FNU,Z) BY MEANS OF THE UNIFORM ASYMPTOTIC
7 C EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.
8 C
9 C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
10 C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
11 C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
12 C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
13 C Y(I)=CZERO FOR I=NLAST+1,N
14 C
15 C***ROUTINES CALLED ZUCHK,ZUNIK,ZUOIK,D1MACH,XZABS
16 C***END PROLOGUE ZUNI1
17 C COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1,
18 C *S2,Y,Z,ZETA1,ZETA2
19  DOUBLE PRECISION alim, aphi, ascle, bry, coner, crsc,
20  * cscl, csrr, cssr, cwrki, cwrkr, c1r, c2i, c2m, c2r, elim, fn,
21  * fnu, fnul, phii, phir, rast, rs1, rzi, rzr, sti, str, sumi,
22  * sumr, s1i, s1r, s2i, s2r, tol, yi, yr, zeroi, zeror, zeta1i,
23  * zeta1r, zeta2i, zeta2r, zi, zr, cyr, cyi, d1mach, xzabs
24  INTEGER i, iflag, init, k, kode, m, n, nd, nlast, nn, nuf, nw, nz
25  dimension bry(3), yr(n), yi(n), cwrkr(16), cwrki(16), cssr(3),
26  * csrr(3), cyr(2), cyi(2)
27  DATA zeror,zeroi,coner / 0.0d0, 0.0d0, 1.0d0 /
28 C
29  nz = 0
30  nd = n
31  nlast = 0
32 C-----------------------------------------------------------------------
33 C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
34 C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
35 C EXP(ALIM)=EXP(ELIM)*TOL
36 C-----------------------------------------------------------------------
37  cscl = 1.0d0/tol
38  crsc = tol
39  cssr(1) = cscl
40  cssr(2) = coner
41  cssr(3) = crsc
42  csrr(1) = crsc
43  csrr(2) = coner
44  csrr(3) = cscl
45  bry(1) = 1.0d+3*d1mach(1)/tol
46 C-----------------------------------------------------------------------
47 C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
48 C-----------------------------------------------------------------------
49  fn = dmax1(fnu,1.0d0)
50  init = 0
51  CALL zunik(zr, zi, fn, 1, 1, tol, init, phir, phii, zeta1r,
52  * zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
53  IF (kode.EQ.1) go to 10
54  str = zr + zeta2r
55  sti = zi + zeta2i
56  rast = fn/xzabs(str,sti)
57  str = str*rast*rast
58  sti = -sti*rast*rast
59  s1r = -zeta1r + str
60  s1i = -zeta1i + sti
61  go to 20
62  10 CONTINUE
63  s1r = -zeta1r + zeta2r
64  s1i = -zeta1i + zeta2i
65  20 CONTINUE
66  rs1 = s1r
67  IF (dabs(rs1).GT.elim) go to 130
68  30 CONTINUE
69  nn = min0(2,nd)
70  DO 80 i=1,nn
71  fn = fnu + dble(float(nd-i))
72  init = 0
73  CALL zunik(zr, zi, fn, 1, 0, tol, init, phir, phii, zeta1r,
74  * zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
75  IF (kode.EQ.1) go to 40
76  str = zr + zeta2r
77  sti = zi + zeta2i
78  rast = fn/xzabs(str,sti)
79  str = str*rast*rast
80  sti = -sti*rast*rast
81  s1r = -zeta1r + str
82  s1i = -zeta1i + sti + zi
83  go to 50
84  40 CONTINUE
85  s1r = -zeta1r + zeta2r
86  s1i = -zeta1i + zeta2i
87  50 CONTINUE
88 C-----------------------------------------------------------------------
89 C TEST FOR UNDERFLOW AND OVERFLOW
90 C-----------------------------------------------------------------------
91  rs1 = s1r
92  IF (dabs(rs1).GT.elim) go to 110
93  IF (i.EQ.1) iflag = 2
94  IF (dabs(rs1).LT.alim) go to 60
95 C-----------------------------------------------------------------------
96 C REFINE TEST AND SCALE
97 C-----------------------------------------------------------------------
98  aphi = xzabs(phir,phii)
99  rs1 = rs1 + dlog(aphi)
100  IF (dabs(rs1).GT.elim) go to 110
101  IF (i.EQ.1) iflag = 1
102  IF (rs1.LT.0.0d0) go to 60
103  IF (i.EQ.1) iflag = 3
104  60 CONTINUE
105 C-----------------------------------------------------------------------
106 C SCALE S1 IF CABS(S1).LT.ASCLE
107 C-----------------------------------------------------------------------
108  s2r = phir*sumr - phii*sumi
109  s2i = phir*sumi + phii*sumr
110  str = dexp(s1r)*cssr(iflag)
111  s1r = str*dcos(s1i)
112  s1i = str*dsin(s1i)
113  str = s2r*s1r - s2i*s1i
114  s2i = s2r*s1i + s2i*s1r
115  s2r = str
116  IF (iflag.NE.1) go to 70
117  CALL zuchk(s2r, s2i, nw, bry(1), tol)
118  IF (nw.NE.0) go to 110
119  70 CONTINUE
120  cyr(i) = s2r
121  cyi(i) = s2i
122  m = nd - i + 1
123  yr(m) = s2r*csrr(iflag)
124  yi(m) = s2i*csrr(iflag)
125  80 CONTINUE
126  IF (nd.LE.2) go to 100
127  rast = 1.0d0/xzabs(zr,zi)
128  str = zr*rast
129  sti = -zi*rast
130  rzr = (str+str)*rast
131  rzi = (sti+sti)*rast
132  bry(2) = 1.0d0/bry(1)
133  bry(3) = d1mach(2)
134  s1r = cyr(1)
135  s1i = cyi(1)
136  s2r = cyr(2)
137  s2i = cyi(2)
138  c1r = csrr(iflag)
139  ascle = bry(iflag)
140  k = nd - 2
141  fn = dble(float(k))
142  DO 90 i=3,nd
143  c2r = s2r
144  c2i = s2i
145  s2r = s1r + (fnu+fn)*(rzr*c2r-rzi*c2i)
146  s2i = s1i + (fnu+fn)*(rzr*c2i+rzi*c2r)
147  s1r = c2r
148  s1i = c2i
149  c2r = s2r*c1r
150  c2i = s2i*c1r
151  yr(k) = c2r
152  yi(k) = c2i
153  k = k - 1
154  fn = fn - 1.0d0
155  IF (iflag.GE.3) go to 90
156  str = dabs(c2r)
157  sti = dabs(c2i)
158  c2m = dmax1(str,sti)
159  IF (c2m.LE.ascle) go to 90
160  iflag = iflag + 1
161  ascle = bry(iflag)
162  s1r = s1r*c1r
163  s1i = s1i*c1r
164  s2r = c2r
165  s2i = c2i
166  s1r = s1r*cssr(iflag)
167  s1i = s1i*cssr(iflag)
168  s2r = s2r*cssr(iflag)
169  s2i = s2i*cssr(iflag)
170  c1r = csrr(iflag)
171  90 CONTINUE
172  100 CONTINUE
173  RETURN
174 C-----------------------------------------------------------------------
175 C SET UNDERFLOW AND UPDATE PARAMETERS
176 C-----------------------------------------------------------------------
177  110 CONTINUE
178  IF (rs1.GT.0.0d0) go to 120
179  yr(nd) = zeror
180  yi(nd) = zeroi
181  nz = nz + 1
182  nd = nd - 1
183  IF (nd.EQ.0) go to 100
184  CALL zuoik(zr, zi, fnu, kode, 1, nd, yr, yi, nuf, tol, elim, alim)
185  IF (nuf.LT.0) go to 120
186  nd = nd - nuf
187  nz = nz + nuf
188  IF (nd.EQ.0) go to 100
189  fn = fnu + dble(float(nd-1))
190  IF (fn.GE.fnul) go to 30
191  nlast = nd
192  RETURN
193  120 CONTINUE
194  nz = -1
195  RETURN
196  130 CONTINUE
197  IF (rs1.GT.0.0d0) go to 120
198  nz = n
199  DO 140 i=1,n
200  yr(i) = zeror
201  yi(i) = zeroi
202  140 CONTINUE
203  RETURN
204  END