pro newsimp, func, A, B, S, EPS=eps, MAX_ITER = max_iter, _EXTRA = _EXTRA ;+ ; NAME: ; NEWSIMP ; ; PURPOSE: ; Integrate using Simpson's rule to specified accuracy. ; ; EXPLANATION: ; Integrate a function of one variable to specified accuracy using the extended ; trapezoidal rule. Adapted from algorithm in Numerical Recipes, ; by Press et al. (1992, 2nd edition), Section 4.2. A earlier version of ; this procedure was made the intrinsic function QSIMP() in IDL V3.5. ; See notes below for distinctions. ; ; CALLING SEQUENCE: ; NEWSIMP, func, A, B, S, [ EPS = , MAX_ITER =, _EXTRA = ] ; ; INPUTS: ; func - Scalar STRING giving the name of the IDL function routine to ; be integrated (one variable only) ; 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 ; NEWSIMP 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 a more robust way of doing integrals that are not ; very smooth. However, if the function has a continuous 3rd derivative ; then NEWSIMP will likely be more efficient at performing the integral. ; ; (2) NEWSIMP can be *much* faster than the intrinsic QSIMP() function (as ; of IDL V5.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> NEWSIMP, 'sin', 0, !PI/3, S & print, S ; ; The value obtained should be cos(!PI/3) = 0.5 ; ; PROCEDURES CALLED: ; TRAPZD, ZPARCHECK ; ; REVISION HISTORY: ; W. Landsman ST Systems Co. August, 1991 ; Continue after max iter warning message W. Landsman March, 1996 ; Converted to IDL V5.0 W. Landsman September 1997 ; Pass keyword to function via _EXTRA facility W. Landsman July 1999 ; Change name to NEWSIMP for IDL Exercises, R. O'Connell, Sept 2003 ;- On_error,2 ;Return to caller if N_params() LT 4 then begin print,'Syntax - newsimp, 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, 'newsimp', func, 1, 7, 0, 'Function name' ;Valid inputs? zparcheck, 'newsimp', A, 2, [1,2,3,4,5], 0, 'Lower limit of Integral' zparcheck, 'newsimp', B, 3, [1,2,3,4,5], 0, 'Upper limit of Integral' If not keyword_set(EPS) then eps = 1.e-6 ;Set defaults if not keyword_set(MAX_ITER) then 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