PRO	TRAPZD,FUNC,A,B,S,N
;
; Compute the Nth stage of refinement of an extended trapezoidal rule
;
;  Input:
;     FUNC - Name of multi-valued function to be integrated
;     A    - Lower limit of integral
;     B    - Upper limit of integral
;     N    - Index for stage of refinement
;
;  Input/Output:
;     S    - Result of integral (array)
;            S should not be modified between calls.
;
	IF N EQ 1 THEN $
	   S = 0.5D0*(B-A)*(CALL_FUNCTION(FUNC,A) + CALL_FUNCTION(FUNC,B)) $
	ELSE BEGIN
	   IT = 2L^(N-2)
	   TNM = DOUBLE(IT)
	   DEL = (B-A)/TNM
	   X = A + 0.5D0*DEL
	   SUM = CALL_FUNCTION(FUNC,X)
	   IF IT GT 1 THEN FOR J=2L,IT DO BEGIN
	      X = X + DEL
	      SUM = SUM + CALL_FUNCTION(FUNC,X)
	   ENDFOR
	   S = 0.5D0*(S + (B-A)*SUM/TNM)
	ENDELSE
	RETURN
	END
FUNCTION QSIMPS,FUNC,A,B,FERR,EPS
;
; Integration of multi-valued function FUNC from A to B using
; Simpson's rule.  For integration of single-valued functions,
; use IDL function QSIMP.
;
;  Input:
;     FUNC - Name of function to be integrated
;     A    - Lower limit of integral
;     B    - Upper limit of integral
;     FERR - Name of error function
;     EPS  - Array of maximum relative errors between interations
;
;  Output:
;     Array returned as function value
;
;  Make sure limits of integration are double precision
	DA = DOUBLE(A)
	DB = DOUBLE(B)
;
;  Initialization
	JMAX = 24
	OST = [-1.D30]
	OS  = [-1.D30]
;
;  Iterate until all integrals converge
	FOR J=1,JMAX DO BEGIN
	   TRAPZD,FUNC,DA,DB,ST,J
	   S = (4.D0*ST - OST)/3.D0
	   STATUS = CALL_FUNCTION(FERR,S,OS,EPS,J)
	   IF STATUS EQ 1 THEN RETURN,S
	   OS = S
	   OST = ST
	ENDFOR
	STOP,'QSIMPS: too many steps'
	END