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