{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"collapsed": true,
"deletable": true,
"editable": true
},
"source": [
" Quantum Scattering from Three Disks, Vpython
\n",
"\n",
"The time dependent Schroedinger\n",
"equation \n",
"$$\n",
"i\\hbar \\frac {\\partial \\psi(x,z)}{\\partial t} = -\\frac{\\hbar^2}{2m} \\Bigl ( \\frac{\\partial^2 \\psi(x,z)}{\\partial x^2} + \\frac{\\partial^2 \\psi(x,z)}{\\partial z^2}\\Bigr )\n",
"+ V(x,z) \\psi \n",
"$$\n",
"is solved for scattering of the Gaussian wave packet\n",
"$$\n",
"\\psi(x,z) = \\exp(i(k_0 x +k_1 z))\\exp(-A(x-x_i)^2-A(z-z_i)^2) \n",
"$$ from three disks fixed to the vertices of an equilateral triangle. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"deletable": true,
"editable": true
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glow.min\");}else{element.textContent = ' ';}"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glowcomm\");}else{element.textContent = ' ';}"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/jquery-ui.custom.min\");}else{element.textContent = ' ';}"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glow.min\"], function(){console.log(\"GLOW LOADED\");});}else{element.textContent = ' ';}"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glowcomm\"], function(){console.log(\"GLOWCOMM LOADED\");});}else{element.textContent = ' ';}"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/javascript": [
"if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/jquery-ui.custom.min\"], function(){console.log(\"JQUERY LOADED\");});}else{element.textContent = ' ';}"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"0\n",
"10\n",
"20\n",
"30\n",
"40\n",
"50\n",
"60\n",
"70\n",
"80\n",
"90\n",
"100\n",
"110\n",
"120\n",
"130\n",
"140\n",
"finito\n"
]
}
],
"source": [
"\"\"\" From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\"\n",
" by RH Landau, MJ Paez, and CC Bordeianu (deceased).\n",
" Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, \n",
" C Bordeianu, Univ Bucharest, 2021. \n",
" Please respect copyright & acknowledge our work.\"\"\"\n",
"\n",
"# 3QMdisksVP.pynb: same as 3QMdisks, but now with Vpython\n",
"\n",
"% matplotlib notebook\n",
"\n",
"from vpython import *\n",
"import numpy as np\n",
"\n",
"scene = canvas(width = 500, height = 500,range = 120,background = color.white,foreground = color.black)\n",
"table = curve(pos = ([vec(-100,0,-100),vec(100,0,-100),vec(100,0,100),\n",
" vec(-100,0,100),vec(-100,0,-100)]))\n",
"x1 = 51 # int (90.*sqrt(3)/2 -30)\n",
"R = 24 # radius of disks\n",
"circ1 = ring(pos = vec(-30,0,45), radius = R,axis = vec(0,1,0),color = color.blue)\n",
"circ2 = ring(pos = vec(-30,0,-45),radius = R,axis = vec(0,1,0),color = color.blue)\n",
"circ3 = ring(pos = vec(x1,0,0),radius = R,axis = vec(0,1,0),color = color.blue)\n",
"scene.forward = vec(0,-1,-1) # angle of scene viewing\n",
"xmax = 100 \n",
"nmax = 101 \n",
"V = np.zeros((nmax,nmax),float) # potential\n",
"dx = 0.1 \n",
"dx2 = dx*dx \n",
"k0 = 20. # wave packet x wave vector\n",
"k1 = 0. # wave packet z vector\n",
"dt = 0.002 \n",
"fc = dt/dx2 \n",
"Re = np.zeros((nmax,nmax),float) # Real psi\n",
"I = np.zeros((nmax,nmax),float) # Imaginary psi\n",
"\n",
"def pot(xa,za): # potential, disks centered at xa,za\n",
" for z in range (za-R,za+R+1): # R : disk radius, \n",
" for x in range(xa-R,xa+R+1): # ends of x axis of disk\n",
" if sqrt((x-xa)**2+(z-za)**2) <= R: # defines potential\n",
" i = int(50./100.*z+50) # convert from x z to i j\n",
" j = int(50./100.*x+50) \n",
" V[i,j] = 5. # equivalent to infinite potential value\n",
"\n",
"def potential():\n",
" pot(-30,45) # center of each ring\n",
" pot(-30,-45)\n",
" pot(x1,0) \n",
"\n",
"def initial(xin,zin): # initial position of wave packet (xin,zin) \n",
" for i in np.arange(0,nmax): z = 200./100.*i-100 \n",
" for j in arange(0, nmax):\n",
" x = 200./100.*j-100 # convert index to coordinate \n",
" Re[i,j] = exp(-0.01*(x-xin)**2- 0.01*(z-zin)**2)*cos(k0*x+k1*z)\n",
" I[i,j] = exp(-0.01*(x-xin)**2-0.01*(z-zin)**2)*sin(k0*x+k1*z)\n",
" \n",
"potential() # potential due to disks\n",
"xin = -50 # initial position of wave packet, -100 <= z, x <= 100\n",
"zin = 24 \n",
"initial(xin,zin)\n",
"\n",
"def plotinitial():\n",
" for i in range(1,nmax-1):\n",
" zp = 200.*i/nmax-100\n",
" for j in range(1,nmax-1): \n",
" if V[i,j]!= 0: # psi = 0 when V != 0\n",
" Re[i,j] = 0\n",
" I[i,j] = 0\n",
" yy = 40*(Re[i,j ]**2 +I[i, j ]**2)\n",
" if yy>.01: # to avoid printing long lines\n",
" xx = 200.*j/nmax-100. # coordinates of neighboring points \n",
" xm1 = 200.*(j-1)/nmax-100.\n",
" yym1 = 40* (Re[i,j-1]**2 +I[i, j-1 ]**2)\n",
" zz = zp \n",
" curve(pos = [vec(xm1,yym1,zz),vec(xx,yy,zz)],color = color.red) \n",
"plotinitial() # plot a short segment\n",
"\n",
"for t in range(0,150): # plot every 10 times\n",
" if t%10 == 0: print(t) \n",
" I[1:-1,1:-1] = I[1:-1,1:-1]+fc*(Re[2: ,1:-1]+Re[:-2 ,1:-1]-4*Re[1:-1,1:-1] + \\\n",
" Re[1:-1,2: ]+Re[1:-1, :-2])+ V[1:-1,1:-1]*dt*Re[1:-1,1:-1]\n",
" Re[1:-1, 1:-1] = Re[1:-1,1:-1]-fc*(I[2: ,1:-1]+I[ :-2,1:-1]-4*I[1:-1,1:-1] \\\n",
" +I[1:-1,2: ]+I[1:-1, :-2])+ V[1:-1,1:-1]*dt*I[1:-1,1:-1]\n",
" for i in range(1, nmax-1):\n",
" zp = 200.*i/nmax-100\n",
" for j in range(1,nmax-1):# wavefunction 0 when potential is diff of 0\n",
" if V[i,j] != 0: # wave function 0 causes reflections\n",
" Re[i,j] = 0\n",
" I[i,j] = 0\n",
" yy = 40*(Re[i,j ]**2 +I[i, j]**2)\n",
" xx = 200.*j/nmax-100. \n",
" xm1 = 200.*(j-1)/nmax-100.\n",
" yym1 = 40*(Re[i,j-1 ]**2 +I[i, j-1 ]**2)\n",
" zz = zp\n",
" if yy>0.1:\n",
" if t%10 == 0: curve(pos =[vec(xm1,yym1,zz),vec(xx,yy,zz)],color = color.red)\n",
"print(\"finito\") "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}