GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
zuni2.f
Go to the documentation of this file.
1 SUBROUTINE zuni2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
2 * TOL, ELIM, ALIM)
3C***BEGIN PROLOGUE ZUNI2
4C***REFER TO ZBESI,ZBESK
5C
6C ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
7C UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
8C OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
9C
10C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
11C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
12C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
13C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
14C Y(I)=CZERO FOR I=NLAST+1,N
15C
16C***ROUTINES CALLED ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,XZABS
17C***END PROLOGUE ZUNI2
18C COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
19C *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN
20 DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI,
21 * argr, ascle, asumi, asumr, bry, bsumi, bsumr, cidi, cipi, cipr,
22 * coner, crsc, cscl, csrr, cssr, c1r, c2i, c2m, c2r, daii,
23 * dair, elim, fn, fnu, fnul, hpi, phii, phir, rast, raz, rs1, rzi,
24 * rzr, sti, str, s1i, s1r, s2i, s2r, tol, yi, yr, zbi, zbr, zeroi,
25 * zeror, zeta1i, zeta1r, zeta2i, zeta2r, zi, zni, znr, zr, cyr,
26 * cyi, d1mach, xzabs, car, sar
27 INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
28 * nn, nuf, nw, nz, idum
29 dimension bry(3), yr(n), yi(n), cipr(4), cipi(4), cssr(3),
30 * csrr(3), cyr(2), cyi(2)
31 DATA zeror,zeroi,coner / 0.0d0, 0.0d0, 1.0d0 /
32 DATA cipr(1),cipi(1),cipr(2),cipi(2),cipr(3),cipi(3),cipr(4),
33 * cipi(4)/ 1.0d0,0.0d0, 0.0d0,1.0d0, -1.0d0,0.0d0, 0.0d0,-1.0d0/
34 DATA hpi, aic /
35 1 1.57079632679489662d+00, 1.265512123484645396d+00/
36C
37 nz = 0
38 nd = n
39 nlast = 0
40C-----------------------------------------------------------------------
41C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
42C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
43C EXP(ALIM)=EXP(ELIM)*TOL
44C-----------------------------------------------------------------------
45 cscl = 1.0d0/tol
46 crsc = tol
47 cssr(1) = cscl
48 cssr(2) = coner
49 cssr(3) = crsc
50 csrr(1) = crsc
51 csrr(2) = coner
52 csrr(3) = cscl
53 bry(1) = 1.0d+3*d1mach(1)/tol
54C-----------------------------------------------------------------------
55C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
56C-----------------------------------------------------------------------
57 znr = zi
58 zni = -zr
59 zbr = zr
60 zbi = zi
61 cidi = -coner
62 inu = int(sngl(fnu))
63 ang = hpi*(fnu-dble(float(inu)))
64 c2r = dcos(ang)
65 c2i = dsin(ang)
66 car = c2r
67 sar = c2i
68 in = inu + n - 1
69 in = mod(in,4) + 1
70 str = c2r*cipr(in) - c2i*cipi(in)
71 c2i = c2r*cipi(in) + c2i*cipr(in)
72 c2r = str
73 IF (zi.GT.0.0d0) GO TO 10
74 znr = -znr
75 zbi = -zbi
76 cidi = -cidi
77 c2i = -c2i
78 10 CONTINUE
79C-----------------------------------------------------------------------
80C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
81C-----------------------------------------------------------------------
82 fn = dmax1(fnu,1.0d0)
83 CALL zunhj(znr, zni, fn, 1, tol, phir, phii, argr, argi, zeta1r,
84 * zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
85 IF (kode.EQ.1) GO TO 20
86 str = zbr + zeta2r
87 sti = zbi + zeta2i
88 rast = fn/xzabs(str,sti)
89 str = str*rast*rast
90 sti = -sti*rast*rast
91 s1r = -zeta1r + str
92 s1i = -zeta1i + sti
93 GO TO 30
94 20 CONTINUE
95 s1r = -zeta1r + zeta2r
96 s1i = -zeta1i + zeta2i
97 30 CONTINUE
98 rs1 = s1r
99 IF (dabs(rs1).GT.elim) GO TO 150
100 40 CONTINUE
101 nn = min0(2,nd)
102 DO 90 i=1,nn
103 fn = fnu + dble(float(nd-i))
104 CALL zunhj(znr, zni, fn, 0, tol, phir, phii, argr, argi,
105 * zeta1r, zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
106 IF (kode.EQ.1) GO TO 50
107 str = zbr + zeta2r
108 sti = zbi + zeta2i
109 rast = fn/xzabs(str,sti)
110 str = str*rast*rast
111 sti = -sti*rast*rast
112 s1r = -zeta1r + str
113 s1i = -zeta1i + sti + dabs(zi)
114 GO TO 60
115 50 CONTINUE
116 s1r = -zeta1r + zeta2r
117 s1i = -zeta1i + zeta2i
118 60 CONTINUE
119C-----------------------------------------------------------------------
120C TEST FOR UNDERFLOW AND OVERFLOW
121C-----------------------------------------------------------------------
122 rs1 = s1r
123 IF (dabs(rs1).GT.elim) GO TO 120
124 IF (i.EQ.1) iflag = 2
125 IF (dabs(rs1).LT.alim) GO TO 70
126C-----------------------------------------------------------------------
127C REFINE TEST AND SCALE
128C-----------------------------------------------------------------------
129C-----------------------------------------------------------------------
130 aphi = xzabs(phir,phii)
131 aarg = xzabs(argr,argi)
132 rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
133 IF (dabs(rs1).GT.elim) GO TO 120
134 IF (i.EQ.1) iflag = 1
135 IF (rs1.LT.0.0d0) GO TO 70
136 IF (i.EQ.1) iflag = 3
137 70 CONTINUE
138C-----------------------------------------------------------------------
139C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
140C EXPONENT EXTREMES
141C-----------------------------------------------------------------------
142 CALL zairy(argr, argi, 0, 2, air, aii, nai, idum)
143 CALL zairy(argr, argi, 1, 2, dair, daii, ndai, idum)
144 str = dair*bsumr - daii*bsumi
145 sti = dair*bsumi + daii*bsumr
146 str = str + (air*asumr-aii*asumi)
147 sti = sti + (air*asumi+aii*asumr)
148 s2r = phir*str - phii*sti
149 s2i = phir*sti + phii*str
150 str = dexp(s1r)*cssr(iflag)
151 s1r = str*dcos(s1i)
152 s1i = str*dsin(s1i)
153 str = s2r*s1r - s2i*s1i
154 s2i = s2r*s1i + s2i*s1r
155 s2r = str
156 IF (iflag.NE.1) GO TO 80
157 CALL zuchk(s2r, s2i, nw, bry(1), tol)
158 IF (nw.NE.0) GO TO 120
159 80 CONTINUE
160 IF (zi.LE.0.0d0) s2i = -s2i
161 str = s2r*c2r - s2i*c2i
162 s2i = s2r*c2i + s2i*c2r
163 s2r = str
164 cyr(i) = s2r
165 cyi(i) = s2i
166 j = nd - i + 1
167 yr(j) = s2r*csrr(iflag)
168 yi(j) = s2i*csrr(iflag)
169 str = -c2i*cidi
170 c2i = c2r*cidi
171 c2r = str
172 90 CONTINUE
173 IF (nd.LE.2) GO TO 110
174 raz = 1.0d0/xzabs(zr,zi)
175 str = zr*raz
176 sti = -zi*raz
177 rzr = (str+str)*raz
178 rzi = (sti+sti)*raz
179 bry(2) = 1.0d0/bry(1)
180 bry(3) = d1mach(2)
181 s1r = cyr(1)
182 s1i = cyi(1)
183 s2r = cyr(2)
184 s2i = cyi(2)
185 c1r = csrr(iflag)
186 ascle = bry(iflag)
187 k = nd - 2
188 fn = dble(float(k))
189 DO 100 i=3,nd
190 c2r = s2r
191 c2i = s2i
192 s2r = s1r + (fnu+fn)*(rzr*c2r-rzi*c2i)
193 s2i = s1i + (fnu+fn)*(rzr*c2i+rzi*c2r)
194 s1r = c2r
195 s1i = c2i
196 c2r = s2r*c1r
197 c2i = s2i*c1r
198 yr(k) = c2r
199 yi(k) = c2i
200 k = k - 1
201 fn = fn - 1.0d0
202 IF (iflag.GE.3) GO TO 100
203 str = dabs(c2r)
204 sti = dabs(c2i)
205 c2m = dmax1(str,sti)
206 IF (c2m.LE.ascle) GO TO 100
207 iflag = iflag + 1
208 ascle = bry(iflag)
209 s1r = s1r*c1r
210 s1i = s1i*c1r
211 s2r = c2r
212 s2i = c2i
213 s1r = s1r*cssr(iflag)
214 s1i = s1i*cssr(iflag)
215 s2r = s2r*cssr(iflag)
216 s2i = s2i*cssr(iflag)
217 c1r = csrr(iflag)
218 100 CONTINUE
219 110 CONTINUE
220 RETURN
221 120 CONTINUE
222 IF (rs1.GT.0.0d0) GO TO 140
223C-----------------------------------------------------------------------
224C SET UNDERFLOW AND UPDATE PARAMETERS
225C-----------------------------------------------------------------------
226 yr(nd) = zeror
227 yi(nd) = zeroi
228 nz = nz + 1
229 nd = nd - 1
230 IF (nd.EQ.0) GO TO 110
231 CALL zuoik(zr, zi, fnu, kode, 1, nd, yr, yi, nuf, tol, elim, alim)
232 IF (nuf.LT.0) GO TO 140
233 nd = nd - nuf
234 nz = nz + nuf
235 IF (nd.EQ.0) GO TO 110
236 fn = fnu + dble(float(nd-1))
237 IF (fn.LT.fnul) GO TO 130
238C FN = CIDI
239C J = NUF + 1
240C K = MOD(J,4) + 1
241C S1R = CIPR(K)
242C S1I = CIPI(K)
243C IF (FN.LT.0.0D0) S1I = -S1I
244C STR = C2R*S1R - C2I*S1I
245C C2I = C2R*S1I + C2I*S1R
246C C2R = STR
247 in = inu + nd - 1
248 in = mod(in,4) + 1
249 c2r = car*cipr(in) - sar*cipi(in)
250 c2i = car*cipi(in) + sar*cipr(in)
251 IF (zi.LE.0.0d0) c2i = -c2i
252 GO TO 40
253 130 CONTINUE
254 nlast = nd
255 RETURN
256 140 CONTINUE
257 nz = -1
258 RETURN
259 150 CONTINUE
260 IF (rs1.GT.0.0d0) GO TO 140
261 nz = n
262 DO 160 i=1,n
263 yr(i) = zeror
264 yi(i) = zeroi
265 160 CONTINUE
266 RETURN
267 END
d1mach
double precision function d1mach(i)
Definition d1mach.f:23
mod
T mod(T x, T y)
Definition lo-mappers.h:294
xzabs
double precision function xzabs(zr, zi)
Definition xzabs.f:2
zairy
subroutine zairy(zr, zi, id, kode, air, aii, nz, ierr)
Definition zairy.f:2
zuchk
subroutine zuchk(yr, yi, nz, ascle, tol)
Definition zuchk.f:2
zunhj
subroutine zunhj(zr, zi, fnu, ipmtr, tol, phir, phii, argr, argi, zeta1r, zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
Definition zunhj.f:3
zuni2
subroutine zuni2(zr, zi, fnu, kode, n, yr, yi, nz, nlast, fnul, tol, elim, alim)
Definition zuni2.f:3
zuoik
subroutine zuoik(zr, zi, fnu, kode, ikflg, n, yr, yi, nuf, tol, elim, alim)
Definition zuoik.f:3