Esope source files

   1 : C ZWGJD     SOURCE    CHAT      05/01/13    04:25:27     5004           
   2 :       SUBROUTINE ZWGJD (Z,W,NP,ALPHA,BETA)
   3 : C--------------------------------------------------------------------
   4 : C
   5 : C     Generate NP GAUSS 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(1),W(1),ALPHA,BETA
  14 : C
  15 :       N     = NP-1
  16 :       DN    = DBLE(N)
  17 :       ONE   = 1.D0
  18 :       TWO   = 2.D0
  19 :       APB   = ALPHA+BETA
  20 : C
  21 :       IF (NP.LE.0) THEN
  22 :          WRITE (6,*) 'Minimum number of Gauss points is 1'
  23 :          STOP
  24 :       ENDIF
  25 :       IF ((ALPHA.LE.-ONE).OR.(BETA.LE.-ONE)) THEN
  26 :          WRITE (6,*) 'Alpha and Beta must be greater than -1'
  27 :          STOP
  28 :       ENDIF
  29 : C
  30 :       IF (NP.EQ.1) THEN
  31 :          Z(1) = (BETA-ALPHA)/(APB+TWO)
  32 :          W(1) = GAMMAF(ALPHA+ONE)*GAMMAF(BETA+ONE)/GAMMAF(APB+TWO)
  33 :      $          * TWO**(APB+ONE)
  34 :          RETURN
  35 :       ENDIF
  36 : C
  36 :       CALL JACG (Z,NP,ALPHA,BETA)
  37 : C
  38 :       NP1   = N+1
  39 :       NP2   = N+2
  40 :       DNP1  = DBLE(NP1)
  41 :       DNP2  = DBLE(NP2)
  42 :       FAC1  = DNP1+ALPHA+BETA+ONE
  43 :       FAC2  = FAC1+DNP1
  44 :       FAC3  = FAC2+ONE
  45 :       FNORM = PNORMJ(NP1,ALPHA,BETA)
  46 :       RCOEF = (FNORM*FAC2*FAC3)/(TWO*FAC1*DNP2)
  47 :       DO 100 I=1,NP
  47 :          CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,NP2,ALPHA,BETA,Z(I))
  48 :          W(I) = -RCOEF/(P*PDM1)
  49 :  100  CONTINUE
  50 :       RETURN
  51 :       END
  52 :  
  53 :