Esope source files

   1 : C ZWGLJD    SOURCE    CHAT      05/01/13    04:25:30     5004           
   2 :       SUBROUTINE ZWGLJD (Z,W,NP,ALPHA,BETA)
   3 : C--------------------------------------------------------------------
   4 : C
   5 : C     Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W)
   6 : C     associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
   7 : C     The polynomial degree N=NP-1.
   8 : C     Double precision version.
   9 : C
  10 : C--------------------------------------------------------------------
  11 :       IMPLICIT INTEGER(I-N)
  12 :       IMPLICIT REAL*8 (A-H,O-Z)
  13 :       REAL*8 Z(*),W(*),ALPHA,BETA
  14 :       N     = NP-1
  15 :       NM1   = N-1
  16 :       ONE   = 1.D0
  17 :       TWO   = 2.D0
  18 : C
  19 :       IF (NP.LE.1) THEN
  20 :          WRITE (6,*) 'Minimum number of Gauss-Lobatto points is 2'
  21 :          STOP
  22 :       ENDIF
  23 :       IF ((ALPHA.LE.-ONE).OR.(BETA.LE.-ONE)) THEN
  24 :          WRITE (6,*) 'Alpha and Beta must be greater than -1'
  25 :          STOP
  26 :       ENDIF
  27 : C
  28 :       IF (NM1.GT.0) THEN
  29 :          ALPG  = ALPHA+ONE
  30 :          BETG  = BETA+ONE
  30 :          CALL ZWGJD (Z(2),W(2),NM1,ALPG,BETG)
  31 :       ENDIF
  32 :       Z(1)  = -ONE
  33 :       Z(NP) =  ONE
  34 :       DO 100 I=2,NP-1
  35 :          W(I) = W(I)/(ONE-Z(I)**2)
  36 :  100  CONTINUE
  36 :       CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,Z(1))
  37 :       W(1)  = ENDW1 (N,ALPHA,BETA)/(TWO*PD)
  37 :       CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,Z(NP))
  38 :       W(NP) = ENDW2 (N,ALPHA,BETA)/(TWO*PD)
  39 :       RETURN
  40 :       END
  41 :  
  42 :