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
cbesi.f
Go to the documentation of this file.
1  SUBROUTINE cbesi(Z, FNU, KODE, N, CY, NZ, IERR)
2 C***BEGIN PROLOGUE CBESI
3 C***DATE WRITTEN 830501 (YYMMDD)
4 C***REVISION DATE 890801 (YYMMDD)
5 C***CATEGORY NO. B5K
6 C***KEYWORDS I-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION,
7 C MODIFIED BESSEL FUNCTION OF THE FIRST KIND
8 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
9 C***PURPOSE TO COMPUTE I-BESSEL FUNCTIONS OF COMPLEX ARGUMENT
10 C***DESCRIPTION
11 C
12 C ON KODE=1, CBESI COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
13 C BESSEL FUNCTIONS CY(J)=I(FNU+J-1,Z) FOR REAL, NONNEGATIVE
14 C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z IN THE CUT PLANE
15 C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESI RETURNS THE SCALED
16 C FUNCTIONS
17 C
18 C CY(J)=EXP(-ABS(X))*I(FNU+J-1,Z) J = 1,...,N , X=REAL(Z)
19 C
20 C WITH THE EXPONENTIAL GROWTH REMOVED IN BOTH THE LEFT AND
21 C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND
22 C NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL
23 C FUNCTIONS (REF.1)
24 C
25 C INPUT
26 C Z - Z=CMPLX(X,Y), -PI.LT.ARG(Z).LE.PI
27 C FNU - ORDER OF INITIAL I FUNCTION, FNU.GE.0.0E0
28 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
29 C KODE= 1 RETURNS
30 C CY(J)=I(FNU+J-1,Z), J=1,...,N
31 C = 2 RETURNS
32 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)), J=1,...,N
33 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
34 C
35 C OUTPUT
36 C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
37 C VALUES FOR THE SEQUENCE
38 C CY(J)=I(FNU+J-1,Z) OR
39 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)) J=1,...,N
40 C DEPENDING ON KODE, X=REAL(Z)
41 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
42 C NZ= 0 , NORMAL RETURN
43 C NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO
44 C DUE TO UNDERFLOW, CY(J)=CMPLX(0.0,0.0),
45 C J = N-NZ+1,...,N
46 C IERR - ERROR FLAG
47 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
48 C IERR=1, INPUT ERROR - NO COMPUTATION
49 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z) TOO
50 C LARGE ON KODE=1
51 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
52 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
53 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE
54 C ACCURACY
55 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
56 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
57 C CANCE BY ARGUMENT REDUCTION
58 C IERR=5, ERROR - NO COMPUTATION,
59 C ALGORITHM TERMINATION CONDITION NOT MET
60 C
61 C***LONG DESCRIPTION
62 C
63 C THE COMPUTATION IS CARRIED OUT BY THE POWER SERIES FOR
64 C SMALL CABS(Z), THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z),
65 C THE MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN AND A
66 C NEUMANN SERIES FOR IMTERMEDIATE MAGNITUDES, AND THE
67 C UNIFORM ASYMPTOTIC EXPANSIONS FOR I(FNU,Z) AND J(FNU,Z)
68 C FOR LARGE ORDERS. BACKWARD RECURRENCE IS USED TO GENERATE
69 C SEQUENCES OR REDUCE ORDERS WHEN NECESSARY.
70 C
71 C THE CALCULATIONS ABOVE ARE DONE IN THE RIGHT HALF PLANE AND
72 C CONTINUED INTO THE LEFT HALF PLANE BY THE FORMULA
73 C
74 C I(FNU,Z*EXP(M*PI)) = EXP(M*PI*FNU)*I(FNU,Z) REAL(Z).GT.0.0
75 C M = +I OR -I, I**2=-1
76 C
77 C FOR NEGATIVE ORDERS,THE FORMULA
78 C
79 C I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z)
80 C
81 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE
82 C THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE
83 C INTEGER,THE MAGNITUDE OF I(-FNU,Z)=I(FNU,Z) IS A LARGE
84 C NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER,
85 C K(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF
86 C TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY
87 C UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN
88 C OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE,
89 C LARGE MEANS FNU.GT.CABS(Z).
90 C
91 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
92 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
93 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
94 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
95 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
96 C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO
97 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
98 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
99 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
100 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
101 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
102 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
103 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
104 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
105 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
106 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
107 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
108 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
109 C
110 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
111 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
112 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
113 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
114 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
115 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
116 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
117 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
118 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
119 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
120 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
121 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
122 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
123 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
124 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
125 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
126 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
127 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
128 C OR -PI/2+P.
129 C
130 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
131 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
132 C COMMERCE, 1955.
133 C
134 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
135 C BY D. E. AMOS, SAND83-0083, MAY, 1983.
136 C
137 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
138 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
139 C
140 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
141 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
142 C 1018, MAY, 1985
143 C
144 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
145 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
146 C MATH. SOFTWARE, 1986
147 C
148 C***ROUTINES CALLED CBINU,I1MACH,R1MACH
149 C***END PROLOGUE CBESI
150  COMPLEX cone, csgn, cy, z, zn
151  REAL aa, alim, arg, dig, elim, fnu, fnul, pi, rl, r1m5, s1, s2,
152  * tol, xx, yy, r1mach, az, fn, bb, ascle, rtol, atol
153  INTEGER i, ierr, inu, k, kode, k1, k2, n, nn, nz, i1mach
154  dimension cy(n)
155  DATA pi /3.14159265358979324e0/
156  DATA cone / (1.0e0,0.0e0) /
157 C
158 C***FIRST EXECUTABLE STATEMENT CBESI
159  ierr = 0
160  nz=0
161  IF (fnu.LT.0.0e0) ierr=1
162  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
163  IF (n.LT.1) ierr=1
164  IF (ierr.NE.0) RETURN
165  xx = REAL(z)
166  yy = aimag(z)
167 C-----------------------------------------------------------------------
168 C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
169 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
170 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
171 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
172 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
173 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
174 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
175 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
176 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU.
177 C-----------------------------------------------------------------------
178  tol = amax1(r1mach(4),1.0e-18)
179  k1 = i1mach(12)
180  k2 = i1mach(13)
181  r1m5 = r1mach(5)
182  k = min0(iabs(k1),iabs(k2))
183  elim = 2.303e0*(float(k)*r1m5-3.0e0)
184  k1 = i1mach(11) - 1
185  aa = r1m5*float(k1)
186  dig = amin1(aa,18.0e0)
187  aa = aa*2.303e0
188  alim = elim + amax1(-aa,-41.45e0)
189  rl = 1.2e0*dig + 3.0e0
190  fnul = 10.0e0 + 6.0e0*(dig-3.0e0)
191  az = cabs(z)
192 C-----------------------------------------------------------------------
193 C TEST FOR RANGE
194 C-----------------------------------------------------------------------
195  aa = 0.5e0/tol
196  bb=float(i1mach(9))*0.5e0
197  aa=amin1(aa,bb)
198  IF(az.GT.aa) go to 140
199  fn=fnu+float(n-1)
200  IF(fn.GT.aa) go to 140
201  aa=sqrt(aa)
202  IF(az.GT.aa) ierr=3
203  IF(fn.GT.aa) ierr=3
204  zn = z
205  csgn = cone
206  IF (xx.GE.0.0e0) go to 40
207  zn = -z
208 C-----------------------------------------------------------------------
209 C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
210 C WHEN FNU IS LARGE
211 C-----------------------------------------------------------------------
212  inu = int(fnu)
213  arg = (fnu-float(inu))*pi
214  IF (yy.LT.0.0e0) arg = -arg
215  s1 = cos(arg)
216  s2 = sin(arg)
217  csgn = cmplx(s1,s2)
218  IF (mod(inu,2).EQ.1) csgn = -csgn
219  40 CONTINUE
220  CALL cbinu(zn, fnu, kode, n, cy, nz, rl, fnul, tol, elim, alim)
221  IF (nz.LT.0) go to 120
222  IF (xx.GE.0.0e0) RETURN
223 C-----------------------------------------------------------------------
224 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE
225 C-----------------------------------------------------------------------
226  nn = n - nz
227  IF (nn.EQ.0) RETURN
228  rtol = 1.0e0/tol
229  ascle = r1mach(1)*rtol*1.0e+3
230  DO 50 i=1,nn
231 C CY(I) = CY(I)*CSGN
232  zn=cy(i)
233  aa=REAL(zn)
234  bb=aimag(zn)
235  atol=1.0e0
236  IF (amax1(abs(aa),abs(bb)).GT.ascle) go to 55
237  zn = zn*cmplx(rtol,0.0e0)
238  atol = tol
239  55 CONTINUE
240  zn = zn*csgn
241  cy(i) = zn*cmplx(atol,0.0e0)
242  csgn = -csgn
243  50 CONTINUE
244  RETURN
245  120 CONTINUE
246  IF(nz.EQ.(-2)) go to 130
247  nz = 0
248  ierr=2
249  RETURN
250  130 CONTINUE
251  nz=0
252  ierr=5
253  RETURN
254  140 CONTINUE
255  nz=0
256  ierr=4
257  RETURN
258  END