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
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