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
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  ufl = d1mach(1)*1.0d+3
209  IF (az.LT.ufl) go to 180
210  IF (fnu.GT.fnul) go to 80
211  IF (fn.LE.1.0d0) go to 60
212  IF (fn.GT.2.0d0) go to 50
213  IF (az.GT.tol) go to 60
214  arg = 0.5d0*az
215  aln = -fn*dlog(arg)
216  IF (aln.GT.elim) go to 180
217  go to 60
218  50 CONTINUE
219  CALL zuoik(zr, zi, fnu, kode, 2, nn, cyr, cyi, nuf, tol, elim,
220  * alim)
221  IF (nuf.LT.0) go to 180
222  nz = nz + nuf
223  nn = nn - nuf
224 C-----------------------------------------------------------------------
225 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
226 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
227 C-----------------------------------------------------------------------
228  IF (nn.EQ.0) go to 100
229  60 CONTINUE
230  IF (zr.LT.0.0d0) go to 70
231 C-----------------------------------------------------------------------
232 C RIGHT HALF PLANE COMPUTATION, REAL(Z).GE.0.
233 C-----------------------------------------------------------------------
234  CALL zbknu(zr, zi, fnu, kode, nn, cyr, cyi, nw, tol, elim, alim)
235  IF (nw.LT.0) go to 200
236  nz=nw
237  RETURN
238 C-----------------------------------------------------------------------
239 C LEFT HALF PLANE COMPUTATION
240 C PI/2.LT.ARG(Z).LE.PI AND -PI.LT.ARG(Z).LT.-PI/2.
241 C-----------------------------------------------------------------------
242  70 CONTINUE
243  IF (nz.NE.0) go to 180
244  mr = 1
245  IF (zi.LT.0.0d0) mr = -1
246  CALL zacon(zr, zi, fnu, kode, mr, nn, cyr, cyi, nw, rl, fnul,
247  * tol, elim, alim)
248  IF (nw.LT.0) go to 200
249  nz=nw
250  RETURN
251 C-----------------------------------------------------------------------
252 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
253 C-----------------------------------------------------------------------
254  80 CONTINUE
255  mr = 0
256  IF (zr.GE.0.0d0) go to 90
257  mr = 1
258  IF (zi.LT.0.0d0) mr = -1
259  90 CONTINUE
260  CALL zbunk(zr, zi, fnu, kode, mr, nn, cyr, cyi, nw, tol, elim,
261  * alim)
262  IF (nw.LT.0) go to 200
263  nz = nz + nw
264  RETURN
265  100 CONTINUE
266  IF (zr.LT.0.0d0) go to 180
267  RETURN
268  180 CONTINUE
269  nz = 0
270  ierr=2
271  RETURN
272  200 CONTINUE
273  IF(nw.EQ.(-1)) go to 180
274  nz=0
275  ierr=5
276  RETURN
277  260 CONTINUE
278  nz=0
279  ierr=4
280  RETURN
281  END
std::string dimension(void) const
subroutine zacon(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: zacon.f:1
subroutine zuoik(ZR, ZI, FNU, KODE, IKFLG, N, YR, YI, NUF, TOL, ELIM, ALIM)
Definition: zuoik.f:1
double precision function d1mach(i)
Definition: d1mach.f:1
F77_RET_T const double const double double * d
subroutine zbunk(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, ALIM)
Definition: zbunk.f:1
subroutine zbknu(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM)
Definition: zbknu.f:1
subroutine zbesk(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)
Definition: zbesk.f:1