GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cbuni.f
Go to the documentation of this file.
1  SUBROUTINE cbuni(Z, FNU, KODE, N, Y, NZ, NUI, NLAST, FNUL, TOL,
2  * ELIM, ALIM)
3 C***BEGIN PROLOGUE CBUNI
4 C***REFER TO CBESI,CBESK
5 C
6 C CBUNI COMPUTES THE I BESSEL FUNCTION FOR LARGE CABS(Z).GT.
7 C FNUL AND FNU+N-1.LT.FNUL. THE ORDER IS INCREASED FROM
8 C FNU+N-1 GREATER THAN FNUL BY ADDING NUI AND COMPUTING
9 C ACCORDING TO THE UNIFORM ASYMPTOTIC EXPANSION FOR I(FNU,Z)
10 C ON IFORM=1 AND THE EXPANSION FOR J(FNU,Z) ON IFORM=2
11 C
12 C***ROUTINES CALLED CUNI1,CUNI2,R1MACH
13 C***END PROLOGUE CBUNI
14  COMPLEX CSCL, CSCR, CY, RZ, ST, S1, S2, Y, Z
15  REAL ALIM, AX, AY, DFNU, ELIM, FNU, FNUI, FNUL, GNU, TOL, XX, YY,
16  * ascle, bry, str, sti, stm, r1mach
17  INTEGER I, IFLAG, IFORM, K, KODE, N, NL, NLAST, NUI, NW, NZ
18  dimension y(n), cy(2), bry(3)
19  nz = 0
20  xx = real(z)
21  yy = aimag(z)
22  ax = abs(xx)*1.7321e0
23  ay = abs(yy)
24  iform = 1
25  IF (ay.GT.ax) iform = 2
26  IF (nui.EQ.0) GO TO 60
27  fnui = float(nui)
28  dfnu = fnu + float(n-1)
29  gnu = dfnu + fnui
30  IF (iform.EQ.2) GO TO 10
31 C-----------------------------------------------------------------------
32 C ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
33 C -PI/3.LE.ARG(Z).LE.PI/3
34 C-----------------------------------------------------------------------
35  CALL cuni1(z, gnu, kode, 2, cy, nw, nlast, fnul, tol, elim, alim)
36  GO TO 20
37  10 CONTINUE
38 C-----------------------------------------------------------------------
39 C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
40 C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
41 C AND HPI=PI/2
42 C-----------------------------------------------------------------------
43  CALL cuni2(z, gnu, kode, 2, cy, nw, nlast, fnul, tol, elim, alim)
44  20 CONTINUE
45  IF (nw.LT.0) GO TO 50
46  IF (nw.NE.0) GO TO 90
47  ay = cabs(cy(1))
48 C----------------------------------------------------------------------
49 C SCALE BACKWARD RECURRENCE, BRY(3) IS DEFINED BUT NEVER USED
50 C----------------------------------------------------------------------
51  bry(1) = 1.0e+3*r1mach(1)/tol
52  bry(2) = 1.0e0/bry(1)
53  bry(3) = bry(2)
54  iflag = 2
55  ascle = bry(2)
56  ax = 1.0e0
57  cscl = cmplx(ax,0.0e0)
58  IF (ay.GT.bry(1)) GO TO 21
59  iflag = 1
60  ascle = bry(1)
61  ax = 1.0e0/tol
62  cscl = cmplx(ax,0.0e0)
63  GO TO 25
64  21 CONTINUE
65  IF (ay.LT.bry(2)) GO TO 25
66  iflag = 3
67  ascle = bry(3)
68  ax = tol
69  cscl = cmplx(ax,0.0e0)
70  25 CONTINUE
71  ay = 1.0e0/ax
72  cscr = cmplx(ay,0.0e0)
73  s1 = cy(2)*cscl
74  s2 = cy(1)*cscl
75  rz = cmplx(2.0e0,0.0e0)/z
76  DO 30 i=1,nui
77  st = s2
78  s2 = cmplx(dfnu+fnui,0.0e0)*rz*s2 + s1
79  s1 = st
80  fnui = fnui - 1.0e0
81  IF (iflag.GE.3) GO TO 30
82  st = s2*cscr
83  str = real(st)
84  sti = aimag(st)
85  str = abs(str)
86  sti = abs(sti)
87  stm = amax1(str,sti)
88  IF (stm.LE.ascle) GO TO 30
89  iflag = iflag+1
90  ascle = bry(iflag)
91  s1 = s1*cscr
92  s2 = st
93  ax = ax*tol
94  ay = 1.0e0/ax
95  cscl = cmplx(ax,0.0e0)
96  cscr = cmplx(ay,0.0e0)
97  s1 = s1*cscl
98  s2 = s2*cscl
99  30 CONTINUE
100  y(n) = s2*cscr
101  IF (n.EQ.1) RETURN
102  nl = n - 1
103  fnui = float(nl)
104  k = nl
105  DO 40 i=1,nl
106  st = s2
107  s2 = cmplx(fnu+fnui,0.0e0)*rz*s2 + s1
108  s1 = st
109  st = s2*cscr
110  y(k) = st
111  fnui = fnui - 1.0e0
112  k = k - 1
113  IF (iflag.GE.3) GO TO 40
114  str = real(st)
115  sti = aimag(st)
116  str = abs(str)
117  sti = abs(sti)
118  stm = amax1(str,sti)
119  IF (stm.LE.ascle) GO TO 40
120  iflag = iflag+1
121  ascle = bry(iflag)
122  s1 = s1*cscr
123  s2 = st
124  ax = ax*tol
125  ay = 1.0e0/ax
126  cscl = cmplx(ax,0.0e0)
127  cscr = cmplx(ay,0.0e0)
128  s1 = s1*cscl
129  s2 = s2*cscl
130  40 CONTINUE
131  RETURN
132  50 CONTINUE
133  nz = -1
134  IF(nw.EQ.(-2)) nz=-2
135  RETURN
136  60 CONTINUE
137  IF (iform.EQ.2) GO TO 70
138 C-----------------------------------------------------------------------
139 C ASYMPTOTIC EXPANSION FOR I(FNU,Z) FOR LARGE FNU APPLIED IN
140 C -PI/3.LE.ARG(Z).LE.PI/3
141 C-----------------------------------------------------------------------
142  CALL cuni1(z, fnu, kode, n, y, nw, nlast, fnul, tol, elim, alim)
143  GO TO 80
144  70 CONTINUE
145 C-----------------------------------------------------------------------
146 C ASYMPTOTIC EXPANSION FOR J(FNU,Z*EXP(M*HPI)) FOR LARGE FNU
147 C APPLIED IN PI/3.LT.ABS(ARG(Z)).LE.PI/2 WHERE M=+I OR -I
148 C AND HPI=PI/2
149 C-----------------------------------------------------------------------
150  CALL cuni2(z, fnu, kode, n, y, nw, nlast, fnul, tol, elim, alim)
151  80 CONTINUE
152  IF (nw.LT.0) GO TO 50
153  nz = nw
154  RETURN
155  90 CONTINUE
156  nlast = n
157  RETURN
158  END
double complex cmplx
Definition: Faddeeva.cc:217
subroutine cbuni(Z, FNU, KODE, N, Y, NZ, NUI, NLAST, FNUL, TOL, ELIM, ALIM)
Definition: cbuni.f:3
subroutine cuni1(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
Definition: cuni1.f:3
subroutine cuni2(Z, FNU, KODE, N, Y, NZ, NLAST, FNUL, TOL, ELIM, ALIM)
Definition: cuni2.f:3
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
static T abs(T x)
Definition: pr-output.cc:1678
real function r1mach(i)
Definition: r1mach.f:23