\input multip.cmm
\advance\hoffset by 0.25in
\advance\voffset by 0.60in
\def\frac#1/#2{\leavevmode\kern.1em\raise.5ex\hbox{\the\scriptfont0 #1}
\kern-.1em/\kern-.15em\lower.25ex\hbox{\the\scriptfont0 #2}}
\head{The Calculation of Spherical Bessel Functions and Coulomb Functions}
\author{A.R.\ts Barnett}
\address{Physics Department, University of Manchester, Manchester, M13 9PL,
England and ~~~Physics Department, University of Auckland,
Auckland, New Zealand}
\titlea{Abstract}{\null}
An account is given of the Steed algorithm for calculating Coulomb functions
and, as a special case, both spherical Bessel and Riccati--Bessel functions.
These functions are needed for boundary-condition matching in scattering
problems in Atomic and Nuclear physics. Central to the technique is the
evaluation of continued fractions and for this calculation Lentz's forward
method (modified by Thompson) is recommended. The FORTRAN~77 programs
SBESJY, RICBES and COUL90 are presented and described, and some test
cases are given. In each program the algorithm returns both the regular and
irregular functions as well as their derivatives. The programs are written
for real arguments and real orders; references guide the reader to more
general codes.
\vfill
\eject
\titlea{1.}{Introduction}
Coulomb wave functions arise in many problems of physical interest when
charged particles scatter from each other. Such scattering is characterised
by a relative angular momentum $L \hbar$ (with $L$ a non-negative integer)
and by a Sommerfeld parameter $\eta = Z\alpha/\beta$ which gives the strength
of the Coulomb interaction. The product of the particle charges is $Ze^2$,
the fine-structure constant $\alpha = e^2/\hbar c(4\pi\epsilon_0)$, and
$\beta c$ is the relative velocity of the particles. With charges of opposite
sign, as is frequently the case in Atomic physics, then $\eta$ is negative,
while in Nuclear physics problems generally the charges have the same sign
and $\eta$
is positive. Positron scattering from a nucleus will also have $\eta > 0$.
For the scattering of a neutral particle ({\it e.g.} a neutron) or by a
neutral target then the parameter $\eta$ is zero, resulting in
Riccati--Bessel functions in which only the angular momentum effects appear.
These functions differ only in a minor way from spherical Bessel functions
and both can properly be regarded as special cases of the Coulomb functions.
This review is mainly concerned with the Atomic physics area and thus
primarily with the cases where $\eta$ is not positive. In atomic units
the relative energy $E$, in Rydberg units, is given by $E = E(eV)/13.605~eV$.
Then $E = k^2$ and $\eta = -Z/k$. Typical energies are between $0.5 - 200
~eV$, with corresponding $k$-values from $(0.2 - 3.8)a_0^{-1}$, and
$\eta$-values between $-5.2$ and $-0.26$. With matching
radii of $r=(10 - 200)a_0$ then the dimensionless variable $x = kr$ lies in
the range $1 - 1000$. The angular momentum quantum number $L$ will
vary from 0 to perhaps 50, while Bessel functions of orders up to 150 may be
required. As an example, the parameters for the calculations in Chapter 1 of
this book, with energy 54.5 $eV$ are
$k =$ 2.00, $\eta = -0.50$ (0.50) for electron (positron) scattering off
hydrogen ions, and $\eta = 0$ for either particle scattering from neutral
hydrogen.
The second-order differential equation satisfied by the Coulomb functions
is,
$$w^{\prime\prime} (x) \,+\, [1 - 2\eta/x - L(L+1)/x^2 ]\,w(x) = 0
\eqno{(1)}$$
(where the primes indicate differentiation with respect to $x$) and it has
two linearly independent solutions. They are chosen to be the regular
solution $F_L(\eta,x)$ which is zero at $x=0$, and the irregular solution
$G_L(\eta,x)$ which is infinite at $x=0$. See also refs$^{1,2,3,9,10}$ for
alternative forms of (1).
The `turning point' for the $L^{th}$ partial wave occurs at a value of $x$
$$x_L = \eta \,+\, (\eta^2 + L^2 + L)^{1/2}\eqno{(2)}$$
which is where the bracket [...] is zero in (1); it is also the first
point of inflexion for the functions. For negative $\eta$, say $\eta=-h$,
the $L=0$ turning point is at the origin, $x_0 = 0$, and
for other $L$-values $x_L$ is always less than $L$; it lies between $L-h$
(for small $h$) and $\frac1/2L(L/h)$ (for small $L$).
It will transpire that the computational task is easier when $x>x_L$, so
that functions with negative $\eta$ are obtained more straightforwardly
than when $\eta$ is positive.
For values of $x$ which are greater than $x_L$ the functions take on an
oscillatory character, although the `period' slowly changes.
Examples of the functions are shown in Fig.1. For $\eta>0$
the regular function magnitude is greater than unity, and it slowly
decreases towards unity as $x$ grows larger. When $\eta<0$ the magnitude of
$F_L(\eta,x)$ is less than unity and it increases steadily for larger $x$.
In the asymptotic region as $x
\rightarrow \infty$ they become circular with
$F \rightarrow \sin \theta_L$ and $G
\rightarrow \cos \theta_L$ , where $\theta_L$ is the asymptotic phase given
by $$\theta_L = x \,-\, \eta ln(2x) \,-\, \frac1/2L\pi \,+\, {\rm arg}
\Gamma(L+1+i\eta)\eqno{(3)}$$ and in the case when $\eta = 0$ this phase
becomes linear in $x$: $$\theta_L = x \,-\,\frac1/2L\pi$$ An additional phase
of $\frac1/2\pi$ will be required for spherical Bessel functions because of
a different definition which is given in eqn(8).
This article deals with solutions to the Coulomb scattering problem in which
the particles have a positive relative energy. Curtis$^{1}$ published a
detailed discussion for $\eta < 0$ in 1964 for $L = 0,1,2$ and for functions
closely related to $F$ and $G$, for use in electron scattering calculations.
More recent work for electron scattering is that of Seaton$^{2}$ and of Bell
and Scott$^{3}$; Seaton dealing with positive, negative and zero energies,
while Bell and Scott treated negative energies. The work of Curtis for
positive energy was verified numerically by Barnett$^{4}$ in 1974. Bardin
{\it et al.}$^{5}$ published a comprehensive suite of programs in 1972
containing codes for all real $\eta$ and $x$. All of these methods involve
some appropriate series expansion, either in powers of $x$, or
asymptotically in powers of $x^{-1}$. Both $F$ and $G$ are calculated
separately, as are their $x$-derivatives. A new approach to the calculation
of Coulomb functions was developed by Steed in the late 1960s and was
published by Barnett {\it et al.}$^{6}$ in 1974. The functions and their
derivatives are calculated {\it together} in an interdependent way,
and two continued fractions are used (\S6). Virtually all previous work,
especially that of Bardin {\it et al.,}$^{5}$ was verified and in some cases
extended in this paper. A detailed description of this algorithm and a
careful comparison with other approaches was given by Barnett$^{7}$ in 1982.
The original program$^{6}$ RCWFN was superseded by a more comprehensive
version,$^{4}$ called COULFG, which includes the calculation of both Bessel
and spherical Bessel functions (of the first and second kind {\it i.e.}
$j_n(x), y_n(x)$) as well as the Coulomb
functions. The version presented here in \S7, COUL90, has been further
improved in one important respect which, however, does not affect the
algorithm, the range of the calculations or the results already given, both
in ref.$^{4}$ and the earlier work. Comments on the methods of COULFG have
been made by Nesbet$^{8}$ with suggestions for improvements for both
Coulomb--function calculations when $x 0$, can be
obtained by using the recent program COULN of Noble and Thompson$^{9}$
which is well suited to electron scattering.
The algorithm of COULFG is not restricted to integer L-values, as
refs$^{4,10}$ demonstrate, and it can hence be used for scattering
solutions to relativistic problems for which the Klein--Gordon equation or
the Dirac equation are appropriate. Here an equivalent non-integer $L$ is
required; for small $L$-values this can be imaginary. Extensions of the
concepts underlying the calculations to {\it complex arguments} are
presented in the program COULCC, which was published in 1985 by Thompson
and Barnett.$^{11}$ The range of each of the three variables has been
extended into the complex plane. The program has been recently been
incorporated into the IMSL SFUN library. A description has been given in
ref.$^{12}$ of the various approaches required in the different parameter
regions, with references to earlier and more restricted work. The code
COULCC also computes Bessel functions with complex arguments and order;
however, for real order, BESSCC for modified Bessel functions is a more
efficient code,$^{11}$ which incorporates similar principles.
\titlea{2.}{Spherical Bessel Functions}
Spherical Bessel functions, and their close relatives the Riccati--Bessel
functions, are required frequently in Atomic physics calculations and in
many other branches of Physics, for example in the calculations of Mie
scattering in Optics.
In Chapter 1 of this book they emerge as
the asymptotic wavefunctions for the problem of charged particles
scattering from neutral atoms or molecules. Similarly, they are required
for the description of the scattering of neutrons from nuclei. Despite
this they rather tend to be the poor relations of the numerical analysis
world. They are covered in the encylopaedic {\it Handbook} of Abramowitz and
Stegun$^{13}$ but they are not
mentioned, for example, in the {\it Numerical Recipes},$^{14}$ chapter 6 on
Special Functions, where Bessel functions ({\it i.e.} with {\it cylindrical}
symmetry) are treated, and few, if any, large-scale libraries (such as NAG,
IMSL {\it etc.}) have suitable subroutines. The reason is
simply that all orders are expressible as sums of polynomials in $x^{-1}$
multiplied by $\sin x$ and $(-\cos x)$, and recurrence
relations connect consecutive orders. The computational task appears
trivial. Nevertheless there are difficulties, shared with the evaluation of
more general Bessel functions and Coulomb functions, which stem directly from
the relevant differential equation. A useful discussion, intended for the
physicist exploring numerical analysis, of some of these
difficulties in the case of cylindrical Bessel functions, is that of
Koonin$^{15}$(\S4.1), which also contains the details of an
fixed-accuracy program to evalute them. Similarly, {\it Numerical Recipes} by
Press {\it et al.}$^{14}$ contains an excellent coverage of the
counter-intuitive aspects
of the evaluation of special functions, of difference and of
differential equations, as well as numerous programs.
Very recently however, Press and Teukolsky, two of the authors of
{\it Numerical Recipes}, have published in their regular column$^{34}$ programs
for $J_{\nu}(x),~Y_{\nu}(x)$ for real cylindrical Bessel functions and for the
modified $I_{\nu}(x),~K_{\nu}(x)$ Bessels of real argument $x$ and real order
$\nu$. These programs
directly adopt Steed's method and the extensions and the techniques
given in the complex-argument, real-order code BESSCC$^{11}$ of Thompson
and Barnett. The Press-Teukolsky codes are thus more suitable for real arguments
and are faster
since they are less general. As examples of their use, calling programs are
provided to find Airy functions and spherical Bessel functions. This route
makes the spherical Bessels conceptually a special case of the usual Bessel
function; it is also evident$^{23}$ that {\it both} are special cases of the
Coulomb functions.
No test cases are given, nor are ranges specified on the parameters, but
in general the same conclusions as those presented in this article will apply.
Two programs have been
recently published for the specific calculation of spherical Bessel functions,
by Gillman and Fiebig$^{16}$ and by Lentz.$^{26}$ They are discussed at the end
of \S5.
For {\bf spherical Bessel functions}, with solutions $j_n(x)$ and
$y_n(x)$, the differential equation is:
$$w^{\prime\prime} \,+\, 2w^\prime/x \,+\, [1 - n(n+1)/x^2]\,w = 0
~\eqno{(4)}$$
In this equation (see ref.$^{13}$ eqn 10.1.1) $x$ is a real variable while
$n$ is
an integer which is positive, negative or zero and which will be identified
with the non-negative angular momentum. Many of the properties of $j_n$ and
$y_n$ follow directly from their identification as Coulomb functions with
$\eta = 0$, as is done below through eqn(11). The results are $$j_n(x) =
x^{-1}F_n(0,x)~~~~~~~~~~y_n(x) = -x^{-1}G_n(0,x)\eqno{(5)}$$
>From this and the properties of $F$ and $G$ we
immediately deduce that $j$ and $y$ will display an oscillatory nature for
$x$-values larger than $\sqrt{(n^2 + n)}$, their magnitude will decrease
as
$x^{-1}$ for these larger $x$-values, and that $y$ will diverge towards
$-\infty$ as $x$ approaches zero. The first two orders of the
{\bf spherical Bessel functions} are:
$$\eqalign{j_0(x) &= x^{-1}\sin x~~~~~~~~~~~~~~~~~~~~~~~~~~~
y_0(x) = ~-~ x^{-1}\cos x \cr
j_1(x) &= x^{-2}\sin x ~-~ x^{-1}\cos x~~~~~~~~~~
y_1(x) = -x^{-2}\cos x ~-~ x^{-1}\sin x \cr}\eqno{(6)}$$
and recurrence relations (\S3) link sets of three consecutive orders. The
derivatives $j'_0(x),~y'_0(x)$ follow from (6).
The two solutions $j_n$ and $y_n$ are called the {\it regular solution} and
the {\it irregular solution} respectively. This describes their behaviour at
the origin of $x$ where the irregular solution diverges to $-\infty$ as
$-(2n-1)!!x^{-(n+1)}$ and the regular one goes to zero as $x^n/(2n+1)!!$,
(and with $j_0(0) = 1)$. The irregular solution, $y_n(x)$, is called the
spherical Neumann function and occasionally given the symbol $n_n(x)$. The
following linear combinations are the spherical Hankel functions,
$$\eqalign{h^{(1)}_n(x) &= j_n(x) \,+\, iy_n(x) \cr
h^{(2)}_n(x) &= j_n(x) \,-\, iy_n(x) \cr}\eqno{(7)}$$
As $x$ becomes large the equations of
ref.$^{13}$ (9.2.1, 9.2.17 and those following 10.1.26) show that
$$\eqalign{j_n(x) \,\rightarrow\, x^{-1}\cos(x - \frac1/2n\pi -
\frac1/2\pi) &~=~~~x^{-1}\sin(\theta_n - \frac1/2\pi),~{\rm and} \cr
y_n(x) \,\rightarrow\, x^{-1}\sin(x -\frac1/2n\pi- \frac1/2\pi)
&~=~ -x^{-1}\cos(\theta_n - \frac1/2\pi) \cr}\eqno{(8)}$$
This difference in phase definition (the additional $\frac1/2\pi$ which
interchanges the role of sine and cosine in (8))
is carried over to the equivalents of the Hankel functions for Coulomb
scattering. For this Coulomb case the expression for the outgoing wave is
chosen to be
$$H^+_L (\eta,x)\,=\, G_L (\eta,x)\,+\,iF_L(\eta,x) \,\rightarrow\, \exp(i\theta_L
)~{\rm as}~~~x \rightarrow \infty\eqno{(9)}$$
and similarly for $H^-$ the incoming wave. For $\eta=0$ the relations are
$$h^{(1)}(x) = -(i/x)H^+ (0,x)~~~~~~~~~~ h^{(2)}(x) =
(i/x)H^-(0,x)\eqno{(10)}$$
with the $\frac1/2\pi$ phase difference remaining. The $H^{\pm}$ scattering
functions must not be confused with the Hankel functions for cylindrical
Bessels, $H^{(1)}_n$ and $H^{(2)}_n$ (see ref.$^{13}$ eqn
9.1.3, 9.1.4, 9.1.6), which are analogous to (7).
To obtain the {\bf Riccati--Bessel functions} we transform (4), by removing
the first-derivative term, into
$$w^{\prime\prime} \,+\, [1 - n(n+1)/x^2]\,w = 0~ \eqno{(11)}$$
whose solutions (ref.$^{13}~\S10.3$) are $xj_n(x)$ and $xy_n(x)$. No
special symbol has been given to them in the mathematical literature
although, for complex argument $z$, a consistent notation is
used$^{17,18,19}$
for Mie scattering:
$$\eqalign{ \psi_n(z) &~=~~~zj_n(z) \cr
\chi_n(z) &~=~-zy_n(z) \cr
\zeta_n(z) &~=~\psi_n(z) ~+~i \chi_n(z)\cr} \eqno{(12)}$$
The Riccati--Bessel properties are briefly treated in the
{\it Handbook}$^{13}$[ \S10.3]; here, however, $\chi_n~=~+zy_n$.
For orders $n~=~0,1$ we then have, adopting ref.$^{13}$
$$\eqalign{xj_n(x) &= \sin x~~~~~~~~~~~~~~~~~~~~~~~~~~~~ xy_0(x)= -\cos x
\cr
xj_1(x) &= x^{-1}\sin x ~-~\cos x~~~~~~~~~~~xy_1(x) = -x^{-1}\cos x ~-~
\sin x \cr}\eqno{(13)}$$
On comparison of (11) with (1) it is evident that the Riccati--Bessel
functions are the same as the Coulomb functions with $\eta = 0$, while (8)
and (9) show that the sign of
the irregular function $chi_n$ is reversed relative to $G_L(0,x)$.
The Riccati--Bessel functions thus behave
like $\sin\theta_n$ and $(-\cos\theta_n)$ as $x$ becomes large, where
the asymptotic phase, $\theta_n$, is given by (3). Alternatively this can be
written $\cos(\theta_n -\frac1/2\pi)$ , $\sin(\theta_n - \frac1/2\pi)$ as in
(8). The subroutine RICBES retains the {\it Handbook} definitions in which
the Riccati--Bessel is the argument $x$ times the corresponding spherical
Bessel function.
\titlea{3.}{Recurrence Relations for Spherical Bessel Functions}
Each of the four spherical Bessel functions $j_n,~y_n,~h^{(1)}_n,~h^{(2)}_n$
obeys recurrence relations [ref.$^{13}$ 10.1.19 -- 10.1.22] which connect
the functions of order $(n-1)$, $n$ and $(n+1)$. Using $g_n$ to
represent any of these four functions the relations are:
$$g_{n-1} \,-\, {{2n + 1}\over{x}}\,g_n \,+\, g_{n+1} = 0\eqno{(14)}$$
and
$$ng_{n-1} \,-\, {{2n + 1}\over{x}}\,g^{\prime}_n \,-\, (n+1)g_{n+1} =
0.\eqno{(15)}$$
These can be rewritten in a form which is suitable for {\bf downward
recurrence}, connecting two successive orders and a derivative,
$$g_{n-1} = S_{n+1}g_n \,+\, g^{\prime}_n\eqno{(16)}$$
$$g^{\prime}_{n-1} = S_{n-1}g_{n-1} \,-\, g_n\eqno{(17)}$$
in which $S_n = n/x$. The equivalent expressions for {\bf upward
recurrence} are equations (17) and (16) rearranged:
$$g_{n+1} = S_ng_n \,-\, g^{\prime}_n \eqno{(18)}$$
and
$$g^{\prime}_{n+1} = g_n \,-\, S_{n+2}g_{n+1} \eqno{(19)}$$
The recurrence relations (14) -- (17) are an alternative way of
expressing the differential equation (4) as difference equations.
It is well known ({\it e.g.} ref.$^{14}$ \S5.4, ref.$^{15}$ \S4.1, and
ref.$^{8}$ \S3) that
a recurrence relation is numerically unstable in the direction
in which the function is decreasing. Successive values are computed as
small differences between nearly equal terms and all accuracy is soon lost.
This occurs for
the regular function $j_L(x)$ once $L>x$ for any fixed value of $x$:
$j_L(x)$ {\it decreases monotonically} as a function of $L$ ,
so that upward recurrence in $n$ of $j_n(x)$ is unstable. Conversely,
since the irregular function $y_L(x)$ {\it increases }
as $L$ increases, once $L>x$, upward recurrence is stable.
Thus for $L>x$ we must use {\it downward recurrence in} $n$
to calculate the values of $j_n(x)$ ($n = L,L-1,...,3,2,1,0)$; similarly {\it
upward recurrence} must be used for $y_n(x)$. In the region where $LFrom (18), which is satisfied by $j_n(x)$, we see that the logarithmic
derivative is given by
$${ {j_n^{\prime}}\over{j_n} } = S_n \,-\, { {j_{n+1}}\over{j_n} }\eqno{(20)}$$
and from (14) it is easy to derive a continued fraction for the ratio of
successive orders:
$${ {j_{n+1}}\over{j_n} } \,=\,
{ {1}\over{(2n+3)/x\,-} }\,\, { {1}\over{(2n+5)/x\,-} }\,\,
{ {1}\over{(2n+7)/x\,-...} }\eqno{(21)}$$
The coefficents in the denominators are $S_k + S_{k+1} \equiv T_k$ for $k$
starting at $n+1$. The equation is implicit in [10, 9.1.73 and 10.1.1];
it was derived (for Coulomb functions) as eqn(2.19) of ref.$^{7}$ and by a
different method in ref.$^{6}$
(It was quoted erroneously in ref.$^{23}$ in equations
(37), (38) and (39) -- in each case the first term should be dropped.)
A most important point is that (21) applies {\it only to the regular
solution.} It
is not apparent from the derivation that the formula {\it does
not hold } for the irregular functions $y_n$ or $h_n$, {\it even though
they, too, satisfy all the recurrence relations}. The explanation of this
remarkable result, and an indication of its generality, appears in
Gautschi$^{24}$ and in
refs$^{20,22}$: only the minimal ({\it i.e. regular}) solution has the
property (21). Also, combining (20) and (21) we find,
$$f \equiv { {j_n^{\prime}}\over{j_n} } \,=\, S_n \,-\, {
1\over{T_{n+1}\,-} }\,\,
{ 1\over{T_{n+2}-\,} }\,\, ...\,\, { 1\over{T_{k}\,-...} }\eqno{(22)}$$
The numerical evaluation of the continued fraction $f$ uses the fact that
eventually the denominator $T_{k}=(2k+1)/x$ will become large
enough so that the value of $f$ can be found by terminating (22) at
the $k^{th}$ step while retaining a chosen accuracy.
The problem of calculating the spherical Bessel functions becomes that of
computing the continued fraction $f$ to sufficient accuracy. Equation (22)
will be referred to as CF1, whereas CF2 will refer to a second continued
fraction: $p~+~iq~=~H^{\prime}/H$ which is discussed in \S6, eqn(34).
It may be noted that the reciprocal of (21) is also used, {\it e.g.} by
Lentz,$^{17}$ and that it is a more complicated expression.
\titlea{4.}{Evaluation of the Continued Fraction}
Continued fractions are mentioned in Abramowitz and Stegun$^{13}$ in
\S3.10 and are covered in Numerical Recipes$^{14}$ in chapter 5.2 .
There is an intimate relation between recurrence relations of the type (14)
and continued fractions: a full discussion and literature guide appears in
Chapter 4 of van der Laan and Temme,$^{20}$ particularly \S4.8,
and much valuable and detailed information appears in Wynn$^{21}$ and in
Wimp.$^{22}$ The forward evaluation suggested in
refs$^{13,14}$ is not to be recommended: a backward recurrence algorithm
originally due to Miller (see ref.$^{14}$ \S5.4 and \S6.4, ref.$^{20}$
\S3.3, or ref.$^{15}$ Chapters 4 and 5) is superior and is very widely used
in the improved form given by Gautschi.$^{24}$
Steed's method$^{6}$ for continued fractions is a stable forward recurrence
and it
has further advantages (discussed in ref.$^{23}$ \S4). It was adopted in
the author's earlier programs.$^{4,6,10,23}$ This method involves a {\it
summation} when updating $f_{n-1}$ to become $f_n$ and was hence subject to
rare numerical cancellation errors when a certain denominator approached
zero [ref.$^{10}$ eqn(11), refs$^{8,25,26}$]. Thus Steed's method for
continued fractions does not compute the result every where to uniform
accuracy. In cases where there are no zeros involved ({\it e.g.} in
calculating CF2 rather than CF1, see end of \S3)
then there can be no objection in principle to using it (as does BESSCC and
part of the COULCC code). The same method, although unnamed, has been used by
Gautschi and Slavik$^{27}$ and is also described and recommended in the
classic Gautschi paper.$^{24}$ Subsequent authors have not been alert to the
manifest advantages of Steed's method. Thompson and Barnett$^{11}$ in 1987
released the code BESSCC, which
treats complex-argument Bessel functions and evaluates the complex
continued fraction involved (CF2) by Steed's method.
The method of choice for continued fractions, however, is none of the above
but the {\bf method of Lentz}$^{17}$ together with the `zero shifts' of
both numerator and denominator discussed by Jaaskelainen and
Ruuskanen.$^{25}$ That paper, however,
proposed an elaborate change to to the algorithm, whereas Thompson
shows in Appendix III of ref.$^{12}$ how to achieve these shifts with minimum
change to the Lentz method.
Explicitly, the Lentz--Thompson algorithm reads:
\bigskip
\noindent {\bf $L$--$T$ algorithm for the forward evaluation of continued
fractions} \hfil\break
Given the $n^{th}$ convergent of a CF, {\it i.e.} the sum to $n$ terms of
$$f_n =\,b_0 \,+\,{ a_1\over{b_1\,+} }\,{ a_2\over{b_2\,+} }\, ...
\,{ a_{n-1}\over{b_{n-1}\,+} }\,{ a_n\over{b_n} } \eqno{(23)}$$
then evaluate $f = lim_{n\rightarrow\infty}\,(f_n)$ to an accuracy $acc$ with
the algorithm:
\bigskip\noindent $f_0:=b_0; ~~~if~ (f_0 = 0)~ f_0=small$
\hfil\break
$C_0:=f_0,~~D_0:=0$ \hfil\break
$for~ n=1,limit~ ~~~do ~begin$ \hfil\break
\indent$C_n:=b_n+a_n~/~C_{n-1};~~~~if~ (C_n=0)~ C_n=small$ \hfil\break
\indent$D_n:=b_n+a_n\times D_{n-1};~~~if~ (D_n=0)~ D_n=small$ \hfil\break
\indent$D_n:=1/D_n$ \hfil\break
\indent$\Delta_n:=C_n\times D_n;~~f_n:=f_{n-1}\Delta_n$ \hfil\break
\indent$if~ (\mid\Delta_n-1\mid<{\it acc})~ exit$ \hfil\break
$ end$
\bigskip{\it Notes:}
\item{1.} If at any stage we represent $f_n$ as $A_n/B_n$(before
cancellation), then the two quotients which are used in the $L$--$T$
algorithm are $C_n=A_n/A_{n-1}$ and $D_n=B_{n-1}/B_n$.
\item{2.} The parameter $small$ should be some non-zero number less
than typical values of the quantity $ acc\times \mid b_n \mid$ {\it e.g.}
10$^{-50}$ for typical double precision calculations in which a sensible
choice of $acc$ is $10^{-14}$.
\item{3.} The zero tests are to a working accuracy,
say $tol$, which can be chosen to avoid divide-by-zero error checks.
It could be taken as equal to $small$.
Thus algebraic conditions such as `$if\, (D_n =0)$' are to
be read as computation conditions:
`$ if\, (abs(D_n)0$,
the {\it do loop} to label 2 implements the downward recurrence (16) and
(17), until unnormalised values for $j_0,~j'_0$ are found. Relative
values of $j_L$ and $j'_L$ for each $L$ from this procedure
are stored ready for normalisation, which is
accomplished by (6). The values of $y_0,~y'_0$ are also
obtained from (6) and the recurrence relations (18), (19) used to find the
remaining $L$-values from $1$ to $Lmax$. It is clear that if a program
without derivatives were required then eqns(16), (17) could be replaced by
(15), with (21) used instead of (22, CF1). SBESJY was tested against
SBESJ,$^{23}$ COULFG,$^{4}$ the values in ref.$^{13}$ and the two programs
described below, DPHRIC$^{16}$ and cfbessel$^{26}$. The range of parameters
was $x=0.01-1000.0$ and $L=0-1000$, and in general, all results agreed to
better than a relative accuracy of $10^{-12}$ when the {\it accur} parameter
is set to $10^{-15}$. Some test data are given in \S8. The first two programs
are unusual in {\it not} using the functions $\cos (x) ,\sin (x)$ in the
calculation. Those which do may suffer from truncation error when a large
value of $x$ is reduced to the range $\pm\frac1/2\pi$ by subtracting $N\pi$,
where $N$ is a suitable integer. This reduction will be compiler dependent
and can be programmed, to some extent, by methods described in Cody and
Waite,$^{33}$ Chapter 8, p 136. The value $\pi$ is written as C1 + C2, where C1
= 3217/1024 = 3.1416 01562 50000 is exactly machine representable and C2 =
--8.9089 10206 76153 73566 E--6.
For the program RICBES both the two recurrences and the continued fraction
are now evaluated as Coulomb functions. First CF1 $=f$ is calculated from
(33) with $\eta=0$ (and hence $R_k^2=1$ and $T_{k}=(2k+1)/x$), and with
$n=Lmax$. It differs from (22) only by $1/x$. This yields a value
$F'_{Lmax}=\,fF_{Lmax}$ when the unnormalised
value of $F_{Lmax}$ is taken as unity. Then, if $Lmax>0$, the {\it do loop}
to label 2 implements the downward recurrence (26) and (27) until
unnormalised values for $F_0,~F'_0$ {\it i.e.} $\psi_0,~\psi'_0$, are found.
Relative values of all the $F_L$ and $F'_L$ are stored ready for
normalisation, which is accomplished by (13). The values of
$\chi_0=-G_0,~\chi'_0=-G'_0$ are also obtained from (13) and the upward
recurrences (29), (30) used to find the remaining $L$-values.
A recent program to compute $j_L(x)$ and $y_L(x)$ --called $n_L(x)$--
is that of Gillman and Fiebig.$^{16}$
Constructive criticisms of its archaic and awkward style have been given by
Welch$^{28}$ who has presented rewritten versions in order to emphasise that
FORTRAN 77 and FORTRAN 90 programmers
have no excuse for producing non-structured programs. The method of Gillman
and Fiebig is to calculate modified functions $u_L(x),~v_L(x)$ in place of
$j_L(x),~y_L(x)$ from which factors of $x^{L}/(2L+1)!!$ and
$-(2L-1)!!/x^{L+1}$ have been extracted. These scaled functions
tend to unity as $x\rightarrow 0$ ({\it cf} \S2 between eqns(6) and (7).)
A version of Miller's method is used to obtain the relative values of $u_L$
for all the $L$-values considered (0--1000), with $x$ in the range
0.01--100.0. The values are normalised by the Wronskian relation for the
$u,~v$ functions.
The reason for choosing $u_L(x)$ and $v_L(x)$ is that computational overflow
and underflow conditions are removed. In the days when representable real
numbers were restricted $10^{\pm38}$ or $10^{\pm70}$ then this was indeed a
problem. It is commonplace for modern FORTRAN 77 compilers to offer
$10^{\pm308}$ in double precision variables (as does the Lahey F77L
compiler, whose version 4.0 was used in this work) which is large enough
to remove the need for special programming in normal physical applications.
A second point is that the $u,v$ recurrences are {\it assumed} to be stable
in the same directions as are the $j,y$ recurrences. This is not an obvious
fact ({\it e.g.} Lentz$^{17,26}$) but the stability was proved
rigorously in 1980 by O'Brien$^{29}$ for real arguments. He claimed, further,
that no Miller's method is necessary, replacing it by two evalutions of `a
rapidly convergent series'. (No numerical details were given and O'Brien's
results have not yet appeared as a program.)
The third point is that tests for overall accuracy in ref.$^{16}$ only
address numerical consistency, and not the method. This is recognised by
Gillman and Fiebig who say that `it is of some value'. It is rather puzzling
to a user who sees tests aimed at a relative accuracy of $2.10^{-8}$
apparently produce results to better than $10^{-15}$, which is about machine
accuracy. One reason may be the tests measure the {\it rms} Wronskian
deviation from unity, which reduces by the square root of the number of
$L$-values, and here by a factor of 30. But the {\it same Wronskian} is
used to normalise the $u_L$-values, so the tests are not independent. There
is no need for this choice of normalisation for the exact value of $u_0(x)$
is $\sin(x)/x$ and should be chosen instead. However, the choice of this
correct
normalisation gives exactly the same result to machine precision. The
correct answer involves a more subtle point. By choosing to compute the rms
deviation {\it over all $L$-values} Gillman and Fiebig minimise the fact
that the downwards recurrence of Miller's method is most at error {\it at the
high $L$-values} which is worst when $x/L$ is largest. An independent
calculation$^{26}$ shows that $j_{1000}(x)$ is in error by a fraction
$10^{-3}$ for $x=100$, which falls to $10^{-7}$ for $x=0.5$. The healing
step in $L$ is small: even for $L=999$ the errors improve to $10^{-5}$ and
$10^{-15}$. The last few $L$-values should not be used, and it would seem
that the choice of 1000 for the starting $L$ is too large and their estimate
of an error $\epsilon \simeq (x/2L)^{6})$ is optimistic for $L$ itself.
Making the correct choice of a starting $L$ is the hard part of Miller's
method. Of course these comments relate only to the regular solution
$j_n(x)$; for the $y_n(x)$-values
depend only on their (accurate) starting values $y_0(x), ~y_0'(x)$ since
the upward recurrence is stable. The value
of $j_{1000}(100)$ is $5.32338~16172 \times 10^{-872}$, while
$j_{1000}(0.5)~=~6.06344~55462 \times 10^{-3172}$, and are clearly of
mathematical rather than of physical interest.
It should be noted that two of the tables in the Gillman and Fiebig article
contain sign errors, which do not occur when the published program is run.
Specifically, in Table III, for $j_L(100)$ the entries for $L=2,3,4,6,11,12$
and $16$ have the wrong sign, as do the $y_L(100)$ entries for $L=19$ and
$L=50$. In Table V, $v_{30}(100)$ also has the incorrect sign.
The approach of ref.$^{16}$ in factorising out the small-$x$ behaviour,
when reinforced by theory,$^{29}$ works well for spherical Bessel functions.
The underflows and overflows in reapplying the normalisation to recover
$j_L,~y_L$ can easily be trapped by extracting, say, powers of $10^{\pm200}$
when appropriate, as is done by Lentz$^{26}$. There is nothing to prevent
the CF1 calculation in SBESJY replacing Miller's method to remove the
inaccuracy near the maximum $L$-value.
The program $cfbessel.for$, published by Lentz,$^{26}$ calculates a single
value of $j_n(z)$, for complex $z$. His algorithm for the forward evaluation
of continued fractions creates an infinite product which is terminated when
a given accuracy is reached. It is a flexible and elegant concept,
applicable to a wide range of complex arguments, although just what the
limits are is not discussed by Lentz. Despite the different appearances it is
the same algorithm which is described in \S4.
As given, the Lentz program only results in $j_n(z)$ and not $y_n(z)$ or the
derivatives, and intermediate $n$-values are not calculated. Thompson's
program SBESJH given in ref.$^{11}$ includes all these features. These two
codes, together with the code {\it sphbes} of Press and Teukolsky,$^{34}$
will be compared in detail in a subsequent publication.
\titlea{6.}{Recurrence Relations for Coulomb Functions}
Each of the four Coulomb functions $F_n,~G_n,~H^{+}_n,~H^{-}_n$ (see
eqn(19) obeys recurrence relations [ref.$^{13}$ 14.2.3, 14.2.1 and 14.2.2]
which connect the functions of order $(n-1)$, $ n$ and $(n+1)$, thus
$$R_nw_{n-1} \,-\, T_nw_n \,+\,R_{n+1}w_{n+1} = 0\eqno{(25)}$$
$$w_{n-1} = [S_nw_n \,+\,w^{\prime}_n]/R_n\eqno{(26)}$$
$$w^{\prime}_{n-1} = S_nw_{n-1} \,-\, R_nw_n\eqno{(27)}$$
The coefficients are:
$$\eqalign{ R_k&\,=\,\sqrt{1+\eta^2/k^2} \cr
S_k&\,=\,k/x +\eta/k \cr
T_k&\,=\,S_k\,+\,S_{k+1}\,=\,(2k + 1)[x^{-1}+{ {\eta}\over{k^2+k} }] }
\eqno{(28)}$$
Equations (26) and (27) connecting two successive orders and a derivative are
in a form suitable for {\bf downward recurrence} in $n$.
The equivalent expressions for {\bf upward recurrence} in $n$ are
rearrangements of (27) and (26),
$$w_{n+1} \,=\, [S_{n+1}w_n \,-\, w^{\prime}_n]/R_{n+1} \eqno{(29)}$$
and
$$w^{\prime}_{n+1} = R_{n+1}w_n \,-\, S_{n+1}w_{n+1} \eqno{(30)}$$
The recurrence relations (25), (26), (27) are an alternative way of
expressing the differential equation (4) as a difference equation. As
Fr\"{o}berg$^{30}$ shows, not even two of these recurrence relations are
independent.
The equation analogous to the spherical Bessel (15) loses its simplicity and
becomes
$$S_{n+1}R_nw_{n-1} \,-\, T_nw^{\prime}_n \,-\, S_nR_{n+1}w_{n+1}
= 0,$$
and this is not particularly useful.
The boundary between solutions of oscillating and monotonic character occurs
at the turning point $x_L$ of (2), when $L$ is fixed, or alternatively, for
fixed $x$ it is found at $L_{TP}$. This is also given by (2) and equivalently
by
$$x^2(R_L^2 -S_L^2) ~=~1 $$
or
$$L_{TP}\,=\,(x^2-2\eta x+\frac1/4)^{1/2} -\frac1/2$$
Since $F_L(x)$ decreases to zero as $L$ increases beyond $L_{TP}$,
the stable direction (the
function must not decrease) is {\it downward recurrence in} $n$
to calculate the values of $F_n(x)$ ($n = L,L-1,...,3,2,1,0)$. Similarly {\it
upward recurrence} must be used for $G_n(x)$. In the region where $L