function nrSimpson, F,A,B,eps=DEL,maxn=MAXN,n=N,arg=arg

@compile_opt.pro			; On error, return to caller

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

if n_elements(DEL ) eq 0 then DEL = 1.E-4
if n_elements(MAXN) eq 0 then MAXN = 13

MAXI = MAXN > 2

N = 1
Fnc = 0.
if B eq A then return, Fnc		; Zero integration range

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

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

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

if I gt MAXI then M = 0			; No convergence
N = 2*M

Fnc = S

return, Fnc  &  end