COULCC: A continued-fraction algorithm for Coulomb functions of complex order with complex arguments

Computer Physics Communications 36 (1985) 363-372 North-Holland, Amsterdam

363

COULCC: A CONTINUED-FRACTION ALGORITHM FOR COULOMB FUNCTIONS OF COMPLEX ORDER WITH COMPLEX ARGUMENTS I.J. THOMPSON Daresbury Laboratory, Warrington WA4 4AD, UK

and A.R. BARNETF Physics Department, The University, Manchester M13 9P1, UK Received 14 June 1984; in final form 5 December 1984

PROGRAM SUMMARY Title of program: COULCC: Complex Coulomb & Bessel Functions Catalogue number: ACDP Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue) Computer: AS/7000

Installation: Daresbury SERC

CRAY4 CDC 7600

ULCC UMRCC

uncharged scattering (spherical Bessels) and cylindrical Bessel functions are special cases which are more easily solved. Two linearly independent solutions are found, in general, to the differential equation f”(x)+ g(x)f(x) = 0, where g(x) has 2 terms, with coefficients 1, —2~and — X(X x°, and x~ + 1),x~ respectively. Method of solution The continued fractions of Steed [1,2J are supplemented, if necessary, by a 1F1 series expansion or by Padé accelerations of a 2F0 asymptotic expansion. Recurrence relations are used for integer steps in A in a stable manner. It should be noted that

Operating system: OS/MVT with VS FORTRAN

the routine will solve for a single arbitrary A without recurrence, if required. For small x a previous restriction on accu-

Programming language used: Fortran 77

racy has been removed by adding a subroutine to evaluate the irregular solution (singular at x = 0) by a suitable combination of series.

High -speed storage required: 822 Kb for compilation; 202 Kb for execution No. of bits in a word: 32 (4 bytes of 8 bits) Peripherals used: card reader or VDU, pnnter No. of lines in combined program and test deck: 1658

Restrictions of the complexity of the problem On the Coulomb bound-state poles the functions are singular (from their boundary conditions). The program returns the residue polynomial but only one solution exists, and it is found. Typical running time A direct comparison of COULCC and its predecessor for real

Keywords: Coulomb, Bessel, Whittaker, hypergeometric, continued fraction, scattering, closed channels, off-shell, resonances, reactions, Regge poles Nature of the physical problem The routine COULCC calculates both the oscillating and the exponentially varying Coulomb wave functions, and their radial derivatives, for complex i~(Sommerfeldparameter), complex energies and complex angular momenta. The functions for

arguments shows an increase by a factor of 2 for the new code. The test deck (comprising 36 test values and excluding the error condition) took 1.14 s for execution on the NAS 7000 and 2.2 s on the CDC 7600. References [1] A.R. Barnett, D.H. Feng, J.W. Steed and L.B.J. Goldfarb, Comput. Phys. Commun. 8 (1974) 377. [21 A.R. Barnett, Comput. Phys. Commun. 27 (1982) 147.

OO1O-4655/85/$03.30 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

364

If. Thompson, AR. Barneu

/

Complex Coulomb function

LONG WRITE-UP 1. Introduction This program, COULCC, is the complex generalisation of a continued-fraction algorithm for calculating real Coulomb wave functions which has been described and evaluated in a series of papers [1—7].It is designed to find linearly independent solutions F, G, H~and H of the equation dX

+

—

2~ A(A + 1) )f(x) —

X

0

subject to the boundary conditions

Steed’s ~method [1] by a series expansion for the regular function near threshold. We also supplement the program by a Padé-accelerated asymp-

F~(i~, x) = 0, =

—~

Hx±(~q, x)

X

—~ -~ 00

~

0) =

~

x) x—.00 —b cos

+

~ 0A’

e±IOA,

~A~rr+ a~ and

—

a~is the

F

—

totic expansion for the Whittaker function as in ref. [15] and find that we are now able to include all of theinabove complex general. physical cases with x, ~ and A 2. Method of calculation

where 0A = x ~1n(2x) Coulomb phase shift F(1 + x + i~)}1/2 e’°”= [r~i + ~. —

—

chosenfor tions instead general to scattering continue calculating problems. We the stanhave dard Coulomb wave functions, but supplementing

X

G~(~j, 0) = ~ 0) and F~(ii,x) x—~00 sin ~

Bell and Scott [11] and for Im( ~)< 0 by Hebbard and Robson [12]. Solutions near zero energy are usually found, on the other hand, by ‘Z-scaled coordinates’ r = qx and c = i/q so that they can be calculated, e.g. by Seaton [8] or Curtis [13], continuously across the energy threshold. The resulting wave functions are not dimensionless as are F, G, etc., and Humblet [14] has proposed similar wave func-

.

For positive-energy scattering problems, x is real and positive, i~and A are real, and the real F and G functions resulting can already be calculated to adequate precision, typically io~12[1—3]. When looking for example for off-shell structures such as resonance poles, x and ~ have small imaginary parts, so F and G are now complex in general. Tamura and Rybicki [9] have considered this case, and it is straightforward to extend COULFG [7] by complex continuation to treat off-shell regions near the physical axis. Another common analytic continuation is to complex angular momentum A, to trace Regge poles, etc., as have been calculated by Takemasa al.~[10]. Negative energy solutions (i.e. x et and imaginary) are also required when seeking bound states, or when closed channels arise in multichannel scattering problems. Traditionally just the decaying Whittaker function is found, for Im(i~)>0 by

The general formulae for defining the Coulomb functions F~(i~,x) and G~(ij, x) for complex A, ~ and x are described in a related paper [16] by the authors, along with a description of the different combinations of continued fractions used in different (x, s~,A) regions. The power of the method derives principally from the evaluation by Steed’s method of following continued fractions F’ ~÷, R~÷, R~~ 2 CF1: = T~+1— T~÷2— —

