1 : C JACG SOURCE CHAT 05/01/13 00:48:05 5004 2 : SUBROUTINE JACG (XJAC,NP,ALPHA,BETA) 3 : C-------------------------------------------------------------------- 4 : C 5 : C Compute NP Gauss points XJAC, which are the zeros of the 6 : C Jacobi polynomial J(NP) with parameters ALPHA and BETA. 7 : C ALPHA and BETA determines the specific type of Gauss points. 8 : C Examples: 9 : C ALPHA = BETA = 0.0 -> Legendre points 10 : C ALPHA = BETA = -0.5 -> Chebyshev points 11 : C 12 : C-------------------------------------------------------------------- 13 : IMPLICIT INTEGER(I-N) 14 : IMPLICIT REAL*8 (A-H,O-Z) 15 : REAL*8 XJAC(1) 16 : DATA KSTOP/10/ 17 : DATA EPS/1.0D-12/ 18 : N = NP-1 19 : DTH = 4.D0*ATAN(1.D0)/(2.D0*N+2.D0) 20 : DO 40 J=1,NP 21 : IF (J.EQ.1) THEN 22 : X = COS((2.D0*(J-1.D0)+1.D0)*DTH) 23 : ELSE 24 : X1 = COS((2.D0*(J-1.D0)+1.D0)*DTH) 25 : X2 = XLAST 26 : X = (X1+X2)/2.D0 27 : ENDIF 28 : DO 30 K=1,KSTOP 28 : CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,NP,ALPHA,BETA,X) 29 : RECSUM = 0.D0 30 : JM = J-1 31 : DO 29 I=1,JM 32 : RECSUM = RECSUM+1.D0/(X-XJAC(NP-I+1)) 33 : 29 CONTINUE 34 : DELX = -P/(PD-RECSUM*P) 35 : X = X+DELX 36 : IF (ABS(DELX) .LT. EPS) GOTO 31 37 : 30 CONTINUE 38 : 31 CONTINUE 39 : XJAC(NP-J+1) = X 40 : XLAST = X 41 : 40 CONTINUE 42 : DO 200 I=1,NP 43 : XMIN = 2.D0 44 : DO 100 J=I,NP 45 : IF (XJAC(J).LT.XMIN) THEN 46 : XMIN = XJAC(J) 47 : JMIN = J 48 : ENDIF 49 : 100 CONTINUE 50 : IF (JMIN.NE.I) THEN 51 : SWAP = XJAC(I) 52 : XJAC(I) = XJAC(JMIN) 53 : XJAC(JMIN) = SWAP 54 : ENDIF 55 : 200 CONTINUE 56 : RETURN 57 : END 58 : 59 :
© Cast3M 2003 - All rights reserved.
Disclaimer