Télécharger jacg.eso

Retour à la liste

Numérotation des lignes :

  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
  29. CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,NP,ALPHA,BETA,X)
  30. RECSUM = 0.D0
  31. JM = J-1
  32. DO 29 I=1,JM
  33. RECSUM = RECSUM+1.D0/(X-XJAC(NP-I+1))
  34. 29 CONTINUE
  35. DELX = -P/(PD-RECSUM*P)
  36. X = X+DELX
  37. IF (ABS(DELX) .LT. EPS) GOTO 31
  38. 30 CONTINUE
  39. 31 CONTINUE
  40. XJAC(NP-J+1) = X
  41. XLAST = X
  42. 40 CONTINUE
  43. DO 200 I=1,NP
  44. XMIN = 2.D0
  45. DO 100 J=I,NP
  46. IF (XJAC(J).LT.XMIN) THEN
  47. XMIN = XJAC(J)
  48. JMIN = J
  49. ENDIF
  50. 100 CONTINUE
  51. IF (JMIN.NE.I) THEN
  52. SWAP = XJAC(I)
  53. XJAC(I) = XJAC(JMIN)
  54. XJAC(JMIN) = SWAP
  55. ENDIF
  56. 200 CONTINUE
  57. RETURN
  58. END
  59.  
  60.  
  61.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales