GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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  35 CONTINUE
206  ufl = r1mach(1)*1.0e+3
207  IF (az.LT.ufl) GO TO 180
208  IF (fnu.GT.fnul) GO TO 80
209  IF (fn.LE.1.0e0) GO TO 60
210  IF (fn.GT.2.0e0) GO TO 50
211  IF (az.GT.tol) GO TO 60
212  arg = 0.5e0*az
213  aln = -fn*alog(arg)
214  IF (aln.GT.elim) GO TO 180
215  GO TO 60
216  50 CONTINUE
217  CALL cuoik(z, fnu, kode, 2, nn, cy, nuf, tol, elim, alim)
218  IF (nuf.LT.0) GO TO 180
219  nz = nz + nuf
220  nn = nn - nuf
221 C-----------------------------------------------------------------------
222 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
223 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
224 C-----------------------------------------------------------------------
225  IF (nn.EQ.0) GO TO 100
226  60 CONTINUE
227  IF (xx.LT.0.0e0) GO TO 70
228 C-----------------------------------------------------------------------
229 C RIGHT HALF PLANE COMPUTATION, REAL(Z).GE.0.
230 C-----------------------------------------------------------------------
231  CALL cbknu(z, fnu, kode, nn, cy, nw, tol, elim, alim)
232  IF (nw.LT.0) GO TO 200
233  nz=nw
234  RETURN
235 C-----------------------------------------------------------------------
236 C LEFT HALF PLANE COMPUTATION
237 C PI/2.LT.ARG(Z).LE.PI AND -PI.LT.ARG(Z).LT.-PI/2.
238 C-----------------------------------------------------------------------
239  70 CONTINUE
240  IF (nz.NE.0) GO TO 180
241  mr = 1
242  IF (yy.LT.0.0e0) mr = -1
243  CALL cacon(z, fnu, kode, mr, nn, cy, nw, rl, fnul, tol, elim,
244  * alim)
245  IF (nw.LT.0) GO TO 200
246  nz=nw
247  RETURN
248 C-----------------------------------------------------------------------
249 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
250 C-----------------------------------------------------------------------
251  80 CONTINUE
252  mr = 0
253  IF (xx.GE.0.0e0) GO TO 90
254  mr = 1
255  IF (yy.LT.0.0e0) mr = -1
256  90 CONTINUE
257  CALL cbunk(z, fnu, kode, mr, nn, cy, nw, tol, elim, alim)
258  IF (nw.LT.0) GO TO 200
259  nz = nz + nw
260  RETURN
261  100 CONTINUE
262  IF (xx.LT.0.0e0) GO TO 180
263  RETURN
264  180 CONTINUE
265  nz = 0
266  ierr=2
267  RETURN
268  200 CONTINUE
269  IF(nw.EQ.(-1)) GO TO 180
270  nz=0
271  ierr=5
272  RETURN
273  210 CONTINUE
274  ierr=4
275  GO TO 35
276  END
subroutine cacon(Z, FNU, KODE, MR, N, Y, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: cacon.f:3
subroutine cbesk(Z, FNU, KODE, N, CY, NZ, IERR)
Definition: cbesk.f:2
subroutine cbknu(Z, FNU, KODE, N, Y, NZ, TOL, ELIM, ALIM)
Definition: cbknu.f:2
subroutine cbunk(Z, FNU, KODE, MR, N, Y, NZ, TOL, ELIM, ALIM)
Definition: cbunk.f:2
subroutine cuoik(Z, FNU, KODE, IKFLG, N, Y, NUF, TOL, ELIM, ALIM)
Definition: cuoik.f:2
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137