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
* quadrature (over automatically adapted (equidistant) radial intervals).
*
*
* 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 ;
|