gauleg.f

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE gauleg
!                       *****************
!
     &( w_leg , x_leg , npoin )
!
!***********************************************************************
! TOMAWAC   V6P1                                   15/06/2011
!***********************************************************************
!
!brief  COMPUTES WEIGHTS AND ABSICSSA FOR THE GAUSS-LEGENDRE QUADRATURE.
!+        THE OUTPUT ABSCISSA ARE INCLUDED BETWEEN -1 AND 1.
!+        SUBROUTINE CALLED BY PRENL3
!
!history  E. GAGNAIRE-RENOU
!+        04/2011
!+        V6P1
!+   CREATED
!
!history  G.MATTAROLO (EDF - LNHE)
!+        15/06/2011
!+        V6P1
!+   Translation of French names of the variables in argument
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| NPOIN          |-->| NUMBER OF INTEGRATION POINTS OVER OMEGA2
!| X_LEG          |<--| ABSICSSAE FOR THE GAUSS-LEGENDRE QUADRATURE
!| W_LEG          |<--| WEIGHTS FOR THE GAUSS-LEGENDRE QUADRATURE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_tomawac, ONLY : pi
!
      USE interface_tomawac, ex_gauleg => gauleg
      IMPLICIT NONE
!
!.....VARIABLES IN ARGUMENT
!     """"""""""""""""""""
      INTEGER ,INTENT(IN)             ::          NPOIN
      DOUBLE PRECISION ,INTENT(INOUT) :: W_LEG(npoin) , X_LEG(npoin)
!
!.....LOCAL VARIABLES
!     """""""""""""""""
      INTEGER           I     , M     , J
      DOUBLE PRECISION  Z     , P1    , P2    , P3    , PP    , Z1
      DOUBLE PRECISION, PARAMETER  :: EPS = 3.d-14
!
      m=(npoin+1)/2
      DO i=1,m
        z=cos(pi*(dble(i)-0.25d0)/(dble(npoin)+0.5d0))
    1   CONTINUE
        p1=1.0d0
        p2=0.0d0
        DO j=1,npoin
          p3=p2
          p2=p1
          p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)
        ENDDO
        pp=dble(npoin)*(z*p1-p2)/(z*z-1.d0)
        z1=z
        z=z-p1/pp
        IF (abs(z-z1).GT.eps) GOTO 1
        x_leg(i)=-z
        x_leg(npoin+1-i)=z
        w_leg(i)=2.d0/((1.d0-z**2)*pp**2)
        w_leg(npoin+1-i)=w_leg(i)
      ENDDO
!
      RETURN
      END