GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
double complex cmplx
Definition: Faddeeva.cc:217
subroutine cuchk(Y, NZ, ASCLE, TOL)
Definition: cuchk.f:2
subroutine cuni1(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
Definition: cuni1.f:3
subroutine cunik(ZR, FNU, IKFLG, IPMTR, TOL, INIT, PHI, ZETA1, ZETA2, SUM, CWRK)
Definition: cunik.f:3
subroutine cuoik(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM)
Definition: cuoik.f:2
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
static T abs(T x)
Definition: pr-output.cc:1678
real function r1mach(i)
Definition: r1mach.f:23