GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zuoik.f
Go to the documentation of this file.
1  SUBROUTINE zuoik(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL,
2  * ELIM, ALIM)
3 C***BEGIN PROLOGUE ZUOIK
4 C***REFER TO ZBESI,ZBESK,ZBESH
5 C
6 C ZUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
7 C EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
8 C (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
9 C WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
10 C EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
11 C THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
12 C MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
13 C EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
14 C EXP(-ELIM)/TOL
15 C
16 C IKFLG=1 MEANS THE I SEQUENCE IS TESTED
17 C =2 MEANS THE K SEQUENCE IS TESTED
18 C NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
19 C =-1 MEANS AN OVERFLOW WOULD OCCUR
20 C IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
21 C THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
22 C IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
23 C IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
24 C ANOTHER ROUTINE
25 C
26 C***ROUTINES CALLED ZUCHK,ZUNHJ,ZUNIK,D1MACH,XZABS,XZLOG
27 C***END PROLOGUE ZUOIK
28 C COMPLEX ARG,ASUM,BSUM,CWRK,CZ,CZERO,PHI,SUM,Y,Z,ZB,ZETA1,ZETA2,ZN,
29 C *ZR
30  DOUBLE PRECISION AARG, AIC, ALIM, APHI, ARGI, ARGR, ASUMI, ASUMR,
31  * ascle, ax, ay, bsumi, bsumr, cwrki, cwrkr, czi, czr, elim, fnn,
32  * fnu, gnn, gnu, phii, phir, rcz, str, sti, sumi, sumr, tol, yi,
33  * yr, zbi, zbr, zeroi, zeror, zeta1i, zeta1r, zeta2i, zeta2r, zi,
34  * zni, znr, zr, zri, zrr, d1mach, xzabs
35  INTEGER I, IDUM, IFORM, IKFLG, INIT, KODE, N, NN, NUF, NW
36  dimension yr(n), yi(n), cwrkr(16), cwrki(16)
37  DATA zeror,zeroi / 0.0d0, 0.0d0 /
38  DATA aic / 1.265512123484645396d+00 /
39  nuf = 0
40  nn = n
41  zrr = zr
42  zri = zi
43  IF (zr.GE.0.0d0) GO TO 10
44  zrr = -zr
45  zri = -zi
46  10 CONTINUE
47  zbr = zrr
48  zbi = zri
49  ax = dabs(zr)*1.7321d0
50  ay = dabs(zi)
51  iform = 1
52  IF (ay.GT.ax) iform = 2
53  gnu = dmax1(fnu,1.0d0)
54  IF (ikflg.EQ.1) GO TO 20
55  fnn = dble(float(nn))
56  gnn = fnu + fnn - 1.0d0
57  gnu = dmax1(gnn,fnn)
58  20 CONTINUE
59 C-----------------------------------------------------------------------
60 C ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
61 C REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
62 C THE SIGN OF THE IMAGINARY PART CORRECT.
63 C-----------------------------------------------------------------------
64  IF (iform.EQ.2) GO TO 30
65  init = 0
66  CALL zunik(zrr, zri, gnu, ikflg, 1, tol, init, phir, phii,
67  * zeta1r, zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
68  czr = -zeta1r + zeta2r
69  czi = -zeta1i + zeta2i
70  GO TO 50
71  30 CONTINUE
72  znr = zri
73  zni = -zrr
74  IF (zi.GT.0.0d0) GO TO 40
75  znr = -znr
76  40 CONTINUE
77  CALL zunhj(znr, zni, gnu, 1, tol, phir, phii, argr, argi, zeta1r,
78  * zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
79  czr = -zeta1r + zeta2r
80  czi = -zeta1i + zeta2i
81  aarg = xzabs(argr,argi)
82  50 CONTINUE
83  IF (kode.EQ.1) GO TO 60
84  czr = czr - zbr
85  czi = czi - zbi
86  60 CONTINUE
87  IF (ikflg.EQ.1) GO TO 70
88  czr = -czr
89  czi = -czi
90  70 CONTINUE
91  aphi = xzabs(phir,phii)
92  rcz = czr
93 C-----------------------------------------------------------------------
94 C OVERFLOW TEST
95 C-----------------------------------------------------------------------
96  IF (rcz.GT.elim) GO TO 210
97  IF (rcz.LT.alim) GO TO 80
98  rcz = rcz + dlog(aphi)
99  IF (iform.EQ.2) rcz = rcz - 0.25d0*dlog(aarg) - aic
100  IF (rcz.GT.elim) GO TO 210
101  GO TO 130
102  80 CONTINUE
103 C-----------------------------------------------------------------------
104 C UNDERFLOW TEST
105 C-----------------------------------------------------------------------
106  IF (rcz.LT.(-elim)) GO TO 90
107  IF (rcz.GT.(-alim)) GO TO 130
108  rcz = rcz + dlog(aphi)
109  IF (iform.EQ.2) rcz = rcz - 0.25d0*dlog(aarg) - aic
110  IF (rcz.GT.(-elim)) GO TO 110
111  90 CONTINUE
112  DO 100 i=1,nn
113  yr(i) = zeror
114  yi(i) = zeroi
115  100 CONTINUE
116  nuf = nn
117  RETURN
118  110 CONTINUE
119  ascle = 1.0d+3*d1mach(1)/tol
120  CALL xzlog(phir, phii, str, sti, idum)
121  czr = czr + str
122  czi = czi + sti
123  IF (iform.EQ.1) GO TO 120
124  CALL xzlog(argr, argi, str, sti, idum)
125  czr = czr - 0.25d0*str - aic
126  czi = czi - 0.25d0*sti
127  120 CONTINUE
128  ax = dexp(rcz)/tol
129  ay = czi
130  czr = ax*dcos(ay)
131  czi = ax*dsin(ay)
132  CALL zuchk(czr, czi, nw, ascle, tol)
133  IF (nw.NE.0) GO TO 90
134  130 CONTINUE
135  IF (ikflg.EQ.2) RETURN
136  IF (n.EQ.1) RETURN
137 C-----------------------------------------------------------------------
138 C SET UNDERFLOWS ON I SEQUENCE
139 C-----------------------------------------------------------------------
140  140 CONTINUE
141  gnu = fnu + dble(float(nn-1))
142  IF (iform.EQ.2) GO TO 150
143  init = 0
144  CALL zunik(zrr, zri, gnu, ikflg, 1, tol, init, phir, phii,
145  * zeta1r, zeta1i, zeta2r, zeta2i, sumr, sumi, cwrkr, cwrki)
146  czr = -zeta1r + zeta2r
147  czi = -zeta1i + zeta2i
148  GO TO 160
149  150 CONTINUE
150  CALL zunhj(znr, zni, gnu, 1, tol, phir, phii, argr, argi, zeta1r,
151  * zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
152  czr = -zeta1r + zeta2r
153  czi = -zeta1i + zeta2i
154  aarg = xzabs(argr,argi)
155  160 CONTINUE
156  IF (kode.EQ.1) GO TO 170
157  czr = czr - zbr
158  czi = czi - zbi
159  170 CONTINUE
160  aphi = xzabs(phir,phii)
161  rcz = czr
162  IF (rcz.LT.(-elim)) GO TO 180
163  IF (rcz.GT.(-alim)) RETURN
164  rcz = rcz + dlog(aphi)
165  IF (iform.EQ.2) rcz = rcz - 0.25d0*dlog(aarg) - aic
166  IF (rcz.GT.(-elim)) GO TO 190
167  180 CONTINUE
168  yr(nn) = zeror
169  yi(nn) = zeroi
170  nn = nn - 1
171  nuf = nuf + 1
172  IF (nn.EQ.0) RETURN
173  GO TO 140
174  190 CONTINUE
175  ascle = 1.0d+3*d1mach(1)/tol
176  CALL xzlog(phir, phii, str, sti, idum)
177  czr = czr + str
178  czi = czi + sti
179  IF (iform.EQ.1) GO TO 200
180  CALL xzlog(argr, argi, str, sti, idum)
181  czr = czr - 0.25d0*str - aic
182  czi = czi - 0.25d0*sti
183  200 CONTINUE
184  ax = dexp(rcz)/tol
185  ay = czi
186  czr = ax*dcos(ay)
187  czi = ax*dsin(ay)
188  CALL zuchk(czr, czi, nw, ascle, tol)
189  IF (nw.NE.0) GO TO 180
190  RETURN
191  210 CONTINUE
192  nuf = -1
193  RETURN
194  END
double precision function d1mach(i)
Definition: d1mach.f:23
double precision function xzabs(ZR, ZI)
Definition: xzabs.f:2
subroutine xzlog(AR, AI, BR, BI, IERR)
Definition: xzlog.f:2
subroutine zuchk(YR, YI, NZ, ASCLE, TOL)
Definition: zuchk.f:2
subroutine zunhj(ZR, ZI, FNU, IPMTR, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
Definition: zunhj.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