BeginPackage["ISP`Integ`"]

AnimateRectangular::usage = "AnimateRectangular[f, low, high, n] shows
how the rectangular method approximates the integral of f between low
and high as the number of rectangles varies from 2 to 2^n"

AnimateTrapezoidal::usage = "AnimateTrapezoidal[f, low, high, n] shows
how the trapezoidal method approximates the integral of f between low
and high as the number of trapezoids varies from 2 to 2^n"

Begin["Private`"]

onerect[h_, n_, f_, lo_] :=
  {{lo+(n-1)*h, 0},
   {lo+(n-1)*h, f[lo+n*h]},
   {lo+n*h, f[lo+n*h]}};

rect[func_, step_, steps_, lo_, hi_, thePlot_] :=
  Module[{n, x, aPlot},
         (aPlot = ListPlot[Join[Apply[Join, Table[onerect[step, n, func, lo],
                                                  {n, 1, steps}]],
                                {{lo+steps*step, 0}}],
                           DisplayFunction->Identity,
                           PlotJoined->True];
          Show[{aPlot, thePlot}, DisplayFunction->$DisplayFunction])]


AnimateRectangular[integrand_, lo_, hi_, nframes_] :=
  Module[{frame, thePlot},
         (thePlot = Plot[integrand[x], {x, lo, hi},
                         DisplayFunction->Identity];
         Table[rect[integrand, (hi-lo)/2^frame, 2^frame, lo, hi,
                    thePlot],
               {frame, 1, nframes}])]

onetrap[h_, n_, f_, lo_] :=
   {{lo+(n-1)*h, 0},
    {lo+(n-1)*h, f[lo+(n-1)*h]},
    {lo+n*h, f[lo+n*h]}}

trap[func_, step_, steps_, lo_, hi_, thePlot_] :=
  Module[{n, x, aPlot},
         (aPlot = ListPlot[Join[Apply[Join, Table[onetrap[step, n, func, lo],
                                                  {n, 1, steps}]],
                                {{lo+steps*step, 0}}],
                           DisplayFunction->Identity,
                           PlotJoined->True];
          Show[{aPlot, thePlot}, DisplayFunction->$DisplayFunction])]

AnimateTrapezoidal[integrand_, lo_, hi_, nframes_] :=
  Module[{frame, thePlot},
         (thePlot = Plot[integrand[x], {x, lo, hi},
                         DisplayFunction->Identity];
         Table[trap[integrand, (hi-lo)/2^frame, 2^frame, lo, hi,
                    thePlot],
               {frame, 1, nframes}])]
    

End[]

EndPackage[]
