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
zkscl.f
Go to the documentation of this file.
1  SUBROUTINE zkscl(ZRR,ZRI,FNU,N,YR,YI,NZ,RZR,RZI,ASCLE,TOL,ELIM)
2 C***BEGIN PROLOGUE ZKSCL
3 C***REFER TO ZBESK
4 C
5 C SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
6 C ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
7 C RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
8 C
9 C***ROUTINES CALLED ZUCHK,XZABS,XZLOG
10 C***END PROLOGUE ZKSCL
11 C COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM
12  DOUBLE PRECISION acs, as, ascle, cki, ckr, csi, csr, cyi,
13  * cyr, elim, fn, fnu, rzi, rzr, str, s1i, s1r, s2i,
14  * s2r, tol, yi, yr, zeroi, zeror, zri, zrr, xzabs,
15  * zdr, zdi, celmr, elm, helim, alas
16  INTEGER i, ic, idum, kk, n, nn, nw, nz
17  dimension yr(n), yi(n), cyr(2), cyi(2)
18  DATA zeror,zeroi / 0.0d0 , 0.0d0 /
19 C
20  nz = 0
21  ic = 0
22  nn = min0(2,n)
23  DO 10 i=1,nn
24  s1r = yr(i)
25  s1i = yi(i)
26  cyr(i) = s1r
27  cyi(i) = s1i
28  as = xzabs(s1r,s1i)
29  acs = -zrr + dlog(as)
30  nz = nz + 1
31  yr(i) = zeror
32  yi(i) = zeroi
33  IF (acs.LT.(-elim)) go to 10
34  CALL xzlog(s1r, s1i, csr, csi, idum)
35  csr = csr - zrr
36  csi = csi - zri
37  str = dexp(csr)/tol
38  csr = str*dcos(csi)
39  csi = str*dsin(csi)
40  CALL zuchk(csr, csi, nw, ascle, tol)
41  IF (nw.NE.0) go to 10
42  yr(i) = csr
43  yi(i) = csi
44  ic = i
45  nz = nz - 1
46  10 CONTINUE
47  IF (n.EQ.1) RETURN
48  IF (ic.GT.1) go to 20
49  yr(1) = zeror
50  yi(1) = zeroi
51  nz = 2
52  20 CONTINUE
53  IF (n.EQ.2) RETURN
54  IF (nz.EQ.0) RETURN
55  fn = fnu + 1.0d0
56  ckr = fn*rzr
57  cki = fn*rzi
58  s1r = cyr(1)
59  s1i = cyi(1)
60  s2r = cyr(2)
61  s2i = cyi(2)
62  helim = 0.5d0*elim
63  elm = dexp(-elim)
64  celmr = elm
65  zdr = zrr
66  zdi = zri
67 C
68 C FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
69 C S2 GETS LARGER THAN EXP(ELIM/2)
70 C
71  DO 30 i=3,n
72  kk = i
73  csr = s2r
74  csi = s2i
75  s2r = ckr*csr - cki*csi + s1r
76  s2i = cki*csr + ckr*csi + s1i
77  s1r = csr
78  s1i = csi
79  ckr = ckr + rzr
80  cki = cki + rzi
81  as = xzabs(s2r,s2i)
82  alas = dlog(as)
83  acs = -zdr + alas
84  nz = nz + 1
85  yr(i) = zeror
86  yi(i) = zeroi
87  IF (acs.LT.(-elim)) go to 25
88  CALL xzlog(s2r, s2i, csr, csi, idum)
89  csr = csr - zdr
90  csi = csi - zdi
91  str = dexp(csr)/tol
92  csr = str*dcos(csi)
93  csi = str*dsin(csi)
94  CALL zuchk(csr, csi, nw, ascle, tol)
95  IF (nw.NE.0) go to 25
96  yr(i) = csr
97  yi(i) = csi
98  nz = nz - 1
99  IF (ic.EQ.kk-1) go to 40
100  ic = kk
101  go to 30
102  25 CONTINUE
103  IF(alas.LT.helim) go to 30
104  zdr = zdr - elim
105  s1r = s1r*celmr
106  s1i = s1i*celmr
107  s2r = s2r*celmr
108  s2i = s2i*celmr
109  30 CONTINUE
110  nz = n
111  IF(ic.EQ.n) nz=n-1
112  go to 45
113  40 CONTINUE
114  nz = kk - 2
115  45 CONTINUE
116  DO 50 i=1,nz
117  yr(i) = zeror
118  yi(i) = zeroi
119  50 CONTINUE
120  RETURN
121  END