## Quantum Harmonic Oscillator

Some trajectories of a harmonic oscillator (a ball attached to a spring) in classical mechanics (A-B) and quantum mechanics (C-H). In quantum mechanics, the position of the ball is represented by a wave (called the wavefunction), with real part shown in blue and imaginary part in red. Some of the trajectories, such as C,D,E,F, are standing waves (or "stationary states"). Each standing-wave frequency is proportional to a possible energy level of the oscillator. This "energy quantization" does not occur in classical physics, where the oscillator can have any energy.

(* Source code written in Mathematica 6.0 by Steve Byrnes, Feb. 2011. This source code is public domain. *)
(* Shows classical and quantum trajectory animations for a harmonic potential. Assume m=w=hbar=1. *)
ClearAll["Global`*"]
(*** Wavefunctions of the energy eigenstates ***)
psi[n_, x_] := (2^n*n!)^(-1/2)*Pi^(-1/4)*Exp[-x^2/2]*HermiteH[n, x];
energy[n_] := n + 1/2;
psit[n_, x_, t_] := psi[n, x] Exp[-I*energy[n]*t];
(*** A random time-dependent state ***)
SeedRandom;
CoefList = Table[Random[]*Exp[2 Pi I Random[]], {n, 0, 4}];
CoefList = CoefList/Norm[CoefList];
Randpsi[x_, t_] := Sum[CoefList[[n + 1]]*psit[n, x, t], {n, 0, 4}];
(*** A coherent state (or "Glauber state") ***)
CoherentState[b_, x_, t_] := Exp[-Abs[b]^2/2] Sum[b^n*(n!)^(-1/2)*psit[n, x, t], {n, 0, 15}];
(*** Make the classical plots...a red ball anchored to the origin by a gray spring. ***)
classical1[t_, max_] := ListPlot[{{max Cos[t], 0}}, PlotStyle -> Directive[Red, AbsolutePointSize]];
zigzag[x_] := Abs[(x + 0.25) - Round[x + 0.25]] - .25;
spring[x_, left_, right_] := (.9 zigzag[3 (x - left)/(right - left)])/(1 + Abs[right - left]);
classical2[t_, max_] := Plot[spring[x, -5, max Cos[t]], {x, -5, max Cos[t]}, PlotStyle -> Directive[Gray, Thick]];
classical3 = ListPlot[{{-5, 0}}, PlotStyle -> Directive[Black, AbsolutePointSize]];
classical[t_, max_, label_] := Show[classical2[t, max], classical1[t, max], classical3,
PlotRange -> {{-5, 5}, {-1, 1}}, Ticks -> None, Axes -> {False, True}, PlotLabel -> label, AxesOrigin -> {0, 0}];
(*** Put all the plots together ***)
SetOptions[Plot, {PlotRange -> {-1, 1}, Ticks -> None, PlotStyle -> {Directive[Thick, Blue], Directive[Thick, Pink]}}];
MakeFrame[t_] := GraphicsGrid[
{{classical[t + 2, 1.5, "A"], classical[t, 3, "B"]},
{Plot[{Re[psit[0, x, t]], Im[psit[0, x, t]]}, {x, -5, 5}, PlotLabel -> "C"],
Plot[{Re[psit[1, x, t]], Im[psit[1, x, t]]}, {x, -5, 5}, PlotLabel -> "D"]},
{Plot[{Re[psit[2, x, t]], Im[psit[2, x, t]]}, {x, -5, 5}, PlotLabel -> "E"],
Plot[{Re[psit[3, x, t]], Im[psit[3, x, t]]}, {x, -5, 5}, PlotLabel -> "F"]},
{Plot[{Re[Randpsi[x, t]], Im[Randpsi[x, t]]}, {x, -5, 5}, PlotLabel -> "G"],
Plot[{Re[CoherentState[1, x, t]], Im[CoherentState[1, x, t]]}, {x, -5, 5}, PlotLabel -> "H"]}
}, Frame -> All, ImageSize -> 300];
output = Table[MakeFrame[t], {t, 0, 4 Pi*96/97, 4 Pi/97}];
SetDirectory["C:\\Users\\Steve\\Desktop"]
Export["test.gif", output]