GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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  35 CONTINUE
205  zn = z
206  csgn = cone
207  IF (xx.GE.0.0e0) GO TO 40
208  zn = -z
209 C-----------------------------------------------------------------------
210 C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
211 C WHEN FNU IS LARGE
212 C-----------------------------------------------------------------------
213  inu = int(fnu)
214  arg = (fnu-float(inu))*pi
215  IF (yy.LT.0.0e0) arg = -arg
216  s1 = cos(arg)
217  s2 = sin(arg)
218  csgn = cmplx(s1,s2)
219  IF (mod(inu,2).EQ.1) csgn = -csgn
220  40 CONTINUE
221  CALL cbinu(zn, fnu, kode, n, cy, nz, rl, fnul, tol, elim, alim)
222  IF (nz.LT.0) GO TO 120
223  IF (xx.GE.0.0e0) RETURN
224 C-----------------------------------------------------------------------
225 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE
226 C-----------------------------------------------------------------------
227  nn = n - nz
228  IF (nn.EQ.0) RETURN
229  rtol = 1.0e0/tol
230  ascle = r1mach(1)*rtol*1.0e+3
231  DO 50 i=1,nn
232 C CY(I) = CY(I)*CSGN
233  zn=cy(i)
234  aa=REAL(ZN)
235  bb=aimag(zn)
236  atol=1.0e0
237  IF (amax1(abs(aa),abs(bb)).GT.ascle) GO TO 55
238  zn = zn*cmplx(rtol,0.0e0)
239  atol = tol
240  55 CONTINUE
241  zn = zn*csgn
242  cy(i) = zn*cmplx(atol,0.0e0)
243  csgn = -csgn
244  50 CONTINUE
245  RETURN
246  120 CONTINUE
247  IF(nz.EQ.(-2)) GO TO 130
248  nz = 0
249  ierr=2
250  RETURN
251  130 CONTINUE
252  nz=0
253  ierr=5
254  RETURN
255  140 CONTINUE
256  ierr=4
257  GO TO 35
258  END
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:860
static T abs(T x)
Definition: pr-output.cc:1696
double complex cmplx
Definition: Faddeeva.cc:217
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)