using Plots using SpecialFunctions # Use the PyPlot interface pyplot() # plotlyjs() # Initialize color C(g::ColorGradient) = RGB[g[z] for z = range(-1, stop=1, length=101)] g = :plasma_r # Draw a traveling plane-wave if !isfile("1-traveling-wave.gif") xrange = range(-2*pi, stop=2*pi, length=1001) wavefn(t,x) = exp(2im*(x - 2*t)) travlingwave = @animate for t = range(0,stop=pi,length=91) p1 = plot(xrange, real.(wavefn.(t,xrange)), label=""); p2 = plot(xrange, imag.(wavefn.(t,xrange)), label=""); p3 = plot(xrange, abs.(wavefn.(t,xrange)), label="", ylims=(0,1.5), lc=cgrad(g), line_z=(real.(wavefn.(t,xrange))), colorbar=:none, clims=(-1,1)); plot(p1, p2, p3, layout=(3,1)); end gif(travlingwave, "1-traveling-wave.gif", fps=15) end # Draw erfc # erfc is defined as 2 / sqrt(\pi) * \int_z^\infty e^{-z^2} dz # Here we use that (1 - im) / sqrt(2) squares to -im, to take care of the signs in the Schrodinger kernel if !isfile("2-errorfcplot.png") xrange = range(-20, stop=20, length=1001) errorfcplotr = plot(xrange, real.(erfc.( (1 - im)/ sqrt(2) * xrange)/2), label="", ylims=(-0.25,1.25)); errorfcploti = plot(xrange, imag.(erfc.( (1 - im)/sqrt(2) * xrange)/2), label="", ylims=(-0.75, 0.75)); errorfcplotm = plot(xrange, abs.(erfc.( (1-im)/sqrt(2) * xrange)/2), label="", lc=cgrad(g), line_z=(cos.(angle.(erfc.((1-im)/sqrt(2)*xrange)))), ylims=(-0.25,1.25), colorbar=:none, clims=(-1,1)); errorfcplot = plot(errorfcplotr, errorfcploti, errorfcplotm, layout=(3,1)); png(errorfcplot, "2-errorfcplot.png") end # Compute the solution for half-filled waves, using the formula (valid when t > 0) # \Phi_k(t,x) = e^{-i \pi/4} e^{ik(x - kt)} / \sqrt{\pi} \cdot # \int_{ k \sqrt{t} - x/\sqrt{4t}}^\infty e^{iz^2} dz # = 1/2 e^{ik(x-kt)} erfc( (1-im)/sqrt(2) * [ k sqrt(t) - x / sqrt(4t) ] ) if !isfile("3-simulation.gif") solution(k,t,x) = 1 / 2 * exp( im*k*(x - k*t )) * erfc( (1 - im)/sqrt(2) * ( k * sqrt(t) - x / sqrt(4*t))) xrange = range(-40 , stop = 40, length = 801) trange = range(-2, stop=25, length=271) # The use of max(s,0.0001) is a bit of a hack; I don't want to # worry about division by 0, and also this way it allows me to # easily have the animation "pause" on the initial data for # a second or so; this helps with the visualization with doing # looping. solnplot = @animate for s = trange t = max(s, 0.0001) soln0abs = abs.( solution.(0,t,xrange)) soln0arg = cos.(angle.(solution.(0,t,xrange))) soln0plot = plot(xrange, soln0abs, label="", lc=cgrad(g), line_z = soln0arg, colorbar=:none, ylims = (-0.25,1.25), clims=(-1,1), title_location=:top, title="time = $t"); soln2pabs = abs.( solution.(2,t,xrange)) soln2parg = cos.(angle.(solution.(2,t,xrange))) soln2pplot = plot(xrange, soln2pabs, label="", lc=cgrad(g), line_z = soln2parg, colorbar=:none, ylims = (-0.25,1.25), clims=(-1,1)); soln2nabs = abs.( solution.(-2,t,xrange)) soln2narg = cos.(angle.(solution.(-2,t,xrange))) soln2nplot = plot(xrange, soln2nabs, label="", lc=cgrad(g), line_z = soln2narg, colorbar=:none, ylims = (-0.25,1.25), clims=(-1,1)); plot(soln0plot, soln2pplot, soln2nplot, layout=(3,1)); soln0abs = 0 soln0arg = 0 soln2pabs = 0 soln2parg = 0 soln2nabs = 0 soln2narg = 0 GC.gc() end gif( solnplot, "3-simulation.gif", fps=15) end # Now, do the same thing, but for much larger times if !isfile("4-simulation.gif") solution(k,t,x) = 1 / 2 * exp( im*k*(x - k*t )) * erfc( (1 - im)/sqrt(2) * ( k * sqrt(t) - x / sqrt(4*t))) xrange = range(-40 , stop = 40, length = 801) trange = range(0, stop=25) solnplot = @animate for s = trange t = 1.5^s soln0abs = abs.( solution.(0,t,xrange)) soln0arg = cos.(angle.( solution.(0,t,xrange))) soln0plot = plot(xrange, soln0abs, label="", lc=cgrad(g,[-1,1]), line_z = soln0arg, colorbar=:none, ylims = (-0.25,1.25), clims=(-1,1), title_location=:top, title="time = $t"); soln2pabs = abs.( solution.(2,t,xrange)) soln2parg = cos.(angle.( solution.(2,t,xrange))) soln2pplot = plot(xrange, soln2pabs, label="", lc=cgrad(g,[-1,1]), line_z = soln2parg, colorbar=:none, ylims = (-0.25,1.25), clims=(-1,1)); soln2nabs = abs.( solution.(-2,t,xrange)) soln2narg = cos.(angle.( solution.(-2,t,xrange))) soln2nplot = plot(xrange, soln2nabs, label="", lc=cgrad(g,[-1,1]), line_z = soln2narg, colorbar=:none, ylims = (-0.25,1.25), clims=(-1,1)); plot(soln0plot, soln2pplot, soln2nplot, layout=(3,1)); soln0abs = 0 soln0arg = 0 soln2pabs = 0 soln2parg = 0 soln2nabs = 0 soln2narg = 0 GC.gc() end gif( solnplot, "4-simulation.gif", fps=2) end