C+
C NAME:
C	Simpson
C PURPOSE:
C	Numerical integration of function of one variable using Simpson rule
C CALLING SEQUENCE:
	function Simpson(F,A,B,DEL,MAXN,N)
C INPUTS:
C	F	real	function; declared external in caller
C	A,B	real	integration range
C	DEL	real	accuracy
C	MAXN	integer	maximum number of iterations
C OUTPUTS:
C	Simpson	real	result of integration
C	N	integer	0: no convergence
C			1: if A=B (then also R=0)
C PROCEDURE:
C >	If DEL <=0 then DEL=1x10^-4 is used
C >	If MAXN< 2 then MAXN=2 is used
C MODIFICATION HISTORY:
C-
	    real	F
	    real	A
	    real	B
	    real	DEL
	    integer	MAXN
	    integer	N

	if (DEL .le. 0.) DEL = 1.E-4
	MAXI = max(2,MAXN)

	N = 1
	Simpson = 0.
	if (B .eq. A) return		! Zero integration range

	X  = (A+B)/2.
	BA = B-A
	M  = 1

	SM = F(X)*BA*2./3.
	S  = SM+(F(A)+F(B))*BA/6.
	SI = S*(1.+2*DEL)		! To get inside DO-loop

	I = 2
	do while (I .le. MAXI .and. abs(S-SI) .gt. abs(DEL*S))
	    SI = S
	    S  = (S-SM/2.)/2.
	    M  = M*2
	    DX = BA/M
	    SM = 0.
	    X0  = A-DX/2.
	    do IX=1,M
		SM = SM+F(X0+IX*DX)
	    end do
	    SM = SM*2.*BA/(3.*M)
	    S = S+SM
	    I = I+1
	end do

	if (I .gt. MAXI) M = 0		! No convergence
	N = 2*M

	Simpson = S

	return
	end