GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cuni2.f
Go to the documentation of this file.
1  SUBROUTINE cuni2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM,
2  * ALIM)
3 C***BEGIN PROLOGUE CUNI2
4 C***REFER TO CBESI,CBESK
5 C
6 C CUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
7 C UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
8 C OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
9 C
10 C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
11 C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
12 C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
13 C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
14 C Y(I)=CZERO FOR I=NLAST+1,N
15 C
16 C***ROUTINES CALLED CAIRY,CUCHK,CUNHJ,CUOIK,R1MACH
17 C***END PROLOGUE CUNI2
18  COMPLEX AI, ARG, ASUM, BSUM, CFN, CI, CID, CIP, CONE, CRSC, CSCL,
19  * csr, css, cy, czero, c1, c2, dai, phi, rz, s1, s2, y, z, zb,
20  * zeta1, zeta2, zn, zar
21  REAL AARG, AIC, ALIM, ANG, APHI, ASCLE, AY, BRY, CAR, C2I, C2M,
22  * c2r, elim, fn, fnu, fnul, hpi, rs1, sar, tol, yy, r1mach
23  INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
24  * nn, nuf, nw, nz, idum
25  dimension bry(3), y(n), cip(4), css(3), csr(3), cy(2)
26  DATA czero,cone,ci/(0.0e0,0.0e0),(1.0e0,0.0e0),(0.0e0,1.0e0)/
27  DATA cip(1),cip(2),cip(3),cip(4)/
28  1 (1.0e0,0.0e0), (0.0e0,1.0e0), (-1.0e0,0.0e0), (0.0e0,-1.0e0)/
29  DATA hpi, aic /
30  1 1.57079632679489662e+00, 1.265512123484645396e+00/
31 C
32  nz = 0
33  nd = n
34  nlast = 0
35 C-----------------------------------------------------------------------
36 C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
37 C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
38 C EXP(ALIM)=EXP(ELIM)*TOL
39 C-----------------------------------------------------------------------
40  cscl = cmplx(1.0e0/tol,0.0e0)
41  crsc = cmplx(tol,0.0e0)
42  css(1) = cscl
43  css(2) = cone
44  css(3) = crsc
45  csr(1) = crsc
46  csr(2) = cone
47  csr(3) = cscl
48  bry(1) = 1.0e+3*r1mach(1)/tol
49  yy = aimag(z)
50 C-----------------------------------------------------------------------
51 C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
52 C-----------------------------------------------------------------------
53  zn = -z*ci
54  zb = z
55  cid = -ci
56  inu = int(fnu)
57  ang = hpi*(fnu-float(inu))
58  car = cos(ang)
59  sar = sin(ang)
60  c2 = cmplx(car,sar)
61  zar = c2
62  in = inu + n - 1
63  in = mod(in,4)
64  c2 = c2*cip(in+1)
65  IF (yy.GT.0.0e0) GO TO 10
66  zn = conjg(-zn)
67  zb = conjg(zb)
68  cid = -cid
69  c2 = conjg(c2)
70  10 CONTINUE
71 C-----------------------------------------------------------------------
72 C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
73 C-----------------------------------------------------------------------
74  fn = amax1(fnu,1.0e0)
75  CALL cunhj(zn, fn, 1, tol, phi, arg, zeta1, zeta2, asum, bsum)
76  IF (kode.EQ.1) GO TO 20
77  cfn = cmplx(fnu,0.0e0)
78  s1 = -zeta1 + cfn*(cfn/(zb+zeta2))
79  GO TO 30
80  20 CONTINUE
81  s1 = -zeta1 + zeta2
82  30 CONTINUE
83  rs1 = real(s1)
84  IF (abs(rs1).GT.elim) GO TO 150
85  40 CONTINUE
86  nn = min0(2,nd)
87  DO 90 i=1,nn
88  fn = fnu + float(nd-i)
89  CALL cunhj(zn, fn, 0, tol, phi, arg, zeta1, zeta2, asum, bsum)
90  IF (kode.EQ.1) GO TO 50
91  cfn = cmplx(fn,0.0e0)
92  ay = abs(yy)
93  s1 = -zeta1 + cfn*(cfn/(zb+zeta2)) + cmplx(0.0e0,ay)
94  GO TO 60
95  50 CONTINUE
96  s1 = -zeta1 + zeta2
97  60 CONTINUE
98 C-----------------------------------------------------------------------
99 C TEST FOR UNDERFLOW AND OVERFLOW
100 C-----------------------------------------------------------------------
101  rs1 = real(s1)
102  IF (abs(rs1).GT.elim) GO TO 120
103  IF (i.EQ.1) iflag = 2
104  IF (abs(rs1).LT.alim) GO TO 70
105 C-----------------------------------------------------------------------
106 C REFINE TEST AND SCALE
107 C-----------------------------------------------------------------------
108 C-----------------------------------------------------------------------
109  aphi = cabs(phi)
110  aarg = cabs(arg)
111  rs1 = rs1 + alog(aphi) - 0.25e0*alog(aarg) - aic
112  IF (abs(rs1).GT.elim) GO TO 120
113  IF (i.EQ.1) iflag = 1
114  IF (rs1.LT.0.0e0) GO TO 70
115  IF (i.EQ.1) iflag = 3
116  70 CONTINUE
117 C-----------------------------------------------------------------------
118 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
119 C EXPONENT EXTREMES
120 C-----------------------------------------------------------------------
121  CALL cairy(arg, 0, 2, ai, nai, idum)
122  CALL cairy(arg, 1, 2, dai, ndai, idum)
123  s2 = phi*(ai*asum+dai*bsum)
124  c2r = real(s1)
125  c2i = aimag(s1)
126  c2m = exp(c2r)*real(css(iflag))
127  s1 = cmplx(c2m,0.0e0)*cmplx(cos(c2i),sin(c2i))
128  s2 = s2*s1
129  IF (iflag.NE.1) GO TO 80
130  CALL cuchk(s2, nw, bry(1), tol)
131  IF (nw.NE.0) GO TO 120
132  80 CONTINUE
133  IF (yy.LE.0.0e0) s2 = conjg(s2)
134  j = nd - i + 1
135  s2 = s2*c2
136  cy(i) = s2
137  y(j) = s2*csr(iflag)
138  c2 = c2*cid
139  90 CONTINUE
140  IF (nd.LE.2) GO TO 110
141  rz = cmplx(2.0e0,0.0e0)/z
142  bry(2) = 1.0e0/bry(1)
143  bry(3) = r1mach(2)
144  s1 = cy(1)
145  s2 = cy(2)
146  c1 = csr(iflag)
147  ascle = bry(iflag)
148  k = nd - 2
149  fn = float(k)
150  DO 100 i=3,nd
151  c2 = s2
152  s2 = s1 + cmplx(fnu+fn,0.0e0)*rz*s2
153  s1 = c2
154  c2 = s2*c1
155  y(k) = c2
156  k = k - 1
157  fn = fn - 1.0e0
158  IF (iflag.GE.3) GO TO 100
159  c2r = real(c2)
160  c2i = aimag(c2)
161  c2r = abs(c2r)
162  c2i = abs(c2i)
163  c2m = amax1(c2r,c2i)
164  IF (c2m.LE.ascle) GO TO 100
165  iflag = iflag + 1
166  ascle = bry(iflag)
167  s1 = s1*c1
168  s2 = c2
169  s1 = s1*css(iflag)
170  s2 = s2*css(iflag)
171  c1 = csr(iflag)
172  100 CONTINUE
173  110 CONTINUE
174  RETURN
175  120 CONTINUE
176  IF (rs1.GT.0.0e0) GO TO 140
177 C-----------------------------------------------------------------------
178 C SET UNDERFLOW AND UPDATE PARAMETERS
179 C-----------------------------------------------------------------------
180  y(nd) = czero
181  nz = nz + 1
182  nd = nd - 1
183  IF (nd.EQ.0) GO TO 110
184  CALL cuoik(z, fnu, kode, 1, nd, y, nuf, tol, elim, alim)
185  IF (nuf.LT.0) GO TO 140
186  nd = nd - nuf
187  nz = nz + nuf
188  IF (nd.EQ.0) GO TO 110
189  fn = fnu + float(nd-1)
190  IF (fn.LT.fnul) GO TO 130
191 C FN = AIMAG(CID)
192 C J = NUF + 1
193 C K = MOD(J,4) + 1
194 C S1 = CIP(K)
195 C IF (FN.LT.0.0E0) S1 = CONJG(S1)
196 C C2 = C2*S1
197  in = inu + nd - 1
198  in = mod(in,4) + 1
199  c2 = zar*cip(in)
200  IF (yy.LE.0.0e0)c2=conjg(c2)
201  GO TO 40
202  130 CONTINUE
203  nlast = nd
204  RETURN
205  140 CONTINUE
206  nz = -1
207  RETURN
208  150 CONTINUE
209  IF (rs1.GT.0.0e0) GO TO 140
210  nz = n
211  DO 160 i=1,n
212  y(i) = czero
213  160 CONTINUE
214  RETURN
215  END
cmplx
double complex cmplx
Definition: Faddeeva.cc:230
cairy
subroutine cairy(Z, ID, KODE, AI, NZ, IERR)
Definition: cairy.f:2
cuchk
subroutine cuchk(Y, NZ, ASCLE, TOL)
Definition: cuchk.f:2
cunhj
subroutine cunhj(Z, FNU, IPMTR, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
Definition: cunhj.f:3
cuni2
subroutine cuni2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
Definition: cuni2.f:3
cuoik
subroutine cuoik(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM)
Definition: cuoik.f:2
real
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
mod
T mod(T x, T y)
Definition: lo-mappers.h:294
r1mach
real function r1mach(i)
Definition: r1mach.f:23