—

where

s~

+

=

~,

T~=

s~+ s~,

X

and 2~ 1 R and

+

CF20: p .

1(0(1

i 2/A2 1 +

iwq ~

=

——~+—

Hr’/Hr

i~

ab

x~’ x 2(x—ij+iw)+

I.J. Thompson, AR. Barnett

(

(a + 1) b + 1) 2(x i~+ 2iw) + —

...

‘

where a = iwij A and b = io~q+ A + 1, with the coefficients given = +1 or w = —1. analytically as indicated and w Steed’s algorithm, derived in ref. [1], is a for—

ward recurrence method for evaluating a continued fraction of the form a1 a2 a,, h = a0 + b+ b+ b,, + ...

by carrying forward just D,,, the ratio between successive denominators B,, in the equivalent polynomial ratio h A,,/B,,, and thus avoiding possible overflows in evaluating A,, and B,, separately. The method is h0 = a0,

=

for n

Complex Coulomb function

365

These analytic continued fractions are supplemented in ‘case 1’ for small x by a direct evaluation of a 1F1 series to give

~

1~C~( F~(n~x) = xA + 1 e ± X 1F1(1 + X ±in; 2X + 2; R2ix). This value of FA together with f from CF1 gives FA’. The upper or lower signs throughout are chosen to given 1 + A ±iij the most negative real part possible, to enable quicker and more accurate convergence. The 1F1 series is absolutely convergent, but to extend the maximum values of n and x which are possible before large partial sums cause large fractional errors in the result, the evaluation is performed in extended precision. Inout COULCC, REAL *16 variables are used, writing the real and imaginary arithmetic separately for use on

1/b1,

machines such as the CRAY-i which has no

a1D1, h0 + ~h1,

COMPLEX * 32-type variables. The corresponding evaluation of GA in case 1 is still achieved with CF2; whereas for yet smaller x

=

h1

/

=

2 upto ~ max do begin

the computation of CF2 becomes too inefficient

=

1/(D,,_1a,, + b,,),

and cases 5 and 6 are used instead. In these cases we re-expand the irregular functions in terms of the regular solution. If 2X is non-integral (case 5),

= =

(b,,D,,

h,,_1

+

—

i)i~h,,._1,

~h,,,

G~is found by means of the 1F1 series for F_A_I, whereas if 2X is an integer (case 6), G~is given by the logarithmic expansion (eq. 13.1.9 of ref. [17]) analogous to the Neumann series for Bessel func-

if I i~h,,I < ~Ih,, then exit end If b1 = 0, then the continued fraction is reformulated as

tions. For large-x values (case 4) asymptotic expansions based on ~ x) = e~°A 2F~(ico~A, icoi~ + A + l;;w/(2ix)) are used, either directly for HA(co), or in combination to give CF1A = (H~’ H~’)/(H~— H). For intermediate-x values, the above asymptotic expansions are accelerated by Padé techniques, a process which has been found [15] to extend inward the region of convergence to limits set only by machine precision and exponent range. We use the F-algorithm of ref. [18] to calculate term-by-term the coefficients of the continued fraction ‘corresponding’ to the original series, in parallel with a term-by-term evaluation of this continued fraction by Steed’s method. The Padé acceleration is only invoked if the asymptotic expansions are not sufficiently accurate on their —

h

a0 + ~ + a4 a2 b3 + b4 + The continued fractions CF1 and CF2 form the heart of the method and above their sole use is labelled as case 2 and case 3 within the program. Case 3 is the real-parameter case and it hence corresponds to the parent program COUFLG [7]. Case 2 is for complex parameters but for the region where the resulting wavefunctions still retam oscillatory character. The absolute normalisations are determined from the Wronskian equation FX’GX GX’FX = 1. For this to hold all parameters need to be close to their real axes. Other regions of the parameters relate to the 4 other cases described below. =

~.

—

...

—

366

I.J. Thompson, A.R. Barnett

own. Since we have an alternative method for smaller x, namely the 1F1 series, we may set the maximum length JMAX of the continued fraction at only 50. We see from fig. 1 of ref. [16] that this includes the region where there is a relatively reliable correlation between the numbers of iterations for CF2~’~~ and for ~ itself, The code is set up to calculate the regular solution FA and one of the irregular solutions GA, H~= G~+ iFA, or H~= GA iFA according to a MODE variable. In general the user should ask for the linear combination with the smallest absolute magnitude, as it is impossible to obtain this solution accurately from the others due to the excessive cancellation. For large A, G and H ± are all large, whereas FA is small so is always returned. —

+

.

For small A and x complex H will usually be small for Im(x)> 0, and H for Im(x) < 0, so the user should choose H ± according to the corresponding half-plane containing x. For small A and x near the axis, G and H~are all oscillating at about the same magnitude, so the choice of MODE is not critical here. The code also calculates the Coulomb phase shifts GA(n), and these can be returned if required. This is strongly advised if either of 1 real + A part, + in or 1n is likely to have a negative as 1then + A the cut in o~(fl) is approached. On the cut all the Coulomb functions change sign together, and there will be a simultaneous change of ii in ~A(fl)• If the Whittaker functions are to be calculated for example by —

/

Complex Coulombfunction

polynomials, and the code uses the zero returned by the log-gamma routine to calculate these polynomials scaled to finite values. A corresponding “G~” is returned such that the Whittaker function calculated by the formula above is still correct, even though the gamma routine printed an error message. The above two paragraphs imply that if there is any choice of 1 + A ±~n approaching the negative real axis (e.g. for Coulomb eigenstates), the user should calculate the Whittaker function from H~ by means of the phase shifts returned at the same time, and use the Whittaker function in subsequent calculations. 3. Code descnphon The calling sequence for the program is (see fig. 1) CALL COULCC(X,ETA,ZLMIN,NL,FC, GC,FCP,GCP,SIG,MODE1,KFN,IFAIL) where the arguments have the following type and meaning. X ETA ZLMIN NL

complex complex complex integer

x ~ 0,

number of A values Amin, Amin + 1, Amin + NL 1 calculated; complex array dimension NL, regular ..

—

FC W_I~.x+l/ 2(—2iwx) GC

1

=expico[~’n(A+icon)—aA(n)JHx(q, x) with co arg X>

FCP

— 11/2,

then the sign of W is only determined if the GA(fl) used in this formula is that returned by the same code which calculated the HA. A consistent treatment of the cut is then maintained. If one of 1 + A ±in is zero or a negative integer, then a gamma function in the expression for a~(n) will be evaluated on its pole. These are just the hydrogemc Coulomb bound-state poles, and with these states &°~,F~and ~ are all, strictly speaking, infinite. The residues, however, are finite .

GCP SIG

solution, FA, complex array dimension NL, irregularsolutionGA orH~, complex array dimension NL, regular derivative F~= d FA/ dx complex array dimension NL, irregular derivative G~or H~, complex array dimension NL, Coulomb phase shifts if KFN = 0

MODEl integer: IMODE1I determines selection of F, G, H ~, and MODEl <0 indicates that exponential scaling is to be used.

If. Thompson, A.R. Barnett

x

/

Complex Coulomb function

SUBROUTINE COULCC(XX,ETA1,ZLMIN,NL, FC,GC,FCP,GCP, SIG, MODE1,KFN,IFAIL)

C

CccccccCccccccccccccccccccccccccccccccccccccccccccccccccccccCccccccccCcc C C

C C

COMPLEX COULOMB WAVEFUNCTION PROGRAM USING STEED’S METHOD

C C C

C A. R. Barnett modified I.J. Thompson

Manchester March 1981 C Daresbury, Sept. 1983 for Complex Function. C

C C C C C C C

original program

RCWFN in + RCWFF in + COULFG in description of real algorithm in description of complex algorithm this version written up in

CPC CPC CPC CPC JCP CPC

C C

8 (1974) 377—395 11 27 21 XX XX

(1976) (1982) (i98i) (1984) (1984)

141—i42 147—166 297—3i4 YYY—ZZZ YY~—ZZZ

C C C C C

C C C C C C

C COULCC returns F,G,G’,G’,SIG for complex XX, ETA1, and ZLMIN, for ML integer—spaced lambda values ZLMIN to ZLMIN+NL—i inclusive, thus giving complex—energy solutions to the Coulomb Schrodinger equation,to the Klein—Gordon equation and to suitable forms of the Dirac equation ,also spherical & cylindrical Bessel equations

C C C C C

C

C

C C C C

if /MODE1/— 1 get F,G,F’,G’ for integer—spaced lambda values C — 2 F,G unused arrays must be dimensioned in C — 3 F, F’ call to at least length (1) C -4 F C — ii get F,H+,F’,H+’ ) if KFN—O, 11+ — G + i.F ) C — 12 F,H+ ) >0, 11+ — J + i.Y — H(i) ) in C — 21 get F,H—,F’,H—’ ) if KFNO, H— — C - i.F ) CC C — 22 F,H— ) >0, H— — .1 — i.Y — H(2) ) C

