BeginPackage["Graphics`OscillatorPlot`"]

Needs["Utilities`FilterOptions`"]

Remove[OscillatorPlot];

OscillatorPlot::usage = "OscillatorPlot[f[x],{x,a,b},n1,n2] plots
the energy representation of a given function in the basis of
harmonic oscillator eigenfunctions.
Numerical region for the determination of the expansion
coefficients is the interval (a,b).
The graph shows the expansion coefficients c_n
for n between n1 and n2. Default values are:
a=-Infinity, b=Infinity, n1=0, n2=10."

Begin["`Private`"]

OscillatorFunction[n_,x_] :=
   HermiteH[n,x]*Exp[-x^2/2]/Sqrt[2^n*n!]/Pi^(1/4)

Options[OscillatorPlot] = Join[Options[Graphics],
Options[NIntegrate],{PlotStyle->Thickness[0.01]}];

SetAttributes[OscillatorPlot,HoldAll];

OscillatorPlot[f_, {x_Symbol, x1_:-Infinity, x2_:Infinity}, opts___Rule] :=
   OscillatorPlot[f, {x, x1, x2}, {0, 10}, opts]

OscillatorPlot[f_, {x_Symbol, x1_:-Infinity, x2_:Infinity},
    {n1_, n2_}, opts___Rule] := 
  Module[{func, coeff, coloredLine, lineTable, up, upper, lower},
    func = Function[{x}, f];
    style = PlotStyle/.{opts}/.Options[OscillatorPlot]; 
    Do[coeff[n] = NIntegrate[OscillatorFunction[n, x]*func[x],
        {x, x1, x2}, Evaluate[FilterOptions[NIntegrate,opts]]],
        {n, n1, n2}]; 
    coloredLine[n_] := {style, 
      Hue[Arg[coeff[n]]/(2*Pi)], 
      Line[{{n + 1/2, Abs[coeff[n]]}, 
        {n + 1/2, 0.}}]};
    lineTable = 
     Graphics[Table[coloredLine[n], {n, n1, n2}]]; 
    up = Max[Table[Abs[coeff[n]], {n, n1, n2}]]; 
    upper = up + up/10; lower = -(up/10); 
    Show[lineTable, Evaluate[FilterOptions[Graphics,opts]], Frame -> True,
        PlotRange -> {{n1, n2 + 1}, {lower, upper}}]
  ]

End[]

EndPackage[]