C	  Integration par la methode de Simpson
C	  Dominique Lefebvre  janvier 2007
C  
C     a = borne inferieure d'integration
C     b = borne superieure d'integration
C     n = nombre de pas
C     aire = surface retournee

		SUBROUTINE IntSimpson(fn,a,b,n,aire)

		REAL a,b,fn,aire
		INTEGER n
		EXTERNAL fn

		REAL h,SommePaire, SommeImpaire
        INTEGER i

C     Boucle d'integration

C	  Initialisation des variables
    	aire = 0
        h = (b-a)/(n*2)
        SommePaire = 0
        SommeImpaire = 0

C     Calcul de la somme des indices impaires
 		DO i=1, n-1
          SommeImpaire = SommeImpaire + fn(a+h*2*i)
        ENDDO

C     Calcul de la somme des indices paires
 		DO i=1, n
          SommePaire = SommePaire + fn(a+h*(2*i-1))
        ENDDO


C     Calcul final de l'aire
        aire = h*(fn(a) + fn(b)+ 2*SommePaire + 4*SommeImpaire)/3
   
		END