/*************************************************************************************
*   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
*    IRA: number of radial points for which the function is to be calculated
*            (dimension of vector RPOINT)
*    RPOINT: Vector containing the radial points (in Bohr radii)
*
*  Output parameter is:
*    PNL:  2-dim array containing the wave function values for all radius points
*           and angular momenta
*
*  Required subroutine: LAGUERRE (returns the corresponding Laguerre polynomials)
*
*   Notes:  for very high principal quantum numbers (i.e. N0 into the hundreds),
*   some of the numerical terms here can become too large/small
*   for a numerical representation in the computer ( I have written a specially adapted
*   alternative routine for high quantum numbers; see the program COULOMBB2).
*
*   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, March 2007
*********************************************************************************/

 COULOMBB: PROCEDURE (Z,N0,IRA,RPOINT,PNL)  ;
 DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
 DCL (RPOINT(*),LL(*),PNL(*,*)) CTL DEC FLOAT(16) ;
 DCL (N0,J,IRA,N) DEC FIXED(7,0) ;
 DCL LAGUERRE ENTRY ;

 ON UNDERFLOW ;

 ALLOCATE LL(N0) ;   /* Allocation of Laguerre polynomial vector */


 ZN=2*Z/N0 ;



 DO J=1 TO IRA ;   /* Loop for all radius points */

     R=RPOINT(J) ;
     RO=ZN*R ;
     ERO2=EXP(-RO/2) ;
     CALL LAGUERRE(RO,N0,LL) ;
     N1=N0 ;
     FAK=RO/SQRT(N1) ;
     PNL(1,J)=LL(1)*FAK ;

     DO N=2 TO N0 ;  /* recursive calculation of wave function for all angular momenta */
     N1=N-1 ;
     FAK=FAK*RO/SQRT((N0+N1)*(N0-N1)) ;
     PNL(N,J)=LL(N)*FAK ;
     END ;

     PNL(*,J)=PNL(*,J)*ERO2 ;

 END ;


 PNL=PNL*SQRT(Z)/N0 ;   /* normalization of array elements */

 FREE LL ;   /* release of allocate storage space */

 END COULOMBB ;





/********************************************************************************
*  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
*    N: principal quantum number (= number of angular momentum values)
*
*  Output parameter is:
*    NLL: vector containing the values of the Laguerre polynomials
*         for the N angular momenta
*****************************************************************************/

 LAGUERRE : PROCEDURE(X,N,NLL) ;
 DEFAULT RANGE(A:Z) DEC FLOAT VALUE (DEC FLOAT(16)) ;
 DCL (L(*,*),NLL(*)) CTL DEC FLOAT(16) ;
 DCL (N,I,J,IJ,J1,J2) DEC FIXED(7,0) ;

 ON UNDERFLOW ;

 ALLOCATE L(N,N) ;

 L=0 ;

 DO I=1 TO N ;
 L(I,I)=1 ;
 END ;

 DO I=1 TO N-1 ;
 L(I,I+1)=2*I-X ;
 END ;

 DO I=1 TO N-2 ;
 DO J=I+2 TO N;
 J1=J-1; J2=J-2; IJ=I+J-2 ;
 L(I,J)=((2*J-2-X)*L(I,J1)-IJ*L(I,J2))/(J-I) ;
 END ;
 END ;

 DO I=1 TO N  ;
 NLL(I)=L(I,N) ;
 END ;

 FREE L ;

 END LAGUERRE ;