Subroutine Int_Simp(f,a,b,s,n)
c=========================================================
c int_trap.f: integration by Simpson rule of f(x) on [a,b]
c method:     Simpson rule
c written by: Alex Godunov (October 2006)
c---------------------------------------------------------
c f     - Function to integrate (supplied by a user)
c a	- Lower limit of integration
c b	- Upper limit of integration
c s     - Result of integration (out)
c n	- number of intervals (should be even)
c=========================================================
        Real*8 a, b, f, s ,dx, x
	Integer*4 n, i
c
c if n is odd we add +1 to make it even
c
        if((n/2)*2.ne.n) n=n+1
        s = 0.0
	dx = (b-a)/float(n)
	Do i=2, n-1, 2
            x = a+float(i)*dx
            s = s + 2.0*f(x) + 4.0*f(x+dx)
	End do
        s = (s + f(a) + f(b) + 4.0*f(a+dx))*dx/3.0
	Return
	End