pro qsimp, func, A, B, S, EPS=eps, MAX_ITER = max_iter, _EXTRA = _EXTRA ;+ ; NAME: ; QSIMP ; PURPOSE: ; Integrate using Simpson's rule to specified accuracy. ; EXPLANATION: ; Integrate a function to specified accuracy using the extended ; trapezoidal rule. Adapted from algorithm in Numerical Recipes, ; by Press et al. (1992, 2nd edition), Section 4.2. This procedure ; has been partly obsolete since IDL V3.5 with the introduction of the ; intrinsic function QSIMP(), but see notes below. ; ; CALLING SEQUENCE: ; QSIMP, func, A, B, S, [ EPS = , MAX_ITER =, _EXTRA = ] ; ; INPUTS: ; func - scalar string giving name of function of one variable to ; be integrated ; A,B - numeric scalars giving the lower and upper bound of the ; integration ; ; OUTPUTS: ; S - Scalar giving the approximation to the integral of the specified ; function between A and B. ; ; OPTIONAL KEYWORD PARAMETERS: ; EPS - scalar specifying the fractional accuracy before ending the ; iteration. Default = 1E-6 ; MAX_ITER - Integer specifying the total number iterations at which ; QSIMP will terminate even if the specified accuracy has not yet ; been met. The maximum number of function evaluations will be ; 2^(MAX_ITER). Default value is MAX_ITER = 20 ; ; Any other keywords are passed directly to the user-supplied function ; via the _EXTRA facility. ; NOTES: ; (1) The function QTRAP is robust way of doing integrals that are not ; very smooth. However, if the function has a continuous 3rd derivative ; then QSIMP will likely be more efficient at performing the integral. ; ; (2) QSIMP can be *much* faster than the intrinsic QSIMP() function (as ; of IDL V8.2.3). This is because the intrinsic QSIMP() function only ; requires that the user supplied function accept a *scalar* variable. ; Thus on the the 16th iteration, the intrinsic QSIMP() makes 32,767 ; calls to the user function, whereas this procedure makes one call ; with a 32,767 element vector. Also, unlike the intrinsic QSIMP(), this ; procedure allows keywords in the user-supplied function. ; ; (3) Since the intrinsic QSIMP() is a function, and this file contains a ; procedure, there should be no name conflict. ; EXAMPLE: ; Compute the integral of sin(x) from 0 to !PI/3. ; ; IDL> QSIMP, 'sin', 0, !PI/3, S & print, S ; ; The value obtained should be cos(!PI/3) = 0.5 ; ; PROCEDURES CALLED: ; SETDEFAULTVALUE, TRAPZD, ZPARCHECK ; ; REVISION HISTORY: ; W. Landsman ST Systems Co. August, 1991 ; Continue after max iter warning message W. Landsman March, 1996 ; Pass keyword to function via _EXTRA facility W. Landsman July 1999 ; Use SETDEFAULTVALUE W. Landsman Aug 2013 ;- On_error,2 ;Return to caller compile_opt idl2 if N_params() LT 4 then begin print,'Syntax - QSIMP, func, A, B, S, [ MAX_ITER = , EPS = ]' print,' func - scalar string giving function name' print,' A,B - endpoints of integration, S - output sum' return endif zparcheck, 'QSIMP', func, 1, 7, 0, 'Function name' ;Valid inputs? zparcheck, 'QSIMP', A, 2, [1,2,3,4,5], 0, 'Lower limit of Integral' zparcheck, 'QSIMP', B, 3, [1,2,3,4,5], 0, 'Upper limit of Integral' setdefaultvalue,eps,1.e-6 ;Typo fixed Oct 2013 setdefaultvalue,max_iter,20 ost = (oS = -1.e30) for i = 0,max_iter - 1 do begin trapzd, func, A,B, st, it, _EXTRA = _EXTRA S = (4.*st - ost)/3. if ( abs(S-oS) LT eps*abs(oS) ) then return os = s ost = st endfor message,/CON, \$ 'WARNING - Sum did not converge after '+ strtrim(max_iter,2) + ' steps' return end