GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
double precision function d1mach(i)
Definition: d1mach.f:23
double precision function xzabs(ZR, ZI)
Definition: xzabs.f:2
subroutine zuchk(YR, YI, NZ, ASCLE, TOL)
Definition: zuchk.f:2
subroutine zuni1(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
Definition: zuni1.f:3
subroutine zunik(ZRR, ZRI, FNU, IKFLG, IPMTR, TOL, INIT, PHIR, PHII, ZETA1R, ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
Definition: zunik.f:3
subroutine zuoik(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL, ELIM, ALIM)
Definition: zuoik.f:3