/*************************************************************************************
*   This subroutine calculates the radial wave functions for a hydrogen-like bound atom.
*   It is written in PL/1,  but should be easily adaptable to other languages.
*
*    It is based on the well known analytical formula
*      Ψn,l(ρ)= √(n-l-1)!/√(n+l)!.√Z /n.(2Z.ρ/n)l+1.Ln-l-12l+1(2Z.ρ/n).exp(-Z.ρ/n) ;
*
*
*  Input parameters are:
*    Z:  nuclear charge
*    N0: principal quantum number
*    L0: angular momentum qunatum number
*    R:  radial distance (in Bohr radii)
*
*  Output parameter is:
*    PNL:  wave function value
*
*
*  Required subroutines: 
*     LAGUERRE2 (returns the corresponding Laguerre polynomials)
*     NUM, MULT (to avoid numerical overflow or underflow for high quantum numbers).
*
*
*   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
*********************************************************************************/

 COULOMBB2: PROCEDURE(Z,N0,L0,R,PNL) ;
 DEFAULT RANGE(A:Z) DEC FLOAT VALUE (DEC FLOAT(16)) ;
 DCL (N0,L0,ERO2EX,FAKEX,FAKEXL,SQNEX,LEX) DEC FIXED(15,0) ;
 DCL (LAGUERRE2,NUM,MULT) ENTRY ; 

 ON UNDERFLOW ; 

 LOG10E=4.342944819032518E-01 ; 
 ZN=2*Z/N0; RO=ZN*R; RLE=RO/2*LOG10E ;
 ERO2EX=-CEIL(RLE) ;
 ERO2=10**(-(RLE+ERO2EX)) ;
 FAKEX=0; SQNEX=0; L=0; LEX=0;

 CALL LAGUERRE2(RO,N0,L0,L,LEX) ; 

 N1=N0; FAK=RO/SQRT(N1) ;
 CALL NUM(FAK,FAKEX) ;
 IF L0=0 THEN GO TO CONT ;

 DO N=1 TO L0 ; 
 N1=N ; 
 SQN=RO/SQRT((N0+N1)*(N0-N1)) ; 
 CALL NUM(SQN,SQNEX) ;
 FAKL=FAK; FAKEXL=FAKEX;
 CALL MULT(FAKL,FAKEXL,SQN,SQNEX,FAK,FAKEX) ;
 END ;

 CONT:LX=FAKEX+ERO2EX+LEX ;
 PNL=FAK*ERO2*L*10**LX ;
 PNL=PNL*SQRT(Z)/N0 ;


 END COULOMBB2 ;	




/********************************************************************************
*  The subroutine for the calculation of the Laguerre polynomials is based on the
*  recursive relationships (Abramowitz and Stegun, Handbook of Mathematical Functions):
*
*  L0q(x)=1
*  L1q(x)=q+1-x
*  Lpq(x)=1/p.[(2p+q-1-x).Lp-1q(x) -(p-1+q).Lp-2q(x)]
*
*  Input parameters are:
*    X: argument of the Laguerre polynomial
*    N0: principal quantum number
*    L0: angular momentum quantum number
*
*  Output parameters are:
*    L: mantissa of Laguerre-polynomial
*    LEX: exponent of Laguerre-polynomial
*
*  Required subroutines: 
*     NUM, ADD, MULT, DIV 
*        (to avoid numerical overflow or underflow for large indices).
*    
*****************************************************************************/

 LAGUERRE2: PROCEDURE(X,N0,L0,L,LEX) ;
 DEFAULT RANGE(A:Z) DEC FLOAT VALUE (DEC FLOAT(16)) ;
 DCL (N0,L0,I,J,IJEX,XJEX,JIEX,LL2EX,LL1EX,L2EX,L1EX,L12EX,LEX) 
 DEC FIXED(15,0) ;
 DCL (NUM,ADD,MULT,DIV) ENTRY ; 

 ON UNDERFLOW ; 
 L=0; LEX=0; LL1EX=0; IJEX=0; XJEX=0; JIEX=0; L1EX=0; L2EX=0; L12EX=0;

 IF L0=N0-1 THEN DO ;
 L=1; LEX=0;
 GO TO ENDE ;
 END ;

 I=L0+1 ;
 LL2=1; LL2EX=0;
 LL1=2*I-X ; CALL NUM(LL1,LL1EX) ;
 IF L0=N0-2 THEN DO ;
 L=LL1; LEX=LL1EX;
 GO TO ENDE ;
 END ;
 DO J=I+2 TO N0 ;
 IJ=I+J-2 ;
 CALL NUM(IJ,IJEX) ;
 XJ=2*J-2-X ;
 CALL NUM(XJ,XJEX) ;
 JI=J-I ;
 CALL NUM(JI,JIEX) ;
 CALL MULT(XJ,XJEX,LL1,LL1EX,L1,L1EX) ; 
 CALL MULT(IJ,IJEX,LL2,LL2EX,L2,L2EX) ; 
 L2=-L2 ;
 CALL ADD(L1,L1EX,L2,L2EX,L12,L12EX) ;
 CALL DIV(L12,L12EX,JI,JIEX,L,LEX) ;
 LL2=LL1; LL2EX=LL1EX;
 LL1=L; LL1EX=LEX;
 END ;

 ENDE: END LAGUERRE2 ;




/***********************************************************************************
*  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 ;