In a previous post I provided a C program for making animations of wavepackets colliding with potential steps. C is a compiled language, so the compiled executable runs very fast, fast enough that the animation will run in real time despite each of the incident, reflected and transmitted waves being 200 component Fourier sums of plane waves calculated at each of 1000 points at each times step, and I had 1000 time steps.

The various “modern” computer languages such as python or julia are nowhere near as fast, but are a bit easier to learn and program in. As far as raw speed is concerned

which is not disparage python, it’s built for comfort, not for speed.

No matter what language you code this program in, computing every wave wave Fourier component at every value should not be done within the time loop. Do it once outside of the time loop, and then just plot time slices within the time loop.

The code below in julia, easily translated to python by removing some “end”s and adding some indentation, draws a png snapshot of the system at each time step, then the png’s can be assembled into an mp4 stream with ffmpeg. To assemble say 500 png’s:

ffmpeg -r 60 -f image2 -s 800x400 -start_number 0 -i pic%04d.png -vframes 200 -vcodec libx264 -crf 25 pulse.mp4
# to make an html5 compatible mp4...
ffmpeg -i pulse.mp4 -c:v libx264 -crf 23 -c:a aac -pix_fmt yuv420p -movflags faststart  pulse2.mp4

Start a julia session, load the code below (download from http://abadonna.physics.wisc.edu/download/stepwave2) and run it for 500 timesteps in an otherwise empty directory. Then quite julia, run the ffmpeg command to make the mp4, then you can remove all of the png’s.

using Plots

stepwave=function(time)
L=20.0
dx=0.04
M=2.0 
C=3.0
sigma=2.0
k0=6.0
V0=4.0
#time = 1000
dk=2.0/(100.0*sigma)
factor=sqrt(2.0)*dk*sigma/sqrt(pi);
Amp=[exp(-4.0*(n-100.0)^2*dk^2/2.0) for n in 1:200]
iexpt=zeros(Complex,200,time);
iexpx=zeros(Complex,200,1000);
Amp=zeros(Complex,200,1);
psi=zeros(Complex,1000,1);
Amp=[exp(-4.0*(n-100.0)^2*dk^2/2.0) for n in 1:200]
X=[n*dx for n in 0:999]

for n in 1:200
k=k0+dk*(n-100)
for T in 1:time
t = T *0.02
iexpt[n,T]=exp(-(0.0+1.0im)*abs(k)*C*t)
end
end

for n in 1:200
k=k0+dk*(n-100)
if 0.5*k*k/M>=V0
kp=sqrt(k*k-2.0*M*V0);
else
kp=(0.0+1.0im)*sqrt(2.0*M*V0-k*k);
end
for p in 1:499
x = p *dx;	     
iexpx[n,p]=exp((0.0+1.0im)*k*x)+((k-kp)/(k+kp))*exp((0.0-1.0im)*k*(x-2.0*L));
end
for p in 500:1000
x = p *dx;
iexpx[n,p]=(2.0*k/(k+kp))*exp((0.0+1.0im)*(kp*x+k*L-kp*L))
end
end

for t in 1:time
for n in 1:1000
psi[n]=factor*sum([iexpx[p,n]*iexpt[p,t]*Amp[p] for p in 1:200])
end
#Y=[abs(psi[n,1]) for n in 1:1000]
plot(X,abs.(psi)[:,1],xlims=(0,40.0), ylims=(0,4.0))
plot!(X,stepV.(X),xlims=(0,40.0), ylims=(0,4.0),linecolor = :red)
savefig("pic"*string(t,pad=4)*".png")
end
end

stepV=function(x)
if x<=20.0
return 0.0
else
return 0.5
end
end

And now my friends, Stepwave, the movie, …