Sources Esope | Cast3M

jacg
C JACG      SOURCE    CHAT      05/01/13    00:48:05     5004
      SUBROUTINE JACG (XJAC,NP,ALPHA,BETA)
C--------------------------------------------------------------------
C
C     Compute NP Gauss points XJAC, which are the zeros of the
C     Jacobi polynomial J(NP) with parameters ALPHA and BETA.
C     ALPHA and BETA determines the specific type of Gauss points.
C     Examples:
C     ALPHA = BETA =  0.0  ->  Legendre points
C     ALPHA = BETA = -0.5  ->  Chebyshev points
C
C--------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 XJAC(1)
      DATA KSTOP/10/
      DATA EPS/1.0D-12/
      N   = NP-1
      DTH = 4.D0*ATAN(1.D0)/(2.D0*N+2.D0)
      DO 40 J=1,NP
         IF (J.EQ.1) THEN
            X = COS((2.D0*(J-1.D0)+1.D0)*DTH)
         ELSE
            X1 = COS((2.D0*(J-1.D0)+1.D0)*DTH)
            X2 = XLAST
            X  = (X1+X2)/2.D0
         ENDIF
         DO 30 K=1,KSTOP
            CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,NP,ALPHA,BETA,X)
            RECSUM = 0.D0
            JM = J-1
            DO 29 I=1,JM
               RECSUM = RECSUM+1.D0/(X-XJAC(NP-I+1))
 29         CONTINUE
            DELX = -P/(PD-RECSUM*P)
            X    = X+DELX
            IF (ABS(DELX) .LT. EPS) GOTO 31
 30      CONTINUE
 31      CONTINUE
         XJAC(NP-J+1) = X
         XLAST        = X
 40   CONTINUE
      DO 200 I=1,NP
         XMIN = 2.D0
         DO 100 J=I,NP
            IF (XJAC(J).LT.XMIN) THEN
               XMIN = XJAC(J)
               JMIN = J
            ENDIF
 100     CONTINUE
         IF (JMIN.NE.I) THEN
            SWAP = XJAC(I)
            XJAC(I) = XJAC(JMIN)
            XJAC(JMIN) = SWAP
         ENDIF
 200  CONTINUE
      RETURN
      END