GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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)
3 C***BEGIN PROLOGUE ZUNI2
4 C***REFER TO ZBESI,ZBESK
5 C
6 C ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
7 C UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
8 C OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
9 C
10 C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
11 C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
12 C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
13 C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
14 C Y(I)=CZERO FOR I=NLAST+1,N
15 C
16 C***ROUTINES CALLED ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,XZABS
17 C***END PROLOGUE ZUNI2
18 C COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
19 C *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/
36 C
37  nz = 0
38  nd = n
39  nlast = 0
40 C-----------------------------------------------------------------------
41 C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
42 C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
43 C EXP(ALIM)=EXP(ELIM)*TOL
44 C-----------------------------------------------------------------------
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
54 C-----------------------------------------------------------------------
55 C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
56 C-----------------------------------------------------------------------
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
79 C-----------------------------------------------------------------------
80 C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
81 C-----------------------------------------------------------------------
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
119 C-----------------------------------------------------------------------
120 C TEST FOR UNDERFLOW AND OVERFLOW
121 C-----------------------------------------------------------------------
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
126 C-----------------------------------------------------------------------
127 C REFINE TEST AND SCALE
128 C-----------------------------------------------------------------------
129 C-----------------------------------------------------------------------
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
138 C-----------------------------------------------------------------------
139 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
140 C EXPONENT EXTREMES
141 C-----------------------------------------------------------------------
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
223 C-----------------------------------------------------------------------
224 C SET UNDERFLOW AND UPDATE PARAMETERS
225 C-----------------------------------------------------------------------
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
238 C FN = CIDI
239 C J = NUF + 1
240 C K = MOD(J,4) + 1
241 C S1R = CIPR(K)
242 C S1I = CIPI(K)
243 C IF (FN.LT.0.0D0) S1I = -S1I
244 C STR = C2R*S1R - C2I*S1I
245 C C2I = C2R*S1I + C2I*S1R
246 C 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
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:932
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