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
zuni2.f
Go to the documentation of this file.
1  SUBROUTINE zuni2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
2  * tol, elim, alim)
3 C***BEGIN PROLOGUE ZUNI2
4 C***REFER TO ZBESI,ZBESK
5 C
6 C ZUNI2 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 ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,XZABS
17 C***END PROLOGUE ZUNI2
18 C COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
19 C *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN
20  DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI,
21  * argr, ascle, asumi, asumr, bry, bsumi, bsumr, cidi, cipi, cipr,
22  * coner, crsc, cscl, csrr, cssr, c1r, c2i, c2m, c2r, daii,
23  * dair, elim, fn, fnu, fnul, hpi, phii, phir, rast, raz, rs1, rzi,
24  * rzr, sti, str, s1i, s1r, s2i, s2r, tol, yi, yr, zbi, zbr, zeroi,
25  * zeror, zeta1i, zeta1r, zeta2i, zeta2r, zi, zni, znr, zr, cyr,
26  * cyi, d1mach, xzabs, car, sar
27  INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
28  * nn, nuf, nw, nz, idum
29  dimension bry(3), yr(n), yi(n), cipr(4), cipi(4), cssr(3),
30  * csrr(3), cyr(2), cyi(2)
31  DATA zeror,zeroi,coner / 0.0d0, 0.0d0, 1.0d0 /
32  DATA cipr(1),cipi(1),cipr(2),cipi(2),cipr(3),cipi(3),cipr(4),
33  * cipi(4)/ 1.0d0,0.0d0, 0.0d0,1.0d0, -1.0d0,0.0d0, 0.0d0,-1.0d0/
34  DATA hpi, aic /
35  1 1.57079632679489662d+00, 1.265512123484645396d+00/
36 C
37  nz = 0
38  nd = n
39  nlast = 0
40 C-----------------------------------------------------------------------
41 C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
42 C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
43 C EXP(ALIM)=EXP(ELIM)*TOL
44 C-----------------------------------------------------------------------
45  cscl = 1.0d0/tol
46  crsc = tol
47  cssr(1) = cscl
48  cssr(2) = coner
49  cssr(3) = crsc
50  csrr(1) = crsc
51  csrr(2) = coner
52  csrr(3) = cscl
53  bry(1) = 1.0d+3*d1mach(1)/tol
54 C-----------------------------------------------------------------------
55 C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
56 C-----------------------------------------------------------------------
57  znr = zi
58  zni = -zr
59  zbr = zr
60  zbi = zi
61  cidi = -coner
62  inu = int(sngl(fnu))
63  ang = hpi*(fnu-dble(float(inu)))
64  c2r = dcos(ang)
65  c2i = dsin(ang)
66  car = c2r
67  sar = c2i
68  in = inu + n - 1
69  in = mod(in,4) + 1
70  str = c2r*cipr(in) - c2i*cipi(in)
71  c2i = c2r*cipi(in) + c2i*cipr(in)
72  c2r = str
73  IF (zi.GT.0.0d0) go to 10
74  znr = -znr
75  zbi = -zbi
76  cidi = -cidi
77  c2i = -c2i
78  10 CONTINUE
79 C-----------------------------------------------------------------------
80 C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
81 C-----------------------------------------------------------------------
82  fn = dmax1(fnu,1.0d0)
83  CALL zunhj(znr, zni, fn, 1, tol, phir, phii, argr, argi, zeta1r,
84  * zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
85  IF (kode.EQ.1) go to 20
86  str = zbr + zeta2r
87  sti = zbi + zeta2i
88  rast = fn/xzabs(str,sti)
89  str = str*rast*rast
90  sti = -sti*rast*rast
91  s1r = -zeta1r + str
92  s1i = -zeta1i + sti
93  go to 30
94  20 CONTINUE
95  s1r = -zeta1r + zeta2r
96  s1i = -zeta1i + zeta2i
97  30 CONTINUE
98  rs1 = s1r
99  IF (dabs(rs1).GT.elim) go to 150
100  40 CONTINUE
101  nn = min0(2,nd)
102  DO 90 i=1,nn
103  fn = fnu + dble(float(nd-i))
104  CALL zunhj(znr, zni, fn, 0, tol, phir, phii, argr, argi,
105  * zeta1r, zeta1i, zeta2r, zeta2i, asumr, asumi, bsumr, bsumi)
106  IF (kode.EQ.1) go to 50
107  str = zbr + zeta2r
108  sti = zbi + zeta2i
109  rast = fn/xzabs(str,sti)
110  str = str*rast*rast
111  sti = -sti*rast*rast
112  s1r = -zeta1r + str
113  s1i = -zeta1i + sti + dabs(zi)
114  go to 60
115  50 CONTINUE
116  s1r = -zeta1r + zeta2r
117  s1i = -zeta1i + zeta2i
118  60 CONTINUE
119 C-----------------------------------------------------------------------
120 C TEST FOR UNDERFLOW AND OVERFLOW
121 C-----------------------------------------------------------------------
122  rs1 = s1r
123  IF (dabs(rs1).GT.elim) go to 120
124  IF (i.EQ.1) iflag = 2
125  IF (dabs(rs1).LT.alim) go to 70
126 C-----------------------------------------------------------------------
127 C REFINE TEST AND SCALE
128 C-----------------------------------------------------------------------
129 C-----------------------------------------------------------------------
130  aphi = xzabs(phir,phii)
131  aarg = xzabs(argr,argi)
132  rs1 = rs1 + dlog(aphi) - 0.25d0*dlog(aarg) - aic
133  IF (dabs(rs1).GT.elim) go to 120
134  IF (i.EQ.1) iflag = 1
135  IF (rs1.LT.0.0d0) go to 70
136  IF (i.EQ.1) iflag = 3
137  70 CONTINUE
138 C-----------------------------------------------------------------------
139 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
140 C EXPONENT EXTREMES
141 C-----------------------------------------------------------------------
142  CALL zairy(argr, argi, 0, 2, air, aii, nai, idum)
143  CALL zairy(argr, argi, 1, 2, dair, daii, ndai, idum)
144  str = dair*bsumr - daii*bsumi
145  sti = dair*bsumi + daii*bsumr
146  str = str + (air*asumr-aii*asumi)
147  sti = sti + (air*asumi+aii*asumr)
148  s2r = phir*str - phii*sti
149  s2i = phir*sti + phii*str
150  str = dexp(s1r)*cssr(iflag)
151  s1r = str*dcos(s1i)
152  s1i = str*dsin(s1i)
153  str = s2r*s1r - s2i*s1i
154  s2i = s2r*s1i + s2i*s1r
155  s2r = str
156  IF (iflag.NE.1) go to 80
157  CALL zuchk(s2r, s2i, nw, bry(1), tol)
158  IF (nw.NE.0) go to 120
159  80 CONTINUE
160  IF (zi.LE.0.0d0) s2i = -s2i
161  str = s2r*c2r - s2i*c2i
162  s2i = s2r*c2i + s2i*c2r
163  s2r = str
164  cyr(i) = s2r
165  cyi(i) = s2i
166  j = nd - i + 1
167  yr(j) = s2r*csrr(iflag)
168  yi(j) = s2i*csrr(iflag)
169  str = -c2i*cidi
170  c2i = c2r*cidi
171  c2r = str
172  90 CONTINUE
173  IF (nd.LE.2) go to 110
174  raz = 1.0d0/xzabs(zr,zi)
175  str = zr*raz
176  sti = -zi*raz
177  rzr = (str+str)*raz
178  rzi = (sti+sti)*raz
179  bry(2) = 1.0d0/bry(1)
180  bry(3) = d1mach(2)
181  s1r = cyr(1)
182  s1i = cyi(1)
183  s2r = cyr(2)
184  s2i = cyi(2)
185  c1r = csrr(iflag)
186  ascle = bry(iflag)
187  k = nd - 2
188  fn = dble(float(k))
189  DO 100 i=3,nd
190  c2r = s2r
191  c2i = s2i
192  s2r = s1r + (fnu+fn)*(rzr*c2r-rzi*c2i)
193  s2i = s1i + (fnu+fn)*(rzr*c2i+rzi*c2r)
194  s1r = c2r
195  s1i = c2i
196  c2r = s2r*c1r
197  c2i = s2i*c1r
198  yr(k) = c2r
199  yi(k) = c2i
200  k = k - 1
201  fn = fn - 1.0d0
202  IF (iflag.GE.3) go to 100
203  str = dabs(c2r)
204  sti = dabs(c2i)
205  c2m = dmax1(str,sti)
206  IF (c2m.LE.ascle) go to 100
207  iflag = iflag + 1
208  ascle = bry(iflag)
209  s1r = s1r*c1r
210  s1i = s1i*c1r
211  s2r = c2r
212  s2i = c2i
213  s1r = s1r*cssr(iflag)
214  s1i = s1i*cssr(iflag)
215  s2r = s2r*cssr(iflag)
216  s2i = s2i*cssr(iflag)
217  c1r = csrr(iflag)
218  100 CONTINUE
219  110 CONTINUE
220  RETURN
221  120 CONTINUE
222  IF (rs1.GT.0.0d0) go to 140
223 C-----------------------------------------------------------------------
224 C SET UNDERFLOW AND UPDATE PARAMETERS
225 C-----------------------------------------------------------------------
226  yr(nd) = zeror
227  yi(nd) = zeroi
228  nz = nz + 1
229  nd = nd - 1
230  IF (nd.EQ.0) go to 110
231  CALL zuoik(zr, zi, fnu, kode, 1, nd, yr, yi, nuf, tol, elim, alim)
232  IF (nuf.LT.0) go to 140
233  nd = nd - nuf
234  nz = nz + nuf
235  IF (nd.EQ.0) go to 110
236  fn = fnu + dble(float(nd-1))
237  IF (fn.LT.fnul) go to 130
238 C FN = CIDI
239 C J = NUF + 1
240 C K = MOD(J,4) + 1
241 C S1R = CIPR(K)
242 C S1I = CIPI(K)
243 C IF (FN.LT.0.0D0) S1I = -S1I
244 C STR = C2R*S1R - C2I*S1I
245 C C2I = C2R*S1I + C2I*S1R
246 C C2R = STR
247  in = inu + nd - 1
248  in = mod(in,4) + 1
249  c2r = car*cipr(in) - sar*cipi(in)
250  c2i = car*cipi(in) + sar*cipr(in)
251  IF (zi.LE.0.0d0) c2i = -c2i
252  go to 40
253  130 CONTINUE
254  nlast = nd
255  RETURN
256  140 CONTINUE
257  nz = -1
258  RETURN
259  150 CONTINUE
260  IF (rs1.GT.0.0d0) go to 140
261  nz = n
262  DO 160 i=1,n
263  yr(i) = zeror
264  yi(i) = zeroi
265  160 CONTINUE
266  RETURN
267  END
dim_vector::str
std::string str(char sep= 'x') const
Definition: dim-vector.cc:63
zairy
subroutine zairy(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
Definition: zairy.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
zuoik
subroutine zuoik(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL, ELIM, ALIM)
Definition: zuoik.f:1
d1mach
double precision function d1mach(i)
Definition: d1mach.f:1
d
F77_RET_T const double const double double * d
Definition: __pchip_deriv__.cc:38
zuni2
subroutine zuni2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
Definition: zuni2.f:1
zunhj
subroutine zunhj(ZR, ZI, FNU, IPMTR, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
Definition: zunhj.f:1
zuchk
subroutine zuchk(YR, YI, NZ, ASCLE, TOL)
Definition: zuchk.f:1
Generated on Thu Jun 4 2015 23:30:23 for GNU Octave by   doxygen 1.8.8