pro trapzd, func, a, b, s, step, _EXTRA = _EXTRA
;+
; NAME:
;       TRAPZD
; PURPOSE:
;       Compute the nth stage of refinement of an extended trapezoidal rule.
; EXPLANATION:
;       This procedure is called by QSIMP and QTRAP.   Algorithm from Numerical
;       Recipes, Section 4.2.   TRAPZD is meant to be called iteratively from
;       a higher level procedure.
;
; CALLING SEQUENCE:
;       TRAPZD, func, A, B, S, step, [ _EXTRA = ]
;
; INPUTS:
;       func - scalar string giving name of function to be integrated.   This
;               must be a function of one variable.
;       A,B -  scalars giving the limits of the integration
;
; INPUT-OUTPUT:
;       S -    scalar giving the total sum from the previous iterations on 
;               input and the refined sum after the current iteration on output.
;
;       step - LONG scalar giving the number of points at which to compute the
;               function for the current iteration.   If step is not defined on
;               input, then S is intialized using the average of the endpoints
;               of limits of integration.
;
; OPTIONAL INPUT KEYWORDS:
;       Any supplied keywords will be passed to the user function via the 
;       _EXTRA facility. 
;
; NOTES:
;       (1) TRAPZD will check for math errors (except for underflow) when 
;       computing the function at the endpoints, but not on subsequent 
;       iterations.
;
;       (2) TRAPZD always uses double precision to sum the function values
;       but the call to the user-supplied function is double precision only if 
;       one of the limits A or B is double precision.
; REVISION HISTORY:
;       Written         W. Landsman                 August, 1991
;       Always use double precision for TOTAL       March, 1996
;       Pass keyword to function via _EXTRA facility  W. Landsman July 1999
;       Don't check for floating underflow  W.Landsman  April 2008
;-
 On_error,2
 compile_opt idl2

 kpresent = keyword_set(_EXTRA)
 if N_elements(step) EQ 0 then begin          ;Initialize?

;If a math error occurs, it is likely to occur at the endpoints
     junk = check_math()                    ;
     if kpresent then s1 = CALL_FUNCTION(func,A, _EXTRA= _EXTRA) $
                 else s1 = CALL_FUNCTION(func,A)
     if check_math(mask=211) NE 0 then $
        message,'ERROR - Illegal lower bound of '+strtrim(A,2)+ $
                ' to function ' + strupcase(func)
     if kpresent then s2 = CALL_FUNCTION(func,B, _EXTRA = _EXTRA) $
                 else s2 = CALL_FUNCTION(func,B)
     if check_math(mask=211) NE 0 then $
        message,'ERROR - Illegal upper bound of '+strtrim(B,2) + $
                ' to function ' + strupcase(func)
     junk= check_math()		
     s = 0.5d * ( double(B)-A ) * ( s1+s2 )    ;First approx is average of endpoints
     step = 1l

 endif else begin

     tnm = float( step )               
     del = ( B - A ) / tnm                    ;Spacing of the points to add
     x = A + 0.5*del + findgen( step ) * del  ;Grid of points @ compute function
     if kpresent then sum = CALL_FUNCTION( func, x, _EXTRA = _EXTRA) $
                 else sum = CALL_FUNCTION( func, x)
     S = 0.5d * ( S + (double(B)-A) * total( sum, /DOUBLE )/tnm )     
     step = 2*step

 endelse

 return
 end