C C C

C C

C

C C C C

if MODEi 0 ( IMAG(XX) if MODEl < 0 4 KFN 4 3 ( REAL(XX) if MODEl < 0 & KEN — 3 then FC — EXP(—ABS(SCALE)) * ( F, j, J, or I) and CC — EXP(—ABS(SCALE)) * ( C, y, or Y ) or EXP(SCALE) * ( H+, H(L), or K) or EXP(—SCALE) * ( H— or H(2) )

C C

C C C

C C C C

if

KFN

— —

C

—

C C

—

C

C

C C

F & C j & y

C C

2 cylindrical Bessel 3 modifIed cyl. Bessel

J & V I 8 K

C C C

“

“

“

“

“

“

and where Coulomb phase shifts put in SIG If KFN—O (not —1) C The use of MODE and KFN is independent (except that for KFN3, H(1) & H(2) are not given)

C C C C

With negative orders lambda, COULCC can still be used but with reduced accuracy as CF1 becomes unstable. The user is thus strongly advised to use reflection formulae based on H+—(ZL,,) — H+—(—ZL—i,,) * exp +—i(.Ig(ZL)—sig(—ZL—i)—(ZL+i/2)pi)

C C C

C

C

C C

C

0,—i complex Coulomb functions are returned 1 spherical Bessel “ “ “

C C C C

C C C C

C C C C C C C

C C C C

C Precision:

results to within 2—3 decimals of ‘machine accuracy’ • but If CF1A fails because X too small or ETA too large the F solution is less accurate if It decreases with decreasing lambda (e.g. for labda.LE. -1 & ETA.NE.0) REEK in COMMON/STEED! traces the main roundoff errors.

C C C C

C C

367

368

I.J. Thompson, A.R. Barnett

/

Complex Coulomb function

C C C

COULCC is coded for real*8 on IBM or equivalent ACCUR >— iO**_14 C with a section of doubled REAL*16 for less roundoff errors. C (If no doubled precision available, increase JMAX to eg 100)C

C C C C C C C C C C

Use IMPLICIT COMPLEX*32 & REAL*16 on VS compiler ACCUR >— 1O**_32 For single precision CDC (48 bits) reassign REAL*8_REAL etc. IFAIL

on input

—

0

:

ne 0 : IFAIL

in output

—

—2

:

—

—i

:

0 : ge 0 : —

C C C

—

—3

:

C

C C no printing of error messages C print error messages on file 6 C argument out of range C one of the continued fractions failed, C or arithmetic check before final recursion C All Calculations satisfactory C results available for orders up to & at C position NL—IFAIL in the output arrays. C values at ZLMIN not found as over/underflowC roundoff errors make results meaningless C

: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Fig. 1. Subroutine COULCC, introductory comment lines. —

—4

If IMODE1I = 1 F, G, F’, G’ = 2 F, G = 3 F, F’ = 4 F = 11 F, H~,F’, H~’\ if KFN ~ 0, = 12 F, H~ H ±= G ±iF = 21 F, H, F’, H’ i~KFN > 0, =22 F, H~ H~=J±iY

1

1

2 3 —

1

3

where F is in the FC array, G, H~or H is in the GC array and where the definition of H ± is fixed by KFN (see below). If MODEl <0, then the values retained are scaled by an exponential factor (dependent only on x), designed to bring nearer unity the functions for large lxi and small A, n: Define s = imag(x) for KFN 2 and s = real(x) for KFN = 3, then FC = e~I5I{F,j, J or I) (e~(G, y or Y) or and GC = e~{H~or K or

}

spherical jA(x) YA(X) Bessels or h ~= j ±iy complex cylindrical J~(x) YA(x) Bessels or H ±= .J ±i Y modified cylindrical IA(X) KA(x) Bessels Coulomb functions but without using the SIG array. complex

The FCP and GCP arrays contain the derivatives of the appropriate Coulomb or Bessel function with respect to x. If KFN = 3, then H ± are not defined, and mod(MODE1, 10) is used in place of MODEl. If KFN 3, the use of KFN and MODEl is independent. IFAIL: integer on input ~‘

=

0 no printing of error messages,

#

0

print error messages on

file 6;

on output = 0 no errors occurred, # 0 error exit. The error exits are further classified: IFAIL = —2: argument out of range;

KFN integer, determines kind of Coulomb or Bessel functions: KFN 0 complex Coulomb functions

FC: GC: FA(n, x) Gx(n, x) or H~=G± iF

=

=

—1:

the continued fractions failed or arithmetic check, before any results calculated; —3: the functions of Amm could not be calculated (but perhaps some at A> X,,,,,,); one of

I..!. Thompson, A.R. Barnett C C C C C C C C C C C C C C C C C C C C

/

Complex Coulomb function

369

C C C ACCUR target bound on relative error (except near 0 crossing.)C (ACCUR should be at least 100 * ACC8) C ACC8 smallest number with 1+ACC8 ,ne.l in REAL*8 arithmetic C 5I6 arithmetic C ACC16 smallest number with i+ACCI6,ne.i in REAL FPMAX magnitude of largest floating point number * ACC8 C FPMIN magnitude of smallest floating point number / ACC8 C FPLMAX LOG(FPMAX) C FPLMIN LOG(FPMIN) C C ROUTINES CALLED : LOGAM/CLOGAM/CDIGAM, C F20, CF1A, RCF, CF1C, CF2, Fli, CF1R C Intrinsic functions : HIN, MAX, SQRT, REAL, IMAG, ABS, LOG, EXP, (Generic names) MINT, MOD, ATAN, ATAN2, COS, SIN, DCMPLX, SIGN, CONJG, INT, TANH C Note: Statement fntn. NINTC — integer nearest to a complex no. C C Parameters determining region of calculations : C

Machine dependent constants :

C

C

C C C C C C C

R2O ASYM XNEAR LIMIT JMAX NDROP

estimate of (2F0 iteration.)/(CF2 iterations) C minimum X/(ETA**2+L) for CFLA to converge easily C minimum ABS(X) for CF2 to converge accurately C maximum no. iterations for CFI, cF2, and iFI series C size of work arrays for Pade accelerations C number of successive decrements to define instabilityC C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC Fig. 2. Machine dependent constants and subroutines used.

>

0:

arithmetic error during final recursion, so that results are available in the FC, etc., arrays up to and including position NL IFAIL; —

=

—4: roundoff errors

4 >

VACCUR.

The criterion for convergence of the internal expansions is consistency to within one half of the value of ‘ACCUR’ (see fig. 2). In the supplied code, ACCUR is set at l0_14, appropriate for REAL * 8 arithmetic on IBM machines. For CDC 6600/7600 systems set ACCUR = l0_12, and set it also to i0_12 on CRAY machines. The ACCUR value may be increased if these full accuracies are not required, and if faster calculations would be useful. If ACCUR is 10—8, for example, then (apart from the limiting cases described in section 4), the relative errors will be less than 40% of this figure except near zero-crossings of the functions. ACCUR is the target upper bound on the relative errors in the results and is deliberately set at about 100 times the unit roundoff errors.

Further information about the performance of COULCC is provided by a named common block /STEED/ containing in order RERR(REAL * 8) and the integers NFP, Nh, NPQ(2), N20, KAS(2). RERR is the estimated maximum relative error in the calculated Coulomb functions FC, GC. Then NFP

=

number of iterations for CF1, or

Nil

=

(numberofofterms terms required for CF1A); number for -

1F1 NPQ(IH)

=

number of terms required for CF2~

N20

=

or for IH CF2 = 1 and 2, respectively; number of terms required for 2F0

KAS(l)

=

expansion; case number, 1, 2, 3, 4, 5 or 6 for

KAS(2)

=

upper range of1,A 2, values; case number, 3, 4, 5 or 6 for

lower range of A values. The meanings of these case numbers are defined in the accompanying paper [16].

370

If. Thompson, A.R. Barnett

/

4. Accuracy and range of arguments As a rule, the results are obtained within 2 or 3 decimals of machine precision, except for the following range limitations: (a) negative A: The code can still be used, but with reduced accuracy as CF1 becomes unstable. The user is thus strongly advised to use A-reflection formulae based on ± HA± =H_A_lexp+i(aA—a .

—

A

1—(A+l/2)11).

—

(b) IA,,~I,JIm(n)I> 50: The convergence of the iF1 series and of the Padé-accelerated 2F0 expansions is now significantly reduced and may become inaccurate. (c) F-unstable Cases: In these cases CF1 cannot be used as the F solution decreases (rather than increases or oscillates) for a range of decreasing A, for example for x = lOOi and n = 50i. These, however, can only be detected if the number of A-values required, NL> 6, and in severely unstable cases only if NL is large enough so that the instability is nearly ended (in the above example F is decreasing below A Once the, instability is found, CF1A is used in preference, but if this occurs for n too large CF1A may not converge; this is the limitation of the program. The code only allows for one reversal of the direction of recurrences because of instabilities, so if negative A are also required the A-reflection formulae will have to be used before and after calling COULCC. The presence of poles (see(e)) also leads to reversals, so the program cannot handle both poles and instabilities, but should be called once for A = A mm to Im(n 1) and then for A = Im n up to

~

—

Am~.

(d) Small or medium x, if no REAL *16 arithmetic is available for the 1F1 series calculation: The partial sums in this series may be much larger than the final result if n is large, or if x is just below where asymptotic expansions may be used. For complex Bessels = 0) the 5error 5 for x(nbetween and would degrade up to i0 15.

Complex Coulomb function

(e) Near poles and zeros: The relative accuracy is determined principally by the accuracy with which the distance to the pole or zero can be calculated in the floatingpoint arithmetic. Poles and zeros are ap. proached whenever 1 + A ~nor 1 + A + in is near a nonpositive integer, and there are the usual zeros when real-valued functions change sign, etc. (I) Large x: The oscillating phase is as sin(x), and the —

accuracy in the results will decrease Just according to the accuracy in mod(x, 211).

5. Subroutines used COULCC calls the two routines CF1A and F20 as complex functions to use asymptotic expansions to calculate the regular logarithmic derivative f = F’/F and 2F0 values, respectively. If the expansions start to diverge, the two routines call RCF to calculate the corresponding continued fraction coefficients. The RCF code is copied from table 9 of ref. [18], but with a modified restart procedure and a change to complex variables. The routines CF1R and CF1C calculate the continued fraction CF1 for real and complex arguments, respectively, while the function CF2 evaluates the complex continued fraction CF2. The complex function Fil evaluates the 1F1 series for cases 1, 5 and 6, either in normal precision (with COMPLEX *16 variables) or in doubled precision (with REAL *16 variables). It can also calculate the 1F1-like series with the digamma factors that appears in the logarithmic solution for the irregular function. The code uses a combined CLOGAM and CDIGAM function to calculate the natural logarithm and the logarithmic derivative of the gamma function with complex argument, and with the cut along the negative real axis. Initially we used Kolbig’s program [19], but subsequently we generalised his methods and his code for arbitrary precision, subject only to machine accuracy and not subject to prestored floating point coefficients. In the test deck this precision is pre-set by a call from COULCC to the ENTRY LOGAM with one

If. Thompson, A.R. Barnett

REAL *8 argument ACCUR. The gamma function results decrease worsen by one or two decimals from machine accuracy in general, because of the number of recurrence steps required for very high accuracy calculations and the fixed collection of Bernoulli numbers. 6. Test deck The test deck contains a main program to call COULCC for a succession of x, i~ and A combinations determined by the data read in. The first 36 data cards are designed to reproduce published results of existing Coulomb, Bessel and Whittaker programs. The final set of 6 cards shows typical error conditions that may result if the input ranges are exceeded. References [1] A.R. Barnett, D.H. Feng, J.W. Steed and L.J.B. Goldfarb, Comput. Phys. Conunun. 8 (1974) 377. [2] A.R. Barnett, Comput. Phys. Commun. 24 (1981) 141. [3] A.R. Barnett, J. Comput. Phys. 46 (1982) 171. [4] A.R. Barnett, Comput. Phys. Comniun. 21 (1981) 297. [5] A.R. Barnett, J. Comput. Appl. Math. 8 (1982) 29.

/

Complex Coulomb function

371

[6] D.H. Feng, A.R. Barnett and L.J.B. Goldfarb, Phys. Rev. 13C (1976) 1151. [7] AR. Barnett, Comput. Phys. Commun. 27 (1982) 147. [8] MJ. Seaton, Comput. Phys. Cominun. 25 (1982) 87. [9] T. Tamura and F. Rybicki, Comput. Phys. Commun. 1 (1969) 25. [10] T. Takemasa, T. Tamura and H.H. Wolter, Comput. Phys. Comniun. 17 (1979) 351. [11] K.L. Bell and N.S. Scott, Comput. Phys. Commun. 20 (1980) 447. [12] D.F. Hebbard and B.A. Robson, Nuci. Phys. 42 (1963) 563. [13] A.R. Wave Functions, vol. 11, RoyalYork, Soc. Math.Curtis, TablesCoulomb (Cambridge Univ. Press, London/New 1964) ~. 9. [14] J. Humblet, Nuci. Phys. 50 (1964) 1; Ann. Phys. 155 (1984) 461. [15] C.J. Noble and I.J. Thompson, Comput. Phys. Commun. 33 (1984) 413. [16] U. Thompson and A.R. Barnett, J. Comput. Phys. (submitted). [17] M. Abramowitz, Handbook of Mathematical Functions. eds. M. Abramowitz and l.A. Stegun (Nat. Bur. Std., New York, 1964) p. 537; Tables of Coulomb Wavefunctions, NBS Applied Math. Ser. 17 (1952) [18] J. Patry and S. Gupta, EIR-Bericht Nr. 247, Eidg. Institut für Reaktorforschung, Wurenlingen, Switzerland (1973). [19] KS. Kolbig, Comput. Phys. Commun. 4 (1972) 221.

372

If. Thompson, A.R. Barnett

/

Complex Coulomb function Lu

U),Cr14 ,,

3. -,

GO

o’U’

)

C -~ 3. C

,‘--O ..‘4 •• 1(10’

C J 3, — C

..‘r

114

C,In1

3. U. 3.

O’0’ 3)114 Ifl1”. ‘O~ Cl’s -Sr..

(‘5 I..

0~ •‘ II

114 II

‘~.

(14

II

C C 3.

III

(IS

C

‘C C 0 ‘C C

0 Ca) C • C C

~

O

Z ~ ~

5(13113 1+ CC ‘0~ ,-C r’4~ p-I”. P-C lIla) 3.3. 01(4 310. ‘GO I-N.

...s

1-1 • CII C

C _4 3.4 $ C~ C C .0 C

(3 0.0 C • C

C 7 —~ 3’..41 (14 S CII C

CU) ‘ICl-I

.51

~.

~.t . . 1(115’. I

tn C)r C U CC .0 0 • 43

C C

— 3. 0 U .55 3. 1.’ U, 0 I” UI 144

—

2 Is. 3.

0-1 • • 0.0 II

3.

II II C 1.3. Is) •3. 0 0 C • C 1-Il .4 Cl” 0< 0115. 01.1 • C.’

CC 1111 00 0414. ‘“ •... CC COO CC CC

00 CC 37(4

II II C 1—0. 14,0. 04 sIX ‘3 C CC C • 1’ I-Il .4 C— DC CI,. 0’-’ • 0”

Or LI 0 .4 110 0 XCI

1111 .40 154(1’S

II

3.

II Ii

3.

Ii II

0 C C C ‘C C

C C00 ‘I • C N.

S

1130’ 0.4 .5(15 .0”. O’I’.l (IllU • ‘ I 1114

C —~ .

.4 £51 CII C

CU) ‘I’ CI—

III 0(4 14.14. ...‘ •... 00 CC CC CC

00 (40 CC

1111 .44) P.JP’S

.J

‘3O’U’(14U) 33 ~ F,t’?’.U).S 0~O.S1 (1J?W? 3..11111

Is. 0 I 3. C 0

Or.5r,J II

(13 II 3. U. 3. 14

II

3.3.3.3.3.

Lu

CCCC

Lu

~

lUll

‘3-3-3-3 +4+4

iii

+114

C000 Ø.a)~..,CU)~-rIflP-.*W)

3.

•

CU)4’I CO>.C 0C,.3

I 3. LI)

lilOGa) I’4CU)r.. 4’I4U).C

U-

O.nJØ..Ø ‘0,-CO

)-

r..C000

U-

~‘)

