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

#  Copyright (c) 1995-2012 by Hamilton Laboratories.  All rights reserved.

#  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