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
cuoik.f
Go to the documentation of this file.
1  SUBROUTINE cuoik(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CUOIK
3 C***REFER TO CBESI,CBESK,CBESH
4 C
5 C CUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
6 C EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
7 C (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
8 C WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
9 C EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
10 C THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
11 C MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
12 C EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
13 C EXP(-ELIM)/TOL
14 C
15 C IKFLG=1 MEANS THE I SEQUENCE IS TESTED
16 C =2 MEANS THE K SEQUENCE IS TESTED
17 C NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
18 C =-1 MEANS AN OVERFLOW WOULD OCCUR
19 C IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
20 C THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
21 C IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
22 C IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
23 C ANOTHER ROUTINE
24 C
25 C***ROUTINES CALLED CUCHK,CUNHJ,CUNIK,R1MACH
26 C***END PROLOGUE CUOIK
27  COMPLEX arg, asum, bsum, cwrk, cz, czero, phi, sum, y, z, zb,
28  * zeta1, zeta2, zn, zr
29  REAL aarg, aic, alim, aphi, ascle, ax, ay, elim, fnn, fnu, gnn,
30  * gnu, rcz, tol, x, yy
31  INTEGER i, iform, ikflg, init, kode, n, nn, nuf, nw
32  dimension y(n), cwrk(16)
33  DATA czero / (0.0e0,0.0e0) /
34  DATA aic / 1.265512123484645396e+00 /
35  nuf = 0
36  nn = n
37  x = REAL(z)
38  zr = z
39  IF (x.LT.0.0e0) zr = -z
40  zb = zr
41  yy = aimag(zr)
42  ax = abs(x)*1.7321e0
43  ay = abs(yy)
44  iform = 1
45  IF (ay.GT.ax) iform = 2
46  gnu = amax1(fnu,1.0e0)
47  IF (ikflg.EQ.1) go to 10
48  fnn = float(nn)
49  gnn = fnu + fnn - 1.0e0
50  gnu = amax1(gnn,fnn)
51  10 CONTINUE
52 C-----------------------------------------------------------------------
53 C ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
54 C REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
55 C THE SIGN OF THE IMAGINARY PART CORRECT.
56 C-----------------------------------------------------------------------
57  IF (iform.EQ.2) go to 20
58  init = 0
59  CALL cunik(zr, gnu, ikflg, 1, tol, init, phi, zeta1, zeta2, sum,
60  * cwrk)
61  cz = -zeta1 + zeta2
62  go to 40
63  20 CONTINUE
64  zn = -zr*cmplx(0.0e0,1.0e0)
65  IF (yy.GT.0.0e0) go to 30
66  zn = conjg(-zn)
67  30 CONTINUE
68  CALL cunhj(zn, gnu, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
69  cz = -zeta1 + zeta2
70  aarg = cabs(arg)
71  40 CONTINUE
72  IF (kode.EQ.2) cz = cz - zb
73  IF (ikflg.EQ.2) cz = -cz
74  aphi = cabs(phi)
75  rcz = REAL(cz)
76 C-----------------------------------------------------------------------
77 C OVERFLOW TEST
78 C-----------------------------------------------------------------------
79  IF (rcz.GT.elim) go to 170
80  IF (rcz.LT.alim) go to 50
81  rcz = rcz + alog(aphi)
82  IF (iform.EQ.2) rcz = rcz - 0.25e0*alog(aarg) - aic
83  IF (rcz.GT.elim) go to 170
84  go to 100
85  50 CONTINUE
86 C-----------------------------------------------------------------------
87 C UNDERFLOW TEST
88 C-----------------------------------------------------------------------
89  IF (rcz.LT.(-elim)) go to 60
90  IF (rcz.GT.(-alim)) go to 100
91  rcz = rcz + alog(aphi)
92  IF (iform.EQ.2) rcz = rcz - 0.25e0*alog(aarg) - aic
93  IF (rcz.GT.(-elim)) go to 80
94  60 CONTINUE
95  DO 70 i=1,nn
96  y(i) = czero
97  70 CONTINUE
98  nuf = nn
99  RETURN
100  80 CONTINUE
101  ascle = 1.0e+3*r1mach(1)/tol
102  cz = cz + clog(phi)
103  IF (iform.EQ.1) go to 90
104  cz = cz - cmplx(0.25e0,0.0e0)*clog(arg) - cmplx(aic,0.0e0)
105  90 CONTINUE
106  ax = exp(rcz)/tol
107  ay = aimag(cz)
108  cz = cmplx(ax,0.0e0)*cmplx(cos(ay),sin(ay))
109  CALL cuchk(cz, nw, ascle, tol)
110  IF (nw.EQ.1) go to 60
111  100 CONTINUE
112  IF (ikflg.EQ.2) RETURN
113  IF (n.EQ.1) RETURN
114 C-----------------------------------------------------------------------
115 C SET UNDERFLOWS ON I SEQUENCE
116 C-----------------------------------------------------------------------
117  110 CONTINUE
118  gnu = fnu + float(nn-1)
119  IF (iform.EQ.2) go to 120
120  init = 0
121  CALL cunik(zr, gnu, ikflg, 1, tol, init, phi, zeta1, zeta2, sum,
122  * cwrk)
123  cz = -zeta1 + zeta2
124  go to 130
125  120 CONTINUE
126  CALL cunhj(zn, gnu, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
127  cz = -zeta1 + zeta2
128  aarg = cabs(arg)
129  130 CONTINUE
130  IF (kode.EQ.2) cz = cz - zb
131  aphi = cabs(phi)
132  rcz = REAL(cz)
133  IF (rcz.LT.(-elim)) go to 140
134  IF (rcz.GT.(-alim)) RETURN
135  rcz = rcz + alog(aphi)
136  IF (iform.EQ.2) rcz = rcz - 0.25e0*alog(aarg) - aic
137  IF (rcz.GT.(-elim)) go to 150
138  140 CONTINUE
139  y(nn) = czero
140  nn = nn - 1
141  nuf = nuf + 1
142  IF (nn.EQ.0) RETURN
143  go to 110
144  150 CONTINUE
145  ascle = 1.0e+3*r1mach(1)/tol
146  cz = cz + clog(phi)
147  IF (iform.EQ.1) go to 160
148  cz = cz - cmplx(0.25e0,0.0e0)*clog(arg) - cmplx(aic,0.0e0)
149  160 CONTINUE
150  ax = exp(rcz)/tol
151  ay = aimag(cz)
152  cz = cmplx(ax,0.0e0)*cmplx(cos(ay),sin(ay))
153  CALL cuchk(cz, nw, ascle, tol)
154  IF (nw.EQ.1) go to 140
155  RETURN
156  170 CONTINUE
157  nuf = -1
158  RETURN
159  END