3~10.0 (1-InCa)

‘3

00-3?’.

“3

— •. C (14 LS~ .4

or-CC p.,545.JU’ OF-SO a)U)~’3.•.. • rYDfl I II

_I

Ce-a),’ Q.*0~ ~fS.(..~CII.,-,‘~)~F,’0(14 • CC-SII~I

‘tU)O.In

(‘4fl4flJr,, 00CC’ jIll 0000 LflIn4.) (13).-,-

000114 +4+4 03)00’ 3.000’ 3),-niU) ‘0a)(13 ((.000’ ‘0093) 3)I-1U’I’4 ~In.0O~ (1’IO’CU’ 0’3).C.?•(13r0C’.(1401(115’ • • •. (‘11-0(0 liii

II

11411411

3. ‘3(43.

00CC

0r40U)

(14

~

00CC

.~ • 3. 3. U)

0. 4J 3. 3. C U

,-,-C14,I

Lu C C 3.

~

~.‘~‘4.3 CN.4--* C0U’~’. U)U’OO1(’ • • (‘J~’Jr’.(1J ~ul

3. , 3.

U)

Or.rO

‘0U)U’S’?

II

II)

U)

~ Cc.*.fl

(14

Lu

U)U)OU)

C4~.Or.CU’O.O.

“3

a)

U)

