/usr/share/gnudatalanguage/astrolib/qsimp.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 | 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
|