!  Sets up integration points for
!  Gaussian quadrature using Laguerre polynomials
SUBROUTINE gauss_laguerre(x,w,n,alf)
  IMPLICIT NONE
  INTEGER n,MAXIT
  REAL(KIND=8) alf,w(n),x(n), EPS
  PARAMETER (EPS=3.D-14,MAXIT=10)
  INTEGER i,its,j
  REAL(KIND=8) ai,gammln
  REAL(KIND=8) p1,p2,p3,pp,z,z1
  DO i=1,n
     IF(i == 1)THEN
        z=(1.+alf)*(3.+.92*alf)/(1.+2.4*n+1.8*alf)
     ELSE IF(i.EQ.2)THEN
        z=z+(15.+6.25*alf)/(1.+.9*alf+2.5*n)
     ELSE
        ai=i-2
        z=z+((1.+2.55*ai)/(1.9*ai)+1.26*ai*alf/(1.+3.5*ai))* &
          (z-x(i-2))/(1.+.3*alf)
     ENDIF
     DO its=1,MAXIT
        p1=1.d0
        p2=0.d0
        DO  j=1,n
           p3=p2
           p2=p1
           p1=((2*j-1+alf-z)*p2-(j-1+alf)*p3)/j
        ENDDO
        pp=(n*p1-(n+alf)*p2)/z
        z1=z
        z=z1-p1/pp
        IF(ABS(z-z1) <= EPS) GOTO 1
     ENDDO
     PAUSE 'too many iterations in gaulag'
1    x(i)=z
     w(i)=-EXP(gammln(alf+n)-gammln(float(n)))/(pp*n*p2)
  ENDDO

END SUBROUTINE gauss_laguerre

REAL(KIND =8) FUNCTION gammln(xx)
  REAL(KIND=8) xx
  INTEGER j
  REAL(KIND=8) ser,stp,tmp,x,y,cof(6)
  SAVE cof,stp
  DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, &
   24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, &
   -.5395239384953d-5,2.5066282746310005d0/
  x=xx
  y=x
  tmp=x+5.5d0
  tmp=(x+0.5d0)*LOG(tmp)-tmp
  ser=1.000000000190015d0
  DO j=1,6
     y=y+1.d0
     ser=ser+cof(j)/y
  ENDDO
  gammln=tmp+LOG(stp*ser/x)

END FUNCTION gammln