.. ,J

rr(’4e 00CC lUll 0000 .U.-O’U) .0(144,>5004 3)InU)(1. 0,”C, .fl’0~C >54O’~ (134(144 InU”OO. ‘Q0”C,’3. ‘N. (‘4 ~C(.I4-i

111111111 Lu C C 3. l~’ 113

II II II II II 000(30 0(7(7(417

C.000N.1’JInlfl

l’JCl’JcO,-I135’)C

C _I Z’3

C C C

0 C CC

.4 C (3

.4 C .0

.4 CS’S. ‘34 1-LI. 0”, • 4..

0

CI

Ci

0

—

11111111 000L) ).)C(4(7

,In0’In14’. O.0a)O 0~(1~CIs ~O’3)CU ,‘-.344 ~tN.3)U) C’O,-’3 -* .OOU’ 000,0 ‘... l~J-4na) II U

Ii 2 .,. 3. l’s

II

11111111 Lu C 0 3.

U).OflsO. ‘OIn1’4U) CCC’?’. aIsnrO.• •‘ r-rrU)

— 3. F-s A. 4

,-000 44-34 +4+4 0000 (0(13,53) ,-InO’U’.-,-

r’.C.4,0I(.

II

ii

3. U. 3.

.)LflCC O0.~C14 (1-4-ICr-s ,I’U)1’15’0

3. Is. 3.

(1J

04-ICr.-

1’..

I

II

II

11111111 U) C 0 3.

3.0.)->’

0)

