{ "cells": [ { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Saved animation to \n", "│ fn = /Users/sandvik/Courses/py502/Quantum/Examples/a.gif\n", "└ @ Plots /Users/sandvik/.julia/packages/Plots/vVVub/src/animation.jl:104\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"/Users/sandvik/Courses/py502/Quantum/Examples/a.gif\")" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ " using FFTW\n", " using Plots\n", "\n", "# Any potential for x between -ll and +ll can be coded here\n", "\n", " function potential(x::Float64,ll::Float64)::Float64\n", " v::Float64=0\n", " v=0.\n", " if abs(abs(x)-ll) < 1. # Walls of width 1 at both ends of the box\n", " v=500. # A large value; effectively a hard-wall box\n", " end\n", " return v\n", " end \n", "\n", "# Initializes the potential and kinetic evolution operators\n", "# - these are diagonal matrices and the dia elements are stored as vectors\n", " function operators(dt::Float64,nn::Int,ll::Float64)\n", " vdt=Vector{ComplexF64}(undef,nn) # potential-energy evolution operator\n", " kdt=Vector{ComplexF64}(undef,nn) # kinetic-energy evolution operator\n", " dx::Float64=2*ll/nn\n", " dk::Float64=2*pi/(2*ll)\n", " for i=-div(nn,2)+1:div(nn,2)\n", " x::Float64=i*dx\n", " k::Float64=i*dk\n", " i>0 ? j=i : j=i+nn # The elements are stored in the\n", " vdt[j]=exp(-dt*im*potential(x,ll)) # order used in the FT\n", " kdt[j]=exp(-dt*im*0.5*k^2)\n", " end\n", " return vdt,kdt\n", " end\n", "\n", "# Initializes a wave packet in momentum space\n", "# Inverse Fourier trandform to real space at the end and normalizes\n", "\n", " function initpacket(nn::Int,ll::Float64,a0::Float64,k0::Float64)\n", " psi=Vector{ComplexF64}(undef,nn)\n", " dk=pi/ll\n", " for i=-div(nn,2)+1:div(nn,2)\n", " k=i*dk\n", " i>0 ? j=i : j=i+nn\n", " psi[j]=exp(-0.5*(a0*(k-k0))^2)\n", " end\n", " psi=ifft(psi)\n", " r::Float64=sum(abs2.(psi))\n", " psi.*=1/sqrt(r*2*ll/nn)\n", " return psi\n", " end\n", "\n", "# Initializes a wave packet in real space (centered at x=x0)\n", "\n", " function xinitpacket(nn::Int,ll::Float64,a0::Float64,k0::Float64,x0::Float64)\n", " psi=Vector{ComplexF64}(undef,nn)\n", " dx=2*ll/nn\n", " for i=-div(nn,2)+1:div(nn,2)\n", " x=i*dx\n", " i>0 ? j=i : j=i+nn\n", " psi[j]=exp(-((x-x0)/a0)^2+im*k0*(x-x0))\n", " end\n", " r::Float64=sum(abs2.(psi))\n", " psi.*=1/sqrt(r*2*ll/nn)\n", " return psi\n", " end\n", "\n", "# Creates vectors x and y, containing the x grid points xi and psi^2(xi)\n", "# Orders the elements in the way needed for later graphical output\n", "\n", "function orderedvectors(nn::Int,ll::Float64,psi)\n", " x=Vector{Float64}(undef,nn)\n", " y=Vector{Float64}(undef,nn)\n", " dx::Float64=2*ll/nn\n", " for i=-div(nn,2)+1:div(nn,2)\n", " i>0 ? j=i : j=i+nn\n", " x[i+div(nn,2)]=i*dx\n", " y[i+div(nn,2)]=abs2(psi[j])\n", " end\n", " return x,y\n", " end\n", "\n", "# Split-operator evolution of the wave function psi by one time step dt\n", "\n", " function evolvestep(nn::Int,vdt,kdt,psi)\n", " psi .*= vdt\n", " psi=fft(psi)\n", " psi .*= kdt\n", " psi=ifft(psi)\n", " return psi\n", " end\n", "\n", "# Does the complete time evolution in bt steps\n", "# Stores the result after every pt steps in vectors x and y (=psi^2)\n", "# - y is appended each time with the next propagated psi^2\n", "\n", " function psiframes(nn::Int,ll::Float64,dt::Float64,nt::Int,pt::Int,psi)\n", " vdt,kdt=operators(dt,nn,ll)\n", " x,y=orderedvectors(nn,ll,psi)\n", " for i=1:nt\n", " psi=evolvestep(nn::Int,vdt,kdt,psi)\n", " if mod(i,pt)==0 \n", " x,z=orderedvectors(nn,ll,psi)\n", " append!(y,z) \n", " end\n", " end\n", " return x,y\n", " end\n", "\n", " nn=1024 # Number of x points in the box\n", " ll=20. # box is located between -ll and ll (L=2*LL)\n", " a0=1. # sets the width of the initial wave packet in momentum space\n", " k0=20. # mean momentum of the initial wave packet\n", " x0=0. # Location of the initial wave pucket if initialized in x-space\n", " dt=0.001 # time step\n", " nt=8000 # Number of time steps\n", " pt=40 # Every pt time step will be included in the animation\n", "\n", " psi=initpacket(nn,ll,a0,k0) # Initialize wave function in momentum space\n", "# psi=xinitpacket(nn,ll,a0,k0,x0) # Initialize in real space\n", "\n", " x,y=psiframes(nn,ll,dt,nt,pt,psi)\n", "\n", "# x contains the real-space grid points of the box\n", "# y is a vector containing all evolved wave functions (psi^2)\n", "# - the segment y[nn*i+1:nn*(i+1)] contains the result after (i-1)*pt steps\n", "# An animation is created using the @animate macro and the gif(0 function)\n", "\n", " anim = @animate for i=0:div(nt,pt)\n", " plot(x,y[nn*i+1:nn*(i+1)])\n", "# plot(x,y[nn*i+1:nn*(i+1)],ylim=(0,maximum(y))\n", " end\n", " gif(anim,\"a.gif\",fps=10) # converts \"anim\" to a gif animation\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.6.1", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.1" } }, "nbformat": 4, "nbformat_minor": 4 }