{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "<center><h2>Quantum Tunneling</h2></center>\n",
    "<p>\n",
    "Two isolated potential wells $|L\\rangle$ and $|R\\rangle$    have eigenenergies:\n",
    "$$\n",
    "H_0|L\\rangle  = E_0 |L\\rangle, \\ \\ \\ H_0 |R\\rangle  =  E_0 |R\\rangle. \n",
    "$$\n",
    "This program examines the time-dependent motion of a wave packet when an interaction $H_{int}$ connects the wells. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true,
    "scrolled": false
   },
   "outputs": [],
   "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",
    "# QuantumTunnel.ipynb: Quantum tunneling of a wave packet between two potential wells. \n",
    "\n",
    "% matplotlib notebook\n",
    "\n",
    "from vpython import *\n",
    "import numpy as np\n",
    "\n",
    "dx  =  0.05;    dx2  =  dx*dx;  k0  =  3.;  dt  =  dx2/15.0;  xmax  =  6.0\n",
    "nmax = 240              # 240 =  2*xmax/dx\n",
    "\n",
    "g = display(width = 500,height = 500,title = 'Wave packet in two wells',background = color.white,\\\n",
    "  foreground = color.black)\n",
    "PlotObj =  curve( color = color.red, radius = 0.02) # to plot packet\n",
    "potential = curve(color = color.yellow, radius = 0.02)  # to plot potential\n",
    "v = np.zeros((nmax+1),float)     # potential\n",
    "psr = np.zeros((nmax+1),float) # real psi\n",
    "psi = np.zeros((nmax+1),float) # imaginary psi\n",
    "i = 0                        \n",
    "for x in np.arange (-xmax,0,dx):\n",
    "    i = i+1\n",
    "    v[i] = 20*(x+2.5)**2    # left hand well\n",
    "i = 119                     # Right hand well\n",
    "for x in np.arange(0,xmax,dx):   # RH well\n",
    "    v[i] = 20*(x-2.5)**2\n",
    "    i = i+1\n",
    "prob = np.zeros((nmax),float)       # initial psi\n",
    "k = 0\n",
    "for xs in np.arange(-xmax,xmax,dx):\n",
    "   psr[k]  =  np.exp(-5.5*((xs+4.5))**2) *np. cos(k0*xs)    # Re Psi\n",
    "   psi[k]  =  np.exp(-5.5*((xs+4.5))**2) * np.sin(k0*xs)    # Im Psi\n",
    "   k += 1 \n",
    "prob = psr*psr + psi*psi                 # probability\n",
    "j = 0     \n",
    "for x in np.arange(-xmax, xmax,dx):    \n",
    "     PlotObj.append(pos = vec(x,5*prob[j]-1,0))   # 5*probability lowered 2\n",
    "     potential.append(pos = vec(x,0.03*v[j]-5,0))  # scaled potential to plot\n",
    "     j  =  j+1\n",
    "for t in range(0,3500):                      \n",
    "   rate(20) \n",
    "   psr[1:-1] = psr[1:-1] - (dt/dx2)*(psi[2:]+psi[:-2]-2*psi[1:-1]) +dt*v[1:-1]*psi[1:-1]\n",
    "   psi[1:-1] = psi[1:-1] + (dt/dx2)*(psr[2:]+psr[:-2]-2*psr[1:-1]) -dt*v[1:-1]*psr[1:-1]\n",
    "   PlotObj.clear()\n",
    "   if t%100 == 0: print(t)\n",
    "   x = -xmax\n",
    "   prob = psr**2 + psi**2\n",
    "   for k in range(0,nmax):\n",
    "      PlotObj.append(pos = vec(x, 5*prob[k]-1,0))  # plot the wave packet with time\n",
    "      x += dx\n",
    "   PlotObj.visible = True    \n",
    "print(\"finished\")"
   ]
  }
 ],
 "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
}