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
cuni1.f
Go to the documentation of this file.
1  SUBROUTINE cuni1(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM,
2  * alim)
3 C***BEGIN PROLOGUE CUNI1
4 C***REFER TO CBESI,CBESK
5 C
6 C CUNI1 COMPUTES I(FNU,Z) BY MEANS OF THE UNIFORM ASYMPTOTIC
7 C EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.
8 C
9 C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
10 C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
11 C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
12 C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
13 C Y(I)=CZERO FOR I=NLAST+1,N
14 C
15 C***ROUTINES CALLED CUCHK,CUNIK,CUOIK,R1MACH
16 C***END PROLOGUE CUNI1
17  COMPLEX cfn, cone, crsc, cscl, csr, css, cwrk, czero, c1, c2,
18  * phi, rz, sum, s1, s2, y, z, zeta1, zeta2, cy
19  REAL alim, aphi, ascle, bry, c2i, c2m, c2r, elim, fn, fnu, fnul,
20  * rs1, tol, yy, r1mach
21  INTEGER i, iflag, init, k, kode, m, n, nd, nlast, nn, nuf, nw, nz
22  dimension bry(3), y(n), cwrk(16), css(3), csr(3), cy(2)
23  DATA czero, cone / (0.0e0,0.0e0), (1.0e0,0.0e0) /
24 C
25  nz = 0
26  nd = n
27  nlast = 0
28 C-----------------------------------------------------------------------
29 C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
30 C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
31 C EXP(ALIM)=EXP(ELIM)*TOL
32 C-----------------------------------------------------------------------
33  cscl = cmplx(1.0e0/tol,0.0e0)
34  crsc = cmplx(tol,0.0e0)
35  css(1) = cscl
36  css(2) = cone
37  css(3) = crsc
38  csr(1) = crsc
39  csr(2) = cone
40  csr(3) = cscl
41  bry(1) = 1.0e+3*r1mach(1)/tol
42 C-----------------------------------------------------------------------
43 C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
44 C-----------------------------------------------------------------------
45  fn = amax1(fnu,1.0e0)
46  init = 0
47  CALL cunik(z, fn, 1, 1, tol, init, phi, zeta1, zeta2, sum, cwrk)
48  IF (kode.EQ.1) go to 10
49  cfn = cmplx(fn,0.0e0)
50  s1 = -zeta1 + cfn*(cfn/(z+zeta2))
51  go to 20
52  10 CONTINUE
53  s1 = -zeta1 + zeta2
54  20 CONTINUE
55  rs1 = REAL(s1)
56  IF (abs(rs1).GT.elim) go to 130
57  30 CONTINUE
58  nn = min0(2,nd)
59  DO 80 i=1,nn
60  fn = fnu + float(nd-i)
61  init = 0
62  CALL cunik(z, fn, 1, 0, tol, init, phi, zeta1, zeta2, sum, cwrk)
63  IF (kode.EQ.1) go to 40
64  cfn = cmplx(fn,0.0e0)
65  yy = aimag(z)
66  s1 = -zeta1 + cfn*(cfn/(z+zeta2)) + cmplx(0.0e0,yy)
67  go to 50
68  40 CONTINUE
69  s1 = -zeta1 + zeta2
70  50 CONTINUE
71 C-----------------------------------------------------------------------
72 C TEST FOR UNDERFLOW AND OVERFLOW
73 C-----------------------------------------------------------------------
74  rs1 = REAL(s1)
75  IF (abs(rs1).GT.elim) go to 110
76  IF (i.EQ.1) iflag = 2
77  IF (abs(rs1).LT.alim) go to 60
78 C-----------------------------------------------------------------------
79 C REFINE TEST AND SCALE
80 C-----------------------------------------------------------------------
81  aphi = cabs(phi)
82  rs1 = rs1 + alog(aphi)
83  IF (abs(rs1).GT.elim) go to 110
84  IF (i.EQ.1) iflag = 1
85  IF (rs1.LT.0.0e0) go to 60
86  IF (i.EQ.1) iflag = 3
87  60 CONTINUE
88 C-----------------------------------------------------------------------
89 C SCALE S1 IF CABS(S1).LT.ASCLE
90 C-----------------------------------------------------------------------
91  s2 = phi*sum
92  c2r = REAL(s1)
93  c2i = aimag(s1)
94  c2m = exp(c2r)*REAL(css(iflag))
95  s1 = cmplx(c2m,0.0e0)*cmplx(cos(c2i),sin(c2i))
96  s2 = s2*s1
97  IF (iflag.NE.1) go to 70
98  CALL cuchk(s2, nw, bry(1), tol)
99  IF (nw.NE.0) go to 110
100  70 CONTINUE
101  m = nd - i + 1
102  cy(i) = s2
103  y(m) = s2*csr(iflag)
104  80 CONTINUE
105  IF (nd.LE.2) go to 100
106  rz = cmplx(2.0e0,0.0e0)/z
107  bry(2) = 1.0e0/bry(1)
108  bry(3) = r1mach(2)
109  s1 = cy(1)
110  s2 = cy(2)
111  c1 = csr(iflag)
112  ascle = bry(iflag)
113  k = nd - 2
114  fn = float(k)
115  DO 90 i=3,nd
116  c2 = s2
117  s2 = s1 + cmplx(fnu+fn,0.0e0)*rz*s2
118  s1 = c2
119  c2 = s2*c1
120  y(k) = c2
121  k = k - 1
122  fn = fn - 1.0e0
123  IF (iflag.GE.3) go to 90
124  c2r = REAL(c2)
125  c2i = aimag(c2)
126  c2r = abs(c2r)
127  c2i = abs(c2i)
128  c2m = amax1(c2r,c2i)
129  IF (c2m.LE.ascle) go to 90
130  iflag = iflag + 1
131  ascle = bry(iflag)
132  s1 = s1*c1
133  s2 = c2
134  s1 = s1*css(iflag)
135  s2 = s2*css(iflag)
136  c1 = csr(iflag)
137  90 CONTINUE
138  100 CONTINUE
139  RETURN
140 C-----------------------------------------------------------------------
141 C SET UNDERFLOW AND UPDATE PARAMETERS
142 C-----------------------------------------------------------------------
143  110 CONTINUE
144  IF (rs1.GT.0.0e0) go to 120
145  y(nd) = czero
146  nz = nz + 1
147  nd = nd - 1
148  IF (nd.EQ.0) go to 100
149  CALL cuoik(z, fnu, kode, 1, nd, y, nuf, tol, elim, alim)
150  IF (nuf.LT.0) go to 120
151  nd = nd - nuf
152  nz = nz + nuf
153  IF (nd.EQ.0) go to 100
154  fn = fnu + float(nd-1)
155  IF (fn.GE.fnul) go to 30
156  nlast = nd
157  RETURN
158  120 CONTINUE
159  nz = -1
160  RETURN
161  130 CONTINUE
162  IF (rs1.GT.0.0e0) go to 120
163  nz = n
164  DO 140 i=1,n
165  y(i) = czero
166  140 CONTINUE
167  RETURN
168  END