idlastro / Math and Statistics: QSIMP

[Source code]

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