#	Trapeziodal rule numeric integration of function f() from a to b,
#	using n intervals.  Release 2.2.49.

#	Sample usage:
#	Integrate the sin() function from 0 to 1, using 500 intervals

#		trapz	sin 0 1 500

#	Compare this result with the "exact" integration:

#		calc cos(0) - cos(1)

#	The area in each interval is (f(x) + f(x + dx))*dx/2.  Summing
#	the areas of all the intervals gives the total area, but a
#	simple optimization is to move the multiplication by dx out of
#	the loop and to observe that the divide by 2 is unnecessary
#	except for the values of f() at the endpoints a and b since,
#	everywhere else, the end of one interval is the beginning of
#	the next.

#	Note the use of eval lets me paste in the name of the function
#	for the call, then cause the C shell to compile and run the
#	resulting string.  The circumflex characters at the end of each
#	line are there just to let me break the string across multiple
#	lines for readability.


proc trapz( f, a, b, n)
	local i, area, dx
	@ dx = (b - a)/n
	eval	"@ area = ($f (a) + $f (b))/2	^
			for i = 1 to n - 1 do			^
				@ area += $f (a + i*dx)		^
			end"
	@ area *= dx
	return area
end

trapz $argv