(‘4,-CO. — •• 3. (14 04 .4

11111111 Lu C C 3.

‘. C(0 C C •‘0 C

0 CC .— • C

C -11”. C.’.4 (‘4 CII C

CU) 1-C) U) C,_ 0 I 00 .0) C. • 3.

0 0 C C ‘C 0

5-” 1(lflS1’JflJ 00530 1111 0000 ~SL(1-153Q. In$.(0’C U’0In$ 001530’

00 (‘SF3 (1• C

r’400$4 .01(1,-C ,-O’C,.01-1(4(13)90.3) 014.51(1113 ~

1’. IIPS flJ(UC5LIflC OCCO.11111 00000 Ql”lI’—4’3 5(10(1(4,-OC 3)3.143310’ (‘IInU’InU’ C’0111In0 004.03.

4(0(1-3.43 (‘.004-0.5 ‘00.1,14,0’0~N.l’( (‘.O.U)NC 0.0300’’3.3)0043 I I 1111111111

Lu C 0 3.

>3.).>

—~ XI’) .4 £55 S CIII C

CU) ‘I’ 0.-’ 5(1 00CI I CC .0 C • 5(1

OU’(’-r“.SflO~’ .30’InAI’I .O’OISJC 11)113(1.51’ ,‘,C)fltfl • • • (‘1,- r.0 I 11141111

.-lS-1I”ll-I..l

Ii

0 C (4

C) C 114

0 C. (‘4

0 .4 (‘4

C

.‘

WI

C r

II

II

II

II

.4 N.

.4 151

.4 (115

.4 (‘.5

‘4 113 • I1-Il .4 I).’, CC CU. 1145-I • (‘5..

C) 0 .4 11= CI >10

IIIIIIIIII

l-Il-l-0-4

00000 04041141(04

00004) CCC CC 00000 000CC

00000 C.OCCC ,-‘,‘,‘rr 01—1000 I-In 1111111111 .4.40.4.4 P’J15dr14F.4N.

11111111 04 -3. 0 0 C (I, • 0’ II

01.4(10 04040404

.4 0-I CC CU. 0”. • “3..

0000

00.4(4 CC CC 000(5

31003110 l’s1’4l’4I’4 (‘-(‘-(1-15.

C C C C .0 C

C 013. ,.

0

3.>.).>

11111111 .44).J.4 154(11155154

11111111

0 ““ a).’.CU’ ~(IIl’44 lUll 0000 U’l’.I~0 In-SO)?”. ((10013.15’. U),-,00 I13U)NII’ ‘3113113(13 C.ON.r’S P-I’s’3r‘O$~U’, (14CIn4 ‘01-P-f’. ~

7 ,Lu14 3’ 4) 155

3.1st C C C C .,C

0 C.C • C C

CII

CU) ‘1” 01.., IA Or Cl CC ‘N. C • C II II C 1._C 4)0 Is) -IX ‘3 C P.C C • C) II .4 C.’. CC ,LI. 5...-’ • .3..

1-~-r0Il S13.’O’3I43 4113.04-

0’r’-,-’C

(1-1(40,U’,.’U’In .—U)OrU . • ‘0(41-1-

11111111

CII

C CU) 1CI’-. LII C’,.. - 0 I 00 ‘C 0 •

an

II

(15(14111(14 0004) lUll 0000 041(13. ,-4In1O ‘O’OflJI,OSCr..00 (143114$

‘OGInC 3.40’,(4.-SIll’.. r’-,’-f-l’G 4(1103) ,0~’-I-,Ø 1130000’ • . • . 41-0 I-Cl I 11111111

II

‘‘~,‘,

.4 11111111

11111111 040 Lu

0000 04Ls.Ls.ts.

C’ , . 540 7.0 0000 0000 (3000 0000’4-, COOP) ,-,-,—,-

11111111 .4.4.4.4 PSPIP’J(1.5

0000 (.11.717(7 II

II

C _a 7.3

3.4

C

5-” ,-000 444.3 114+4+4 C 0000 C ‘3003. C 00(13400 .1” 5(10.3.,-

0 • C I II .4 CI~, CC CU, • 0):.

2004 0000 0000 0000

0000 01.4.00 ‘40(40

(1~In0l14 InII3IIs’C 401430’ P-In00I13 CInOP-

C

I

11111111 .4.4,4.4 (‘.5(1.5(1,1154

C CC.) 0 C

U

C II

7 .-115) 1~ 4)1 (1-s

01-1-0 344-5 +44+ 0000 ~.a)r’,.. ~4U).-l’.

2’ ‘.,00 3.1 .4 (14

• CII

..(IlP-(4 P-CCU’

CII

C CII) ‘4C.’,

3.0.41,) CC-Sr.. (‘k’ CO. ‘34313.00 ‘(143(01-

.0 CII) ‘4014

51’. Cr

‘CCC)’. • . . . “rr II U

F.’. Cr C I CC ‘0 0 ‘ 114

CC • 0 0 • 0) II (I C (‘‘(0 55,1)5 154 5(1’ C 0 (40 0) • C ‘311 I .4 Cl’CC (CU. 0.-’ • .7.:

01-1,10, 0 0 .4 110 0 >1.0

114 0 C C .1-

.551340 C 00 C • C

I

II ‘3C’3C

0000 0414.14.114

(111111

II

3. 1-10. (“1(4 .4 155

OrVl-,) LI 0 .4 tIC 0 XLI

“5InU’U)U) 0000 lUll C.000 15(3.03. 0900”. 3.1,3,-C (1-U’ ,“Lnl13U’ (154113(1-5-504 P-COIL’. r0~’C (14113(1—5(1 ~4InC ~P-00

II 1N’3)C 000113 1111 COCCI (1-0414.CO’In C0’.OU’ 0.431(400 43000.0.

,—

0000 ‘.317(017

N.

0.5-00 CI U i IC 0 3(0

,~

0000 (41414(7

II Irrrp’l 00CC liii 0000 (1.-OIL’,— (1-In.-?” I’.InC’.( InN.’3433.311r-.-OlsJ

II II

I’~C4 IA.3. Lu sIX C.

11111111

II

0 ,.J 3.Lfl

III

C 1411111111111111 5 5 00000000 is.U. Ll.Ls.I~.

..— —. II

C “5” (0114(151,10 0000,11111 000CC 3.5(1904 0(1.0044 00CCG 03.0’ .1-3 (0N.N.04 .30.5(100.0 ,-.C1’141,I~03.15-GIn CI’-I~n3’3 .GIn044 .-p-CIn,-

II —

,‘ ,— II

II rCIflr0”CrU’ CCICCCOIIlI’.L 1+111111 00000000 41-F’I.OlIIO’O’U’ 0(1-1130.9 U)’0U)flJf’44~~ 0.0.In0’4-I’0l43C U’0’3.00C~,-C P-’OI43C.01’,43,Ir

II C 1.3. 5,4CC 54 sOS “5’ Cs CC’ C • C’ II

0 0 .4 110 0 ((U

OO’.C

—

U)-3-t~3nJ

,-,II

3-

.. ,J

C

II CC

(4

liii

U)U)U)C ~ r43’CU’

r,fl50(1CCCCLC 111+4 00000 InIn-a)0’ 0’U)O.O.’~? 41”1.3.3 11530’44 41C’O’ U)NOFII’.4 .U((1?’.I13l’.. 0”CC..OIfl 03.04-00’ 0(1.5,01(1 ,-‘~r13r- .3 (01-1,100

((1

11111111

1-001-4343411303 L3?’-000’0’3.In U’,’0rc’0113r.. UI P-CU) U-tl?’—N-U’ C1 CI11’lP,~,~IAr 5110) CC .0 C’ • II II II II 3’

0 0 0 ._4 HO 0 310

5-555-” U)4U”0,(1~U)(0 ~~r~(4(044 11111111 00000000 InU’O.-3OCU)(0 CInr-GU)InoOr’. (1OCOCN.9(’JN. 43)5(10 01(1(1—In ..t00C3.’S1,IO’ ,-O,-00,00,0 (‘5ON.Q’3..CF.I0’ 4r~OU’I,4-54(14113.000. 0’,O,-COl’3O”O ~r(1-31CU’ $~O(0CIn(0

If 1113. CC +1 CC 0-1”( (‘.00 13.0 1’4I’) I-I’) P-PS

3. U’ a) 3, C)

~CCC

,-a).tr140 o’CC-*-, t(140.U

II

C _J 3.15’.

I

•. .J J Lu

O~O(U

II Ii

S 00000000

U’In CC 1+ CC (1(1,1 (10 ON. rO. 0’, 1131(4 ,-O’ InC P-i’) 113,’ 1143.

‘4,U).3

,-.

3. U. 3.

~,SF,0 ,.•••..•. LflrNlflt.JO. (4(0 I I I I II II II II II (55 C C 3. ~“

~

,,-

(14a)O’.O

U~rcO’C,O.I’4~.3

II

II

‘C~ .‘-C (1.1, 434

1111

InInr~O.0aO’ ~ O~.l14 ~CI’4O.C’O 0’0I13r’.U)3)0 ,-rlI3,-In(S5r,4I(1

II

00 (317

CU) ‘ICI-,

Or 0 I CC .0 0 • In

C1’Jfld.4.O’Zr 000000L.In 4.+++++++ COO CC .OU)CU)44.. ~,-CC~J U~4..4U)U’C

((1

Ill

U) II) LA. >5

3. O — (1 0 a) 3. 0. 00. C

-Sa) ‘sC I CC 41’., 1140 -Sr..

II (11$ 0,— +1 CC ‘0,‘03. -SC 31.0 (‘153. ‘03.

-*,-r-,Dn3sn

C(,’-IF,O C-’0.S ?flr ICIflIC I I

II

U

II 3. S.C 3.04

(14?~t,’.O. .00) 04C..(U’

3. ~ C 3. (4

~ U’~ 00 I

‘CCC

fla)4.fl

,~njr...ojU)fl~.C’~

.. _I _J Lu

II

3.4

~

00’

=3.

01..) (.317 II

-~

III~

‘~‘~1

==

lIllI~++

CC 0.4 .31(1 4114 ‘34 p-a)

3. U. 3.

,-,-efljNCCCCC ,llI+

.~C,-.-

~4 ,~‘.3 ~‘~3

till

U) 3. 0 4 C C 0

~-r14,~-3O.r.. OU)O,I~.. U’~’41’.”I’G ~,U)Or~

CC

I

4) Lu

I’-r~J .-‘3

3.

rII) (1150 II 00

I’) U) 3. — I— C 0 IX

~CCC

114-S

,J

(Is

CC U’C

CC CO. .30

O.0

II

$irdt’4CN-21’3,..

C’. • -, ,S.’

_i

144 C C 3.

a),$ Cr14 II

11111111

II ~ 11111111 ‘.000 Is.U,Ls.U.

4.400 0000 0(400 00CC

0000 00040 0C7C~ 0,-lIla

CI 0 .4 IC .0 XC)

11111111 .4.4.4.4 r’sISIP’IVsJ

C 4-U’ (54(51. 04 3. C 0 PlC Ci ‘ C)

.4 C-’ CC CIA, 0’.’. (C:: 0 0 0 .4 ISO 0 3<14