function integrate, x, f
;;; Deal with the edge effects by forcing the edge values to have
;;; zero weight in the sum. In other words, dx_0 and dx_(N-1) are
;;; zero. 
xpad = [x[0], x, max(x)]
;;; Make sure f has the same length as the new x value array
fpad = [0., f, 0.]

dx = shift(xpad, -1) - xpad
arg = 0.5 * (shift(fpad, -1) + fpad)
arg *= dx
;;; Trim fat
a = arg[1:n_elements(xpad)-3]
return, total(a)
end