CBYS
Evaluates a sequence of Bessel functions of the second kind with real order and complex arguments.
Required Arguments
XNU — Real argument which is the lowest order desired. (Input)
XNU must be greater than –1/2.
Z — Complex argument for which the sequence of Bessel functions is to be evaluated. (Input)
N — Number of elements in the sequence. (Input)
CBS — Vector of length N containing the values of the function through the series. (Output)
CBS(I) contains the value of the Bessel function of order XNU + I ‑ 1 at Z for I = 1 to N.
FORTRAN 90 Interface
Generic: CALL CBYS (XNU, Z, N, CBS)
Specific: The specific interface names are S_CBYS and D_CBYS.
FORTRAN 77 Interface
Single: CALL CBYS (XNU, Z, N, CBS)
Double: The double precision name is DCBYS.
Description
The Bessel function Yν(z) is defined to be
This code is based on the code BESSEC of Barnett (1981) and Thompson and Barnett (1987).
This code computes
Yν(
z) from the modified Bessel functions
Iν(
z) and
Kν(
z),
CBIS and
CBKS, using the following relation:
CBYS implements the Yousif and Melka (Y&M) algorithm (
Yousif and Melka(2003)) for approximating
Yν (
z) with a Taylor series expansion when
x ~ 0 or
y ~ 0, where complex argument
z =
x +
iy and “
x ~ 0” == “
∣x∣ <
amach(4)”. To be consistent with the existing
CBYS argument definitions, the original Y&M algorithm, which was limited to integral order and to (
x ~ 0 and
y ≥ 0) or (
y ~ 0 and
x ≥ 0), has been generalized to also work for integral and real order
ν >
-1, and for (
x ~ 0 and
y < 0) and (
y ~ 0 and
x < 0).
To deal with the Bessel function discontinuity occurring at the negative x axis, the following procedure is used for calculating the Y&M approximation of Yν (z) with argument z = x + iy when ((x ~ 0 and y < 0) or (y ~ 0 and x < 0)):
1. Calculate the Y&M approximation of Yν (-z).
2. If (y > 0), use forward rotation, otherwise use backward rotation, to calculate the Bessel function Yν (z), where the “forward” and “backward” rotation transformations are defined as:
forward: Yν (z) = i‑2ν Yν (-z) + 2i cos(νπ)Jν (-z)
backward: Yν (z) = i2ν Yν (-z) - 2i cos(νπ)Jν (-z)
These definitions are based on
Abromowitz and Stegun (1972), eq. 9.1.36:
Yν (
zemπ i ) =
e‑mνπ i Yν (
z) + 2
i sin(
mνπ) cot(
νπ)
Jν (
z), where
m = 1 represents forward transformation and
m =
-1 represents backward transformation. These specified rotations insure that the continuous rotation transformation
Yν (
-z)
→ Yν (
z) does not cross the negative x axis, so no discontinuity is encountered.
Comments
1. Workspace may be explicitly provided, if desired, by use of C2YS/DC2Y. The reference is:
CALL C2YS (XNU, Z, N, CBS, FK)
The additional argument is:
FK — complex work vector of length N.
2. Informational errors
Type | Code | Description |
|---|
3 | 1 | One of the continued fractions failed. |
4 | 2 | Only the first several entries in CBS are valid. |
Example
In this example, Y0.3+k−1(1.2 + 0.5i), k = 1, …, 4 is computed and printed.
USE CBYS_INT
USE UMACH_INT
IMPLICIT NONE
! Declare variables
INTEGER N
PARAMETER (N=4)
!
INTEGER K, NOUT
REAL XNU
COMPLEX CBS(N), Z
! Compute
XNU = 0.3
Z = (1.2, 0.5)
CALL CBYS (XNU, Z, N, CBS)
! Print the results
CALL UMACH (2, NOUT)
DO 10 K=1, N
WRITE (NOUT,99999) XNU+K-1, Z, CBS(K)
10 CONTINUE
99999 FORMAT (' Y sub ', F6.3, ' ((', F6.3, ',', F6.3, &
')) = (', F9.3, ',', F9.3, ')')
END
Output
Y sub 0.300 (( 1.200, 0.500)) = ( -0.013, 0.380)
Y sub 1.300 (( 1.200, 0.500)) = ( -0.716, 0.338)
Y sub 2.300 (( 1.200, 0.500)) = ( -1.048, 0.795)
Y sub 3.300 (( 1.200, 0.500)) = ( -1.625, 3.684)
