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