Program Code for Coulomb Wave Function
(specially adapted for very small energies)
/*************************************************************************************
* 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.
* It is based on a program that was part of the PL/1 scientific subroutine
* package at the University of Bonn, which I modified so that it can be used
* for very small energies and high angular momenta (the modification
* consists of using special subroutines for arithmetic operations that
* could otherwise not be completed within the normal numerical range).
* This routine should in principle be faster than my routine COULOMBC
* as it does not involve numerical integration but uses a recursion technique,
* and also it requires only the usual double precision.
*
* The program calculates the wave function for a given radial interval
* and stores it on file. Subsequent runs can append to this file if required.
* There is in principle no limit to how far the radial range can be extended
* this way, so for the calculation of the dipole integral with the wave
* function of a given bound state, one should know beforehand what about the
* maximum radius required is (for instance, by means of my subroutine RSMAX).
*
*
* Input parameters are:
* CO: if CO=0: output file is newly created in this run;
otherwise: output is appended to existing file
* MOD: if ='TEST': print out of recursion index
* Z1: electron charge number (<0(=-1))
* Z2: nucleus charge number (>0)
* E: energy of continuum electron (Rydberg)
* N: highest angular momentum to be computed
* DL: angular momentum step width for output
* LMAX: highest angular momentum for output
* URM: (only relevant for CO=0): maximum radius point for which the wave
* function will be stored (including subsequent extensions of the file)
* (in Bohr radii)
* LR: (only relevant if LRMOD not ='FILE') smallest radius point for this run
* (in Bohr radii)
* UR: (only relevant if LRMOD not ='FILE') largest radius point for this run
* (in Bohr radii; <URM)
* LRMOD: if ='FILE': LR is taken from external file
* DRF: (only relevant if LRMOD ='FILE'): radius step width for this run;
* if CO not=0 : if this run terminates regularly, UR=LR+DRF is stored;
if not, a negative value is stored (preventing further runs);
* if CO=0: >4.5.π.DRFAK.√E
* (in Bohr radii)
* DRFAK: multiplicator for changing the density of radius points
* D: desired relative accuracy of results
* LIMFAK: multiplicator for length of recursion array vectors (=1.5)
* NUMAX: if >0: internally computed maximum recursion index will be
* replaced by this value
* DNU: (only relevant if NUMAX>0): step width for increase of recursion index
* FNDIG: number of relevant digits after the decimal point for output
*
*
* Required subroutines: NUM, ADD, MULT, DIV
* (replacing basic arithmetic operations where these could lead
* to a numerical overflow or underflow).
*
*
* 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.
*
* 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, May 2007
*********************************************************************************/
COULOMBC2: PROCEDURE OPTIONS(MAIN) REORDER ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
DCL (ETA,RHO,EPS,D1,R,S,OMEGA,RHO1,A,B,SUM,PI2) DEC FLOAT(16) ;
DCL (F(*),RR(*),FAP(*),LAM(*)) CTL DEC FLOAT(16) ;
DCL (LAMEX(*),FEX(*),RREX(*),FAPEX(*)) CTL DEC FIXED(15,0);
DCL (N,D,LIM,CTR,K,NU,DNU,K1,NU1,DL,LMAX) DEC FIXED(15,0) ;
DCL (XEX,YEX,XYEX,REX,SEX,RLEX,SLEX,S1EX,NULL,SUMEX) DEC FIXED(15,0);
DCL(AEX,R1EX,R2EX,SQBEX,KKEX,BEX,NUMAX,CO) DEC FIXED(15,0) ;
DCL MMEX DEC FIXED(15,0) ;
DCL RPOINT(*) CTL DEC FLOAT(16) ;
DCL QG(16) INIT
(4.947004674958250E-01, 1.357622970587705E-02,
4.722875115366163E-01, 3.112676196932395E-02,
4.328156011939159E-01, 4.757925584124639E-02,
3.777022041775015E-01, 6.231448562776694E-02,
3.089381222013219E-01, 7.479799440828837E-02,
2.290083888286137E-01, 8.457825969750127E-02,
1.408017753896295E-01, 9.130170752246179E-02,
4.750625491881872E-02, 9.472530522753425E-02);
DCL QGF(16) DEC FLOAT(16) ;
DCL (MOD,LRMOD) CHARACTER(4) ;
DCL (LR,UR,DR,LR1,UR1,URM,DRF,URF,LR0) DEC FIXED (12,3) ;
DCL (IR,IRA,J,JJ,I,I2,JJ1,FNDIG) DEC FIXED(15,0) ;
DCL OUTFILE FILE STREAM OUTPUT ;
DCL ILRFILE FILE STREAM INPUT ;
DCL OLRFILE FILE STREAM OUTPUT ;
DCL (NUM,ADD,MULT,DIV) ENTRY ;
ON UNDERFLOW ;
GET DATA COPY ; /* Data input */
ALLOCATE F(0:N),RR(0:N),FAP(0:N),FEX(0:N),
RREX(0:N),FAPEX(0:N) ;
KE=SQRT(E) ;
ETA=Z1*Z2/KE ;
OPEN FILE(OUTFILE) ;
MEINS=-1 ;
PI=ACOS(MEINS) ;
DR=1.5*PI/KE*DRFAK; DRR=DR;
IF LRMOD='FILE' THEN DO ;
IF CO>0|CO<0 THEN DO ;
OPEN FILE(ILRFILE); GET FILE(ILRFILE) EDIT(LR) (SKIP,E(22,15)) ;
IF LR<0 THEN GO TO ENDE ;
END ;
ELSE LR=0 ;
LR0=-1E07 ;
OPEN FILE(OLRFILE) ;
PUT FILE(OLRFILE) EDIT(LR0) (SKIP,E(22,15)) ;
CLOSE FILE(OLRFILE) ;
URF=LR+DRF; UR=URF ;
END ;
LR=CEIL(LR/DR)*DR ;
UR=CEIL(UR/DR)*DR ;UR1=URM ;
UR1=CEIL(UR1/DR)*DR ;LR1=0 ;
IF CO=0 THEN IR=FLOOR((UR-LR)/DR)+7 ;
ELSE IR=FLOOR((UR-LR)/DR)+1 ;
IRA=(IR-1)*16 ;
ALLOCATE RPOINT(IRA) ;
KVV=1E-04; IF LMAX>=N THEN LMAX=N-DL ;
DO J=1 TO IR-1 ; /* Calculation of radius points */
IF J<=9&CO=0 THEN DO ;
DRJ=DRR/3 ;
R=(J-1)*DRJ+LR ;
END ;
ELSE DO ;
DRJ=DRR ;
IF CO=0 THEN R=3*DRR+(J-10)*DRJ ;
ELSE R=LR+(J-1)*DRJ ;
END ;
DR2=0.5*DRJ ;
JJ=16*(J-1) ; JJ1=JJ+17 ;
DO I=1 TO 8 ;
I2=2*I-1 ;
RPOINT(JJ+I)=R+DR2-QG(I2)*DRJ;
RPOINT(JJ1-I)=R+DR2+QG(I2)*DRJ;
END ;
END ;
IF CO=0 THEN DO ; /* Output of file header */
PUT FILE(OUTFILE) EDIT(KVV,E,Z1,Z2,N,LMAX,DL,LR1,UR1,DR,QG)
(SKIP,26 E(FNDIG+7,FNDIG)) ;
END ;
DO J=1 TO IRA ; /* Radius point loop */
RHO=RPOINT(J)*KE ;
F=0; RR=0; FAP=0; RLEX=0; SLEX=0; S1EX=0; NULL=0;
FEX=0; RREX=0; FAPEX=0; AEX=0; R1EX=0; R2EX=0; SQBEX=0; KKEX=0;
BEX=0;
IF RHO<0|N<0 THEN GO TO ENDE ;
IF RHO=0 THEN GO TO C3 ;
IF D<2 THEN D=2; DEPS=D;
EPS=0.5/10**DEPS ;
RHO1=1/RHO ;
PI2=1.570796326794897E00 ; SQPI2=SQRT(1/PI2) ;
LRR0=1E-14; REX=0; SEX=0; EINS=1 ;
DO K=0 TO N ;
FAP(K)=0 ;
END ;
R=2*ETA*RHO1 ;
IF R<0 THEN OMEGA=0 ;
ELSE DO ;
OMEGA=PI2*R ;
IF R>1 THEN DO ;
S=ATAN(SQRT(SQRT(R*R-1))) ;
OMEGA=OMEGA-S*R+SQRT(R-1) ;
END ;
END ;
D1=2.3026E00*D+1.3863E00 ;
IF OMEGA=0 THEN DO ;
A=PI2*ABS(ETA) ;
END ;
ELSE A=EXP(-ATAN(1/OMEGA)*ETA) ;
S=SQRT(1+OMEGA*OMEGA) ; /* Calculation of maximum recursion index */
B=OMEGA+S ;
S=D1-OMEGA*RHO+A+LOG(SQRT(B/S)) ;
B=1.3591E00*RHO*B ;
S=S/B ;
IF S<-0.3679E00 THEN S=1 ;
ELSE DO ;
TS=S ;
CALL T ;
S=B*TT ;
END ;
IF N=0 THEN R=1 ;
ELSE DO ;
TS=0.5*D1/N ;
CALL T ;
R=N*TT ;
END ;
IF R<=S THEN R=S ;
NU=FLOOR(R)+2 ;
IF NUMAX>0 THEN NU=NUMAX ; ELSE NU=CEIL(NU*8E-01) ;
NU1=1 ; IF NUMAX=0 THEN DNU=CEIL(NU*1E-01) ;
LIM=CEIL(NU*LIMFAK); ALLOCATE LAM(0:LIM),LAMEX(0:LIM) ;
LAM=0; LAMEX=0; LAM(0)=1; LAM(1)=OMEGA-ETA ;
CALL NUM(LAM(1),LAMEX(1)) ;
SUM=RHO*EXP(OMEGA*RHO) ;
CALL NUM(SUM,SUMEX) ;
C1: IF NU>LIM THEN DO ; /* Iteration loop */
PUT DATA(NU); PUT SKIP;
GO TO ENDE ;
END ;
IF MOD='TEST' THEN DO ;
PUT DATA(NU); PUT SKIP; /* Print recursion index in test mode */
END ;
DO K=NU1 TO NU-1 ; /* Recursive calculation of recursion coefficient */
A=1+1E0/K ;
B=ETA/K ;
R=OMEGA*(1+1/A) ;
IF ABS(R)>0 THEN CALL NUM(R,REX) ;
S=(1+B*B)/A ;
CALL NUM(S,SEX) ;
CALL MULT(R,REX,LAM(K),LAMEX(K),RL,RLEX) ;
CALL MULT(S,SEX,LAM(K-1),LAMEX(K-1),SL,SLEX) ;
CALL ADD(RL,RLEX,SL,SLEX,LAM(K+1),LAMEX(K+1)) ;
END ;
R=0; REX=0; S=0; SEX=0;
DO K=NU TO 1 BY -1 ; /* Recursive calculation of reduced wave function */
A=ETA/(K+1) ;
B=A/K+RHO1 ;
A=(1+A*A)/(2*K+3) ;
CALL NUM(A,AEX) ;
CALL NUM(B,BEX) ;
MM=R; MMEX=REX ;
CALL MULT(A,AEX,MM,MMEX,R,REX) ;
MM=-R; MMEX=REX;
CALL ADD(B,BEX,MM,MMEX,R,REX) ;
KK=2*K-1;
CALL NUM(KK,KKEX) ;
MM=R; MMEX=REX;
CALL MULT(KK,KKEX,MM,MMEX,R,REX);
MM=R; MMEX=REX;
CALL DIV(EINS,NULL,MM,MMEX,R,REX) ;
CALL MULT(R,REX,LAM(K),LAMEX(K),RL,RLEX) ;
CALL MULT(R,REX,S,SEX,SL,SLEX) ;
CALL ADD(RL,RLEX,SL,SLEX,S,SEX) ;
IF K<=N THEN DO ;
RR(K-1)=R; RREX(K-1)=REX ;
END ;
END ;
CALL ADD(EINS,NULL,S,SEX,S1,S1EX) ;
CALL DIV(SUM,SUMEX,S1,S1EX,F(0),FEX(0)) ;
IF N>0 THEN DO ;
DO K=0 TO N-1 ;
MM=F(K); MMEX=FEX(K);
CALL MULT(RR(K),RREX(K),MM,MMEX,F(K+1),FEX(K+1)) ;
END ;
END ;
DO K=0 TO N ; /* Check for convergence and repeat recursion if required */
IF ABS(F(K)-FAP(K))<EPS*ABS(F(K))&FEX(K)=FAPEX(K) THEN GO TO CC1;
ELSE DO ;
DO K1=0 TO N ;
FAP(K1)=F(K1) ; FAPEX(K1)=FEX(K1) ;
END ;
NU1=NU; NU=NU+DNU ;
GO TO C1 ;
END ;
CC1: END ; /* End iteration loop */
A=4*PI2*ETA ; /* Calculation of wave function from reduced wave function */
EPS=EPS/10 ; D2=D+1 ;
IF ABS(A)<0.75 THEN DO ;
CALL NUM(A,AEX) ;
SUM=1; SUMEX=0; K=1; R=1; REX=0;
C2: K=K+1; KK=K ;
CALL NUM(KK,KKEX) ;
MM=R; MMEX=REX ;
CALL MULT(MM,MMEX,A,AEX,R,REX) ;
MM=R; MMEX=REX;
CALL DIV(MM,MMEX,KK,KKEX,R,REX) ;
MM=SUM; MMEX=SUMEX;
CALL ADD(MM,MMEX,R,REX,SUM,SUMEX) ;
IF REX>SUMEX-D2 THEN GO TO C2;
SUME=SUMEX ;
SUM=SUM*10**SUME ;
END ;
ELSE SUM=(EXP(A)-1)/A ;
A=SQRT(1/SUM) ;
CALL NUM(A,AEX) ;
MM=F(0); MMEX=FEX(0) ;
CALL MULT(A,AEX,MM,MMEX,F(0),FEX(0)) ;
IF N>0 THEN DO ;
DO K=1 TO N ;
R=1E0/K ;
B=ETA*R ;
R1=2-R; R2=2+R;
SQB=SQRT(1+B*B);
CALL NUM(SQB,SQBEX) ;
MM=A; MMEX=AEX;
CALL MULT(MM,MMEX,R1,R1EX,A,AEX);
MM=A; MMEX=AEX;
CALL MULT(MM,MMEX,SQB,SQBEX,A,AEX);
MM=A; MMEX=AEX;
CALL DIV(MM,MMEX,R2,R2EX,A,AEX);
MM=F(K); MMEX=FEX(K);
CALL MULT(A,AEX,MM,MMEX,F(K),FEX(K)) ;
END ;
END ;
GO TO C4 ;
C3: DO K=0 TO N ;
F(K)=0 ;
END ;
C4: DO K=0 TO LMAX+1 ;
FEXK=FEX(K) ;
F(K)=SQPI2*F(K)*10**FEXK ;
END ;
DO K=0 TO LMAX BY DL ; /* Store results for desired angular momenta */
IF K=0 THEN FK1=0; ELSE FK1=F(K-1) ;
PUT FILE(OUTFILE) EDIT(FK1,F(K+1)) (SKIP,2 E(FNDIG+7,FNDIG)) ;
END ;
FREE LAM,LAMEX ;
END ; /* End radius point loop */
IF LRMOD='FILE' THEN DO ; /* Store maximum radius point for next run */
OPEN FILE(OLRFILE); PUT FILE(OLRFILE) EDIT(URF) (SKIP,E(22,15)) ;
END ;
/* Internal Subroutines */
T: PROCEDURE ; /* For calculation of maximum recursion index */
IF TS<=0 THEN TT=0.36788+1.0422*SQRT(TS+0.36788) ;
ELSE DO ;
IF TS<=10 THEN DO ;
TP=(0.1*TS-0.00176148)*TS+0.0208645 ;
TP=(TP*TS-0.129013)*TS+0.85777 ;
TT=TP*TS+1.0125 ;
END ;
ELSE DO ;
TZ=LOG(TS)-0.775 ;
TP=(0.775-LOG(TZ))/(1+TZ) ;
TT=TS/(TZ*(1+TP)) ;
END ;
END ;
END T ;
ENDE: END COULOMBC2 ;
/***********************************************************************************
* The subroutines NUM, ADD, MULT, DIV enable basic arithmetic operations of numbers
* with exponents up to ±1015.
* NUM separates usual floating point numbers (X) into a mantissa (X) and exponent (XEX),
* and ADD, MULT, DIV provide the basic arithmetic operations of these representations,
* where XM,XEX,YM,YEX are the two input numbers and XYM, XYEX the result.
***********************************************************************************/
NUM: PROCEDURE(X,XEX) ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
DCL X DEC FLOAT(16) ;
DCL XEX DEC FIXED(15,0) ;
LAX=LOG10(ABS(X)) ;
XEX=FLOOR(LAX) ;
X=10**(LAX-XEX)*SIGN(X) ;
END NUM ;
ADD: PROCEDURE(XM,XEX,YM,YEX,XYM,XYEX) ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
DCL (XM,YM,XYM) DEC FLOAT(16) ;
DCL (XEX,YEX,XYEX) DEC FIXED(15,0) ;
IF XM=0&YM=0 THEN DO;
XYM=0; XYEX=0;
RETURN;
END ;
IF XM=0 THEN DO;
XYM=YM; XYEX=YEX;
RETURN;
END ;
IF YM=0 THEN DO;
XYM=XM; XYEX=XEX;
RETURN;
END ;
IF XEX>=YEX THEN DO ;
XYEX=XEX ;
XYM=YEX-XEX; XYME=XYM;
IF XYM<-60 THEN XYM=0; ELSE XYM=10**XYME ;
XYM=XM+YM*XYM ;
END ;
ELSE DO ;
XYEX=YEX;
XYM=XEX-YEX; XYME=XYM;
IF XYM<-60 THEN XYM=0; ELSE XYM=10**XYME ;
XYM=XM*XYM+YM ;
END ;
IF ABS(XYM)>=10 THEN DO ;
XYM=XYM/10; XYEX=XYEX+1 ;
END ;
IF ABS(XYM)<1 THEN DO;
XYM=XYM*10; XYEX=XYEX-1;
END ;
END ADD ;
MULT: PROCEDURE(XM,XEX,YM,YEX,XYM,XYEX) ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
DCL (XM,YM,XYM) DEC FLOAT(16) ;
DCL (XEX,YEX,XYEX) DEC FIXED(15,0) ;
IF XM=0|YM=0 THEN DO ;
XYM=0; XYEX=0;
RETURN;
END ;
XYM=XM*YM ;
XYEX=XEX+YEX ;
IF ABS(XYM)>=10 THEN DO ;
XYM=XYM/10; XYEX=XYEX+1 ;
END ;
END MULT ;
DIV: PROCEDURE(XM,XEX,YM,YEX,XYM,XYEX) ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
DCL (XM,YM,XYM) DEC FLOAT(16) ;
DCL (XEX,YEX,XYEX) DEC FIXED(15,0) ;
IF XM=0 THEN DO;
XYM=0; XYEX=0;
RETURN;
END ;
XYM=XM/YM ;
XYEX=XEX-YEX ;
IF ABS(XYM)<1 THEN DO ;
XYM=XYM*10; XYEX=XYEX-1 ;
END ;
END DIV ;
|