### Program Code for Regular Coulomb Wave Function

```/*************************************************************************************
*   This program calculates the regular Coulomb wave function (the radial part of the
*   continuum wave function for a hydrogen-like atom).
*   It is written in PL/1,  but should be easily adaptable to other languages.
*   (note however that it requires extended numerical accuracy in order for the numerical
*   integration to work correctly (due to the oscillating integrand)).
*
*    It is based on a numerical calculation of the integral representation
*    Flη(ρ) = e-πη/2.ρl+1/2l/|Γ(l+1+iη)|.0∫∞dt (1-tanh2(t))l+1.cos(ρ.tanh(t)-2ηt)
*
*    and the recursive formula
*    Fl-1η(ρ) = [ (2l+1).(η +l.(l+1)/ρ). Flη(ρ) - l.√((l+1)2+η2). Fl+1η(ρ) ] /
*                                                          /[(l+1).√(l2+η2) ]
*
*  (see Abramowitz and Stegun, Handbook of Mathematical Functions, for both relationships)
*
*  This program can work in two output modes: either the results are printed (at equidistant
*  radius points), or they are stored on file at the quadrature point of a 16-point Gauss
*
*
*  Input parameters are:
*    LEN:  lowest energy for which the wave functions should be computed (in Rydberg)
*    UEN: highest energy for which the wave functions should be computed (in Rydberg)
*    DEN: energy step width (in Rydberg)
*    Z1:  electron charge number (<0(=-1))
*    Z2:  nucleus charge number (>0)
*    L0:  maximum angular momentum
*    NMIN:  lowest bound state to be used later for the calculation of the dipole integrals
*    KV:  desired relative accuray of the results (≥10-12-10-7 depending on UEN and L0)
*    NDIG:  number of significant digits after the decimal point for output
*    XD00:  integration stepwidth for Gauss Quadrature for l=0 (recommended =0.5)
*    KVG:  desired relative accuracy for Gamma function (should be <=KV)
*    XDG:  integration stepwidth for Gamma function (recommended =5)
*    OUTP  :   =1: print results;    =2 store results on file
*    LR: (only relevant for OUTP=1) smallest radius point (in Bohr radii)
*    UR: (only relevant for OUTP=1) largest radius point (in Bohr radii)
*    DR: (only relevant for OUTP=1) radius point step width (in Bohr radii)
*    KVV: (only relevant for OUTP=2) desired relative accuracy of the dipole integrals
to be computed later from the results (determines the maximum radius point here)
*    XDS0: (only relevant for OUTP=2) determines the density of radius points
*                                        for small energies (recommended =5)
*    PRSTOP:  (only relevant for OUTP=2) if=1: wave function is only computed up to radius URMAX
*                   (otherwise the maximum radius is automatically determined by the routine RSMAX)
*    URMAX:  (only relevant for OUTP=2 and PRSTOP=1) maximum radius for which wave function
*                                        should be computed (in Bohr radii):
*
*    Required subroutines: GAMMA (returns Gamma function),
*                          RSMAX (returns maximum radius point
*                                to be used for dipole integral calculations).
*
*
*   This program has worked well and correct for me in the past (using IBM PL/1 compilers),
*   but I can't guarantee this for all compilers, especially if the code is translated
*   into a different language.
*   Note also that this program is quite CPU intensive and may not be the fastest way
*   to calculate the Coulomb wave functions (especially for small energies and
*   high angular momenta, for which my alternative program COULOMBC2 is the better
*   option).
*
*   If there are any problems, please let me know and I mention the issue here.
*   (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
*
*                                               Thomas Smid, April 2007
*********************************************************************************/

COULOMBC: PROCEDURE OPTIONS(MAIN) REORDER ;

DEFAULT RANGE(A:Z) DEC FLOAT VALUE (DEC FLOAT(33)) ;
DCL QG(16) INIT
(9.894009349916499325961541735E-01, 2.71524594117540948517805725E-02,
9.445750230732325760779884155E-01, 6.22535239386478928628438370E-02,
8.656312023878317438804678977E-01, 9.51585116824927848099251076E-02,
7.554044083550030338951011948E-01, 1.246289712555338720524762822E-01,
6.178762444026437484466717640E-01, 1.495959888165767320815017305E-01,
4.580167776572273863424194430E-01, 1.691565193950025381893120790E-01,
2.816035507792589132304605015E-01, 1.826034150449235888667636680E-01,
9.50125098376374401853193354E-02, 1.894506104550684962853967232E-01)
REAL DEC FLOAT(33) ;
DCL QG2(16) DEC FLOAT(16) ;
DCL CPOINT INIT(2.886751345948128E-01) DEC FLOAT(16) ;
DCL(I,IR,IL,J,IR2,IRA,CO,OUTP,PRSTOP,IE,JE,DI,NDIG) DEC FIXED(7,0);
DCL PI INIT(3.141592653589793E00) DEC FLOAT(33) ;
DCL (LE,GLE) COMPLEX DEC FLOAT(33) ;
DCL (KVG,XDG) REAL DEC FLOAT(33) ;
DCL (RPOINT(*),FLN(*),FLL(*),FLU(*)) CTL DEC FLOAT(16) ;
DCL (LR,UR,DR,LEN,UEN,DEN) DEC FIXED(12,6) ;
DCL OUTFILE FILE STREAM OUTPUT ;
DCL (A,B,C,XU,XL,XD,LY,QY,CHX,CHX2,RO,ETA2,L1,FL) DEC FLOAT(33) ;
DCL GAMMA ENTRY (REAL FLOAT(33),REAL FLOAT(33),CPLX FLOAT(33),
CPLX FLOAT(33)) RETURNS (CPLX FLOAT(33)) ;
DCL RSMAX ENTRY ;

ON UNDERFLOW ;

PI20=PI/2 ;  SQPI2=SQRT(PI20) ;
QG=QG/2 ;
QG2=QG ;

GET DATA COPY ;      /* Read input data */

ULIM0=1E28*KV ;
IE=FLOOR((UEN-LEN)/DEN)+1 ;
NMAX=L0 ;
IF OUTP=2 THEN UR=0 ;

DO JE=1 TO IE ;                           /* Energy loop */
E=(JE-1)*DEN+LEN ;
K=SQRT(E) ;
ETA=Z1*Z2/K ;
ETA2=2*ETA;  ETA2A=ABS(ETA2) ;

IF OUTP=2 THEN DO ;
LR=0; URS=0;
DR0=MIN(1.5*PI/K,XDS0*NMIN/Z2) ;  DR=DR0 ;  DRR=DR ;
IF PRSTOP=1 THEN UR=URMAX ;
ELSE
CALL RSMAX(Z2,L0,E,DRR,KVV,URS) ;  /* Determine maximum radius point */
UR=MAX(UR,URS) ;
UR=CEIL(UR/DR)*DR ;
END ;

IF OUTP=2 THEN IR=FLOOR((UR-LR)/DR)+7 ; ELSE IR=FLOOR((UR-LR)/DR)+1;
IR2=(IR-1)*16 ;
IF OUTP=2 THEN IRA=IR2 ;  ELSE IRA=IR ;
ALLOCATE FLL(IRA),FLU(IRA),FLN(IRA) ;
PUT DATA(IE,IR,IR2,IRA) ;  PUT SKIP ;

IF OUTP=2 THEN DO ;

ALLOCATE RPOINT(IRA) ;

DO J=1 TO IR-1 ;      /* Determination of radius points */
IF J<=9 THEN DRR=DR0/3; ELSE DRR=DR0;
IF J<=9 THEN R=(J-1)*DRR+LR ; ELSE R=3*DR0+(J-10)*DRR ;
DO I=1 TO 8 ;
RPOINT(16*(J-1)+I)=R+0.5*DRR-QG(2*I-1)*DRR ;
RPOINT(16*(J-1)+17-I)=R+0.5*DRR+QG(2*I-1)*DRR ;
END ;
END ;
PUT FILE(OUTFILE) EDIT(KV,E,Z1,Z2,L0,LR,UR,DR0,QG2)
(SKIP,24 E(NDIG+7,NDIG)) ;

END ;

FLL=0; FLU=0; FLN=0; FLN1=0; FLN2=0;

DO IL=1 TO 2 ;       /* Loop for two highest angular momenta */

IF IL=1 THEN L=L0 ;   ELSE L=L0-1 ;
L1=L+1 ;
LE=CPLX(L1,ETA) ;
CALL GAMMA(KVG,XDG,LE,GLE) ;
GAE=EXP(-PI20*ETA)/2**L/SQRT(GLE*CONJG(GLE))/SQPI2 ;
XDL=XD00/SQRT(L1) ;

DO J=1 TO IRA ;   /* Radius points loop */

IF OUTP=2 THEN R=RPOINT(J) ;
ELSE R=(J-1)*DR+LR ;
RO=R*K    ;
FL=0; XU=0;
CONVGC=1; JJ=1; QY=0; DI=1; PI2=PI20*2; XD0=XDL;
QYU=0 ;

INT: DO WHILE(CONVGC>=KV) ;  /* Numerical integration loop */
QYL=QY ;
XL=XU ;
XD=MIN(XD0,PI2/(RO/COSH(XL)**2+ETA2A)) ;
IF XL=0 THEN XD=XD/2 ;
XU=XL+XD ;
CALL DQG16 ;
QYU=QYU+ABS(QY) ;
FL=FL+QY ;
QYM=ABS(FL)/JJ  ;
IF QYM>0 THEN CONVGC=(ABS(QY)+ABS(QYL))/2/QYM ;
JJ=JJ+1 ;
END ;

IF DI=1 THEN DO ;  /* Repeat last integration step with half stepwidth */
DI=-1; CONVGC=1;
FL=FL-QY; JJ=JJ-1; XU=XU-XD; QYU=QYU-ABS(QY) ;
PI2=PI20; XD0=XDL/2;
GO TO INT ;
END ;

CFL=GAE*RO**L1 ;
FL=FL*CFL ;
IF (QYU/2*CFL>ULIM0)|(QYU/2*CFL/FL>1E28) THEN CALL ERPRINT ;
IF IL=1 THEN FLU(J)=FL ;
ELSE FLL(J)=FL        ;

END ;   /* End  'Radius points loop' */

END ;  /* End  'Loop for two highest angular momenta' */

FLN=FLU; L=L0;
CALL OUTPUT ;
FLN=FLL; L=L-1;
CALL OUTPUT ;

FLL=FLU ;
DO IL=L0-1 TO 1 BY -1 ;  /* Recursive calculation of other angular momenta */
L=IL ;
L1=L+1 ;
FLU=FLL; FLL=FLN;
DO J=1 TO IRA ;
IF OUTP=2 THEN R=RPOINT(J) ;
ELSE R=(J-1)*DR+LR ;
RO=R*K ;
FLN1=(2*L+1)*(ETA+L*L1/RO)*FLL(J) ;
FLN2=L*SQRT(L1**2+ETA**2)*FLU(J);
IF ((ABS(FLN1)+ABS(FLN2))/ABS(FLN1-FLN2))>ULIM0 THEN CALL ERPRINT ;
FLN(J)=(FLN1-FLN2)/L1/SQRT(L**2+ETA**2) ;
END ;
L=L-1 ;
CALL OUTPUT ;
END ;

FREE FLL,FLU,FLN ;
IF OUTP=2 THEN DO ;
FREE RPOINT ;
END ;

END ;    /* End 'Energy loop' */

/* Internal Subroutines */

OUTPUT: PROCEDURE ;
IF OUTP=1 THEN DO ;
PUT EDIT(L) (SKIP(5),E(NDIG+7,NDIG),SKIP) ;
END ;
DO J=1 TO IRA ;
IF OUTP=1 THEN DO ;
R=(J-1)*DR+LR ;
PUT EDIT(R,FLN(J)) (SKIP,2 E(NDIG+7,NDIG)) ;
END ;
ELSE DO ;
PUT FILE(OUTFILE) EDIT(FLN(J)) (SKIP,E(NDIG+7,NDIG)) ;
END ;
END ;
END OUTPUT ;

ERPRINT:PROCEDURE ;
PUT EDIT('DESIRED ACCURACY CAN NOT BE REACHED') (A(35)) ;
PUT SKIP; PUT DATA(K,R,L1,CFL,QYU,FL,FLN1,FLN2); PUT SKIP;
GO TO ENDE ;
END ERPRINT ;

FCT:PROCEDURE(X) ;
DCL X REAL    DEC FLOAT(33) ;
CHX=COSH(X) ;CHX2=CHX**2 ;
RETURN(COS(RO*SQRT(CHX2-1)/CHX-ETA2*X)/CHX2**L1)     ;
END FCT ;

DQG16: PROCEDURE ;
A=0.5*(XU+XL) ;
B=XU-XL ;
LY=0 ;
DO I=1 TO 15 BY 2 ;
C=QG(I)*B ;
LY=LY+QG(I+1)*(FCT(A+C)+FCT(A-C)) ;
END ;
QY=LY*B ;
END DQG16 ;

ENDE: END COULOMBC ;

/********************************************************************************
*  The subroutine for the calculation of the Gamma function uses a direct
*  numerical integration (adaptive stepwise 16-point Gauss quadrature) of Euler's integral
*
*   Γ(z)=0∫∞dt tz-1.e-t
*
*
*  Input parameters are:
*    KV: desired relative accuracy
*    XD: integration stepwidth (recommended =5)
*    Z: argument of gamma function (complex)
*
*  Output parameter is:
*    G: value of Gamma function (complex)
*
*   Note that this program uses extended accuracy because of the possibility
*   that the integrand oscillates rapidly and the integral becomes very small
*   as a result. Normal (double or single) accuracy may not produce correct
*   results in these cases.
*   Note also that this may not be the computationally most efficient method
*   to calculate the Gamma function. Other methods may be much faster, but
*   are not as flexible with regard to the desired accuracy for instance.
*****************************************************************************/

GAMMA: PROCEDURE (KV,XD,Z,G) ;
DEFAULT RANGE(A:Z) COMPLEX FLOAT VALUE (DEC FLOAT(33)) ;
DCL QG(16) INIT
(9.894009349916499325961541735E-01, 2.71524594117540948517805725E-02,
9.445750230732325760779884155E-01, 6.22535239386478928628438370E-02,
8.656312023878317438804678977E-01, 9.51585116824927848099251076E-02,
7.554044083550030338951011948E-01, 1.246289712555338720524762822E-01,
6.178762444026437484466717640E-01, 1.495959888165767320815017305E-01,
4.580167776572273863424194430E-01, 1.691565193950025381893120790E-01,
2.816035507792589132304605015E-01, 1.826034150449235888667636680E-01,
9.50125098376374401853193354E-02, 1.894506104550684962853967232E-01)
REAL DEC FLOAT(33) ;
DCL I REAL DEC FIXED(5,0) ;
DCL (Z,G) COMPLEX DEC FLOAT(33) ;
DCL (CONVGC,KV,RZ, JJ,XL,XU,XD,EIZ,PI2 INIT(3.141592653589793E00),IZ,
QYMR,QYMI,GR,GI,CONVGCI,CONVGCR,QYR,QYI,DI,XL0) REAL DEC FLOAT(33) ;

ON UNDERFLOW ;
IZ=ABS(IMAG(Z)) ; RZ=REAL(Z) ;
IF IZ>0 THEN EIZ=EXP(PI2/IZ); ELSE EIZ=1E10/KV ;
QG=QG/2 ;

XL0=MAX(RZ-1,1) ;  XU=XL0 ;
XD=XD/(FLOOR(RZ/32)+1) ;

G=0  ;
CONVGC=1 ; JJ=1 ;DI=1;

INT:
DO WHILE(CONVGC>=KV) ;    /* Integration loop */
IF DI=1 THEN DO ;   /* Integration into positive direction from real maximum  */
XL=XU ;
XU=XL+MIN(XL*EIZ-XL,XD) ;
END ;
ELSE DO ;   /* Integration into negative direction from real maximum  */
XU=XL ;
XL=MAX(XU/EIZ,XU-XD) ;
END ;
CALL DQG16 ;
G=G+QY ;
GR=REAL(G) ;  QYR=REAL(QY) ;
GI=IMAG(G) ;  QYI=IMAG(QY) ;
QYMR=ABS(GR)/JJ ;
QYMI=ABS(GI)/JJ ;
IF QYMR>0 THEN CONVGCR=ABS(QYR/QYMR);
IF QYMI>0 THEN CONVGCI=ABS(QYI/QYMI);
CONVGC=MAX(CONVGCR,CONVGCI) ;
JJ=JJ+1 ;
END ;

IF DI=1 THEN DO ;
CONVGC=1 ; DI=-1 ; XL=XL0 ;
GO TO INT ;
END ;

/* Internal Subroutines */

FCT:PROCEDURE(X) ;
DCL X REAL DEC FLOAT(33) ;
RETURN(X**(Z-1)*EXP(-X))   ;
END FCT ;

DQG16: PROCEDURE ;
A=0.5*(XU+XL) ;
B=XU-XL ;
LY=0 ;
DO I=1 TO 15 BY 2 ;
C=QG(I)*B ;
LY=LY+QG(I+1)*(FCT(A+C)+FCT(A-C)) ;
END ;
QY=LY*B ;
END DQG16 ;

ENDE: END GAMMA ;

/*******************************************************************************
*  The subroutine RSMAX returns the maximum radius point for which the continuous
*  Coulomb wave function needs to be calculated in COULOMBC in order to achieve
*  convergence of the corresponding bound-free dipole integral. Essentially, it is a
*  stepwise approximate calculation of the dipole integral using a sine-function for
*  the continuous wave function and assuming the highest angular monentum for the
*  bound state.
*
*  Input parameters are:
*    Z: nuclear charge number
*    N0: principal quantum number of the highest bound state in the dipole intergral
*    E: continuum energy (in Rydberg)
*    DR: radial stepwidth for 16-point Gauss quadrature
*    KV: desired relative accuracy of result
*
*  Output parameter is:
*    RMAX: maximum radius point (in Bohr radii)
*****************************************************************************/

RSMAX: PROCEDURE(Z,N0,E,DR,KV,RMAX) ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(33)) ;
DCL QG(16) INIT
(9.894009349916499325961541735E-01, 2.71524594117540948517805725E-02,
9.445750230732325760779884155E-01, 6.22535239386478928628438370E-02,
8.656312023878317438804678977E-01, 9.51585116824927848099251076E-02,
7.554044083550030338951011948E-01, 1.246289712555338720524762822E-01,
6.178762444026437484466717640E-01, 1.495959888165767320815017305E-01,
4.580167776572273863424194430E-01, 1.691565193950025381893120790E-01,
2.816035507792589132304605015E-01, 1.826034150449235888667636680E-01,
9.50125098376374401853193354E-02, 1.894506104550684962853967232E-01)
REAL DEC FLOAT(33) ;
DCL (I,DI) DEC FIXED(7,0) ;
DCL PI INIT (3.141592653589793E00) DEC FLOAT(16) ;
ON UNDERFLOW ;
QG=QG/2 ;
K=SQRT(E) ;
K2=2*K ;
ZN=Z/N0 ;

CONVGC=1; DI=1; JJ=1; INT=0; XU=0;
LOOP: DO WHILE(CONVGC>=KV) ;      /* Integration loop */
XL=XU ;
XU=XL+DR ;
CALL DQG16 ;
INT=INT+QY ;
QYM=ABS(INT)/JJ ;
IF QYM>0 THEN CONVGC=ABS(QY/QYM) ;
JJ=JJ+1 ;

END ;

IF DI=1 THEN DO ;  /* Additional test step after convergence is achieved */
DI=-1; CONVGC=1;
GO TO LOOP ;
END ;

RMAX=XU ;

/* Internal Subroutines */

FCT: PROCEDURE(X) ;
DCL X REAL DEC FLOAT(33) ;
RETURN(X**N0*EXP(-ZN*X)*SIN(K*X)) ;
END FCT ;

DQG16: PROCEDURE ;
A=0.5*(XU+XL) ;
B=XU-XL ;
LY=0 ;
DO I=1 TO 15 BY 2 ;
C=QG(I)*B ;
LY=LY+QG(I+1)*(FCT(A+C)+FCT(A-C)) ;
END ;
QY=LY*B ;
END DQG16 ;

ENDE : END RSMAX  ;

```