GNU Octave  4.0.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
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
cairy
subroutine cairy(Z, ID, KODE, AI, NZ, IERR)
Definition: cairy.f:1
workspace_element::dimension
std::string dimension(void) const
Definition: workspace-element.h:74
mod
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:959
cuni2
subroutine cuni2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
Definition: cuni2.f:1
base_det::exp
int exp(void) const
Definition: DET.h:64
cuchk
subroutine cuchk(Y, NZ, ASCLE, TOL)
Definition: cuchk.f:1
cmplx
double complex cmplx
Definition: Faddeeva.cc:226
octave_value::sin
octave_value sin(void) const
Definition: ov.h:1198
cuoik
subroutine cuoik(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM)
Definition: cuoik.f:1
cunhj
subroutine cunhj(Z, FNU, IPMTR, TOL, PHI, ARG, ZETA1, ZETA2, ASUM, BSUM)
Definition: cunhj.f:1
octave_value::cos
octave_value cos(void) const
Definition: ov.h:1170
abs
T abs(T x)
Definition: pr-output.cc:3062
Generated on Thu Jun 4 2015 23:30:23 for GNU Octave by   doxygen 1.8.8