Fourier Series Movie

1.0 Set constants.

The constants are set for an electron (m=0.511 MeV/(c^2)) in an infinite square well 1.0 angstroms in width with the initial wave packet a square wave located between 'a' and 'b'.

a0=1.0;
x0=0.3;
x1=0.5;
d0=x1-x0;
hbarc=1970 ;(* eV-angstroms *)
hbar = 6.5821*(10^(-16)) (* eV-s *);
mass= 0.55*(10^6)(* electron mass in eV. *);
E0 = (((hbarc^2)*(π^2))/(2*mass*(a0^2))) (* energy scale in eV. *);

2.0 Square wave definition.

Define the initial wave function and create a plot of it for later use.

wf[x_] := Which[x<x0,0,
			x0≤x≤x1,(1/(Sqrt[d0])),
			x>x1,0
]
wf1 = Plot[wf[x]*Conjugate[wf[x]],{x,0,a0},
PlotRange->{{0,1},{-0.2,8.0}},
Frame->True,
FrameLabel->{"x (angstroms)","Probability Density","Electron in a Box"},
PlotStyle->RGBColor[1,0,0],
FrameStyle->Black,
BaseStyle->FontSize->16,Exclusions->None
]

3.0 Make Movies.

Make the movie. The radial frequency, ω, is 

                ω= ((Subscript[E, n])/ℏ) = (n^2)5.2902×(10^16)(s^(-1))
               
for the constants entered in part 1.0. To get reasonably-sized numbers enter the time in units of (10^(-16)) s.

(* parameters for the plot and making the movie. *)
nterms=50;
nsteps=200;
xmin=0.0;
xmax=2.0;
tmax=0.060;
tstep=(tmax/nsteps);
tconvert=(10^(-16));
Energy[n_]:= (n^2)*E0
psi[x_,time_] := Sum[(Sqrt[((2*a0)/((n^2)*(π^2)*d0))])*(Cos[((n*π*x0)/a0)]-Cos[((n*π*x1)/a0)])*(Sqrt[(2/a0)])*Sin[((n*π*x)/a0)]*Exp[-I*((Energy[n])/hbar)*time*tconvert],{n,1,nterms}]
Animate[
wf2=Plot[Evaluate[Re[psi[x,time]*Conjugate[psi[x,time]]]],{x,xmin,xmax},PlotRange->{{0,1},{0,8.0}},Frame->True,FrameLabel->{"x","Probability Density",StringForm["time=``\[Cross](10^(-16))s.",NumberForm[time,{2,4}]]," "},
FrameStyle->Black,
BaseStyle->FontSize->16];
Show[wf2,wf1,Graphics[Text["a=1.0 Å",{0.7,6.0},{-1,0}]],Graphics[Text["electron",{0.7,7.0},{-1,0}]]],
{time,0,tmax,tstep}
]

4.0 Make individual frames at specific times.

time=0;
wf3=Plot[Evaluate[Re[psi[x,0]*Conjugate[psi[x,0]]]],{x,xmin,xmax},PlotRange->{{0,1},{0,8.0}},Frame->True,FrameLabel->{"x","Probability Density",StringForm["time=``\[Cross](10^(-16))s.",NumberForm[time,{2,4}]]," "},
FrameStyle->Black,
BaseStyle->FontSize->16];
wf4=Show[wf3,wf1]

Export["/run/media/gilfoyle/KINGSTON/309f17/examples/TimeDevelopment/dt1.pdf",wf4]

time=0.004;
wf5=Plot[Evaluate[Re[psi[x,time]*Conjugate[psi[x,time]]]],{x,xmin,xmax},PlotRange->{{0,1},{0,8.0}},Frame->True,FrameLabel->{"x","Probability Density",StringForm["time=``\[Cross](10^(-16))s.",NumberForm[time,{2,4}]]," "},
FrameStyle->Black,
BaseStyle->FontSize->16];
wf6=Show[wf5,wf1]

Export["/run/media/gilfoyle/KINGSTON/309f17/examples/TimeDevelopment/dt2.pdf",wf6]

time=0.008;
wf7=Plot[Evaluate[Re[psi[x,time]*Conjugate[psi[x,time]]]],{x,xmin,xmax},PlotRange->{{0,1},{0,8.0}},Frame->True,FrameLabel->{"x","Probability Density",StringForm["time=``\[Cross](10^(-16))s.",NumberForm[time,{2,4}]]," "},
FrameStyle->Black,
BaseStyle->FontSize->16];
wf8=Show[wf7,wf1]

Export["/run/media/gilfoyle/KINGSTON/309f17/examples/TimeDevelopment/dt3.pdf",wf8]

time=0.012;
wf9=Plot[Evaluate[Re[psi[x,time]*Conjugate[psi[x,time]]]],{x,xmin,xmax},PlotRange->{{0,1},{0,8.0}},Frame->True,FrameLabel->{"x","Probability Density",StringForm["time=``\[Cross](10^(-16))s.",NumberForm[time,{2,4}]]," "},
FrameStyle->Black,
BaseStyle->FontSize->16];
wf10=Show[wf9,wf1]

Export["/run/media/gilfoyle/KINGSTON/309f17/examples/TimeDevelopment/dt4.pdf",wf10]