Sources Esope | Cast3M

zwgljd
C ZWGLJD    SOURCE    CHAT      05/01/13    04:25:30     5004
      SUBROUTINE ZWGLJD (Z,W,NP,ALPHA,BETA)
C--------------------------------------------------------------------
C
C     Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W)
C     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
C     The polynomial degree N=NP-1.
C     Double precision version.
C
C--------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Z(*),W(*),ALPHA,BETA
      N     = NP-1
      NM1   = N-1
      ONE   = 1.D0
      TWO   = 2.D0
C
      IF (NP.LE.1) THEN
         WRITE (6,*) 'Minimum number of Gauss-Lobatto points is 2'
         STOP
      ENDIF
      IF ((ALPHA.LE.-ONE).OR.(BETA.LE.-ONE)) THEN
         WRITE (6,*) 'Alpha and Beta must be greater than -1'
         STOP
      ENDIF
C
      IF (NM1.GT.0) THEN
         ALPG  = ALPHA+ONE
         BETG  = BETA+ONE
         CALL ZWGJD (Z(2),W(2),NM1,ALPG,BETG)
      ENDIF
      Z(1)  = -ONE
      Z(NP) =  ONE
      DO 100 I=2,NP-1
         W(I) = W(I)/(ONE-Z(I)**2)
 100  CONTINUE
      CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,Z(1))
      W(1)  = ENDW1 (N,ALPHA,BETA)/(TWO*PD)
      CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,Z(NP))
      W(NP) = ENDW2 (N,ALPHA,BETA)/(TWO*PD)
      RETURN
      END