/usr/share/gnudatalanguage/astrolib/trapzd.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 | 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
|