/********************************************************************************
*  This PL/1 subroutine calculates of the Laguerre polynomials, 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)]
*
*  Since the program was written for calculating the bound atomic wave functions,
*  it returns in fact not just one value but all N values  L02N-1, L12N-2,...LN-11
*
*  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
*
*
*   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
*****************************************************************************/

 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 ;