gauleg.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\tomawac\gauleg.f
00002 !
00047                         SUBROUTINE GAULEG
00048 !                       *****************
00049 !
00050      &( W_LEG , X_LEG , NPOIN )
00051 !
00052 !***********************************************************************
00053 ! TOMAWAC   V6P1                                   15/06/2011
00054 !***********************************************************************
00055 !
00056 !
00057 !
00058 !
00059 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00060 !| NPOIN          |-->| NUMBER OF INTEGRATION POINTS OVER OMEGA2
00061 !| X_LEG          |<--| ABSICSSAE FOR THE GAUSS-LEGENDRE QUADRATURE
00062 !| W_LEG          |<--| WEIGHTS FOR THE GAUSS-LEGENDRE QUADRATURE
00063 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00064 !
00065       USE DECLARATIONS_TOMAWAC, ONLY : PI
00066 !
00067       IMPLICIT NONE
00068 !
00069 !.....VARIABLES IN ARGUMENT
00070 !     """"""""""""""""""""
00071       INTEGER           NPOIN
00072       DOUBLE PRECISION  W_LEG(NPOIN) , X_LEG(NPOIN)
00073 !
00074 !.....LOCAL VARIABLES
00075 !     """""""""""""""""
00076       INTEGER           I     , M     , J
00077       DOUBLE PRECISION  EPS   , Z     , P1    , P2    , P3    ,
00078      &                  PP    , Z1
00079       PARAMETER        (EPS=3.D-14)
00080 !
00081 !
00082       M=(NPOIN+1)/2
00083       DO I=1,M
00084         Z=COS(PI*(DBLE(I)-0.25D0)/(DBLE(NPOIN)+0.5D0))
00085     1   CONTINUE
00086         P1=1.0D0
00087         P2=0.0D0
00088         DO J=1,NPOIN
00089           P3=P2
00090           P2=P1
00091           P1=((2.D0*DBLE(J)-1.D0)*Z*P2-(DBLE(J)-1.D0)*P3)/DBLE(J)
00092         ENDDO
00093         PP=DBLE(NPOIN)*(Z*P1-P2)/(Z*Z-1.D0)
00094         Z1=Z
00095         Z=Z-P1/PP
00096         IF (ABS(Z-Z1).GT.EPS) GOTO 1
00097         X_LEG(I)=-Z
00098         X_LEG(NPOIN+1-I)=Z
00099         W_LEG(I)=2.D0/((1.D0-Z**2)*PP**2)
00100         W_LEG(NPOIN+1-I)=W_LEG(I)
00101       ENDDO
00102 !
00103       RETURN
00104       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0