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
zbinu.f
Go to the documentation of this file.
1  SUBROUTINE zbinu(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL,
2  * tol, elim, alim)
3 C***BEGIN PROLOGUE ZBINU
4 C***REFER TO ZBESH,ZBESI,ZBESJ,ZBESK,ZAIRY,ZBIRY
5 C
6 C ZBINU COMPUTES THE I FUNCTION IN THE RIGHT HALF Z PLANE
7 C
8 C***ROUTINES CALLED XZABS,ZASYI,ZBUNI,ZMLRI,ZSERI,ZUOIK,ZWRSK
9 C***END PROLOGUE ZBINU
10  DOUBLE PRECISION alim, az, cwi, cwr, cyi, cyr, dfnu, elim, fnu,
11  * fnul, rl, tol, zeroi, zeror, zi, zr, xzabs
12  INTEGER i, inw, kode, n, nlast, nn, nui, nw, nz
13  dimension cyr(n), cyi(n), cwr(2), cwi(2)
14  DATA zeror,zeroi / 0.0d0, 0.0d0 /
15 C
16  nz = 0
17  az = xzabs(zr,zi)
18  nn = n
19  dfnu = fnu + dble(float(n-1))
20  IF (az.LE.2.0d0) go to 10
21  IF (az*az*0.25d0.GT.dfnu+1.0d0) go to 20
22  10 CONTINUE
23 C-----------------------------------------------------------------------
24 C POWER SERIES
25 C-----------------------------------------------------------------------
26  CALL zseri(zr, zi, fnu, kode, nn, cyr, cyi, nw, tol, elim, alim)
27  inw = iabs(nw)
28  nz = nz + inw
29  nn = nn - inw
30  IF (nn.EQ.0) RETURN
31  IF (nw.GE.0) go to 120
32  dfnu = fnu + dble(float(nn-1))
33  20 CONTINUE
34  IF (az.LT.rl) go to 40
35  IF (dfnu.LE.1.0d0) go to 30
36  IF (az+az.LT.dfnu*dfnu) go to 50
37 C-----------------------------------------------------------------------
38 C ASYMPTOTIC EXPANSION FOR LARGE Z
39 C-----------------------------------------------------------------------
40  30 CONTINUE
41  CALL zasyi(zr, zi, fnu, kode, nn, cyr, cyi, nw, rl, tol, elim,
42  * alim)
43  IF (nw.LT.0) go to 130
44  go to 120
45  40 CONTINUE
46  IF (dfnu.LE.1.0d0) go to 70
47  50 CONTINUE
48 C-----------------------------------------------------------------------
49 C OVERFLOW AND UNDERFLOW TEST ON I SEQUENCE FOR MILLER ALGORITHM
50 C-----------------------------------------------------------------------
51  CALL zuoik(zr, zi, fnu, kode, 1, nn, cyr, cyi, nw, tol, elim,
52  * alim)
53  IF (nw.LT.0) go to 130
54  nz = nz + nw
55  nn = nn - nw
56  IF (nn.EQ.0) RETURN
57  dfnu = fnu+dble(float(nn-1))
58  IF (dfnu.GT.fnul) go to 110
59  IF (az.GT.fnul) go to 110
60  60 CONTINUE
61  IF (az.GT.rl) go to 80
62  70 CONTINUE
63 C-----------------------------------------------------------------------
64 C MILLER ALGORITHM NORMALIZED BY THE SERIES
65 C-----------------------------------------------------------------------
66  CALL zmlri(zr, zi, fnu, kode, nn, cyr, cyi, nw, tol)
67  IF(nw.LT.0) go to 130
68  go to 120
69  80 CONTINUE
70 C-----------------------------------------------------------------------
71 C MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN
72 C-----------------------------------------------------------------------
73 C-----------------------------------------------------------------------
74 C OVERFLOW TEST ON K FUNCTIONS USED IN WRONSKIAN
75 C-----------------------------------------------------------------------
76  CALL zuoik(zr, zi, fnu, kode, 2, 2, cwr, cwi, nw, tol, elim,
77  * alim)
78  IF (nw.GE.0) go to 100
79  nz = nn
80  DO 90 i=1,nn
81  cyr(i) = zeror
82  cyi(i) = zeroi
83  90 CONTINUE
84  RETURN
85  100 CONTINUE
86  IF (nw.GT.0) go to 130
87  CALL zwrsk(zr, zi, fnu, kode, nn, cyr, cyi, nw, cwr, cwi, tol,
88  * elim, alim)
89  IF (nw.LT.0) go to 130
90  go to 120
91  110 CONTINUE
92 C-----------------------------------------------------------------------
93 C INCREMENT FNU+NN-1 UP TO FNUL, COMPUTE AND RECUR BACKWARD
94 C-----------------------------------------------------------------------
95  nui = int(sngl(fnul-dfnu)) + 1
96  nui = max0(nui,0)
97  CALL zbuni(zr, zi, fnu, kode, nn, cyr, cyi, nw, nui, nlast, fnul,
98  * tol, elim, alim)
99  IF (nw.LT.0) go to 130
100  nz = nz + nw
101  IF (nlast.EQ.0) go to 120
102  nn = nlast
103  go to 60
104  120 CONTINUE
105  RETURN
106  130 CONTINUE
107  nz = -1
108  IF(nw.EQ.(-2)) nz=-2
109  RETURN
110  END