{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Investigating the $\\alpha%$ decay using the split-step algorithm\n",
"\n",
"A simple model for $\\alpha$ decay involves a particle in a square well incident on a Coulomb barrier with an energy such that it must tunnel through this barrier.\n",
"\n",
"A gaussian wavepacket can be used in order to model the behaviour of a particle and, to avoid the periodic boundary conditions, a near-hard wall is set up near x = 0.\n",
"\n",
"For an explaination of the split-step method see the [previous notebook](split_step_schrodinger.ipynb)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setting up potential and wavepacket functions\n",
"\n",
"Starting with a Gaussian wavepacket incident on the potential barrier which is a Coulomb barrier which starts a third of the way along the x-domain."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"import pycav.pde as pde\n",
"import pycav.display as display\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.animation as anim\n",
"\n",
"def oneD_gaussian(x,mean,std,k0):\n",
" return np.exp(-((x-mean)**2)/(4*std**2)+ 1j*x*k0)/(2*np.pi*std**2)**0.25\n",
"\n",
"def V(x):\n",
" V_x = np.zeros_like(x)\n",
" mid = int(len(x)/3)\n",
" V_x[mid:] = 250/(x[mid:]-x[0])\n",
" V_x[:50] = 10**2\n",
" return V_x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We must define the spatial domain as well as the time step size and number of steps. Then the initial wavefunction is passed to the split step algorithm."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"N_x = 2**11\n",
"dx = 0.1\n",
"x = dx * (np.arange(N_x) - 0.5 * N_x)\n",
"\n",
"dt = 0.01\n",
"N_t = 8000\n",
"\n",
"p0 = 2.5\n",
"d = np.sqrt(N_t*dt/2.)\n",
"\n",
"psi_0 = oneD_gaussian(x,x[0]+5*d,d,p0)\n",
"\n",
"psi_x,psi_k,k = pde.split_step_schrodinger(psi_0, dx, dt, V, N_t, x_0 = x[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Energy of $\\alpha$ particle and max potential barrier height\n",
"\n",
"The ratio of the particle energy within the well compared to the max height of the Coulomb barrier. As long as this is less than unity then the decay is classically forbidden and must proceed via tunneling."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.8525\n"
]
}
],
"source": [
"E_particle = p0**2/2\n",
"E_barrier = V(x)[50:].max()\n",
"print(E_particle/E_barrier)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The time evolution can then be animated showing the expected behaviour. The momentum space wavefunction can also be visualised (here shown in the second subplot)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"real_psi = np.real(psi_x)\n",
"imag_psi = np.imag(psi_x)\n",
"absl_psi = np.absolute(psi_x)\n",
"abs_psik = np.absolute(psi_k)\n",
"\n",
"fig = plt.figure(figsize = (10,10))\n",
"ax1 = plt.subplot(211)\n",
"line1_R = ax1.plot(x,real_psi[:,0],'b')[0]\n",
"line1_I = ax1.plot(x,imag_psi[:,0],'r')[0]\n",
"line1_A = ax1.plot(x,absl_psi[:,0],'k')[0]\n",
"line_V = ax1.plot(x,0.12*(V(x)-p0**2/2),'k',alpha=0.5)[0]\n",
"ax1.set_ylim((real_psi.min(),real_psi.max()))\n",
"ax1.set_xlim((x.min(),x.max()))\n",
"\n",
"ax2 = plt.subplot(212)\n",
"line2_A = ax2.plot(k,abs_psik[:,1],'k')[0]\n",
"ax2.set_ylim((abs_psik.min(),abs_psik.max()))\n",
"ax2.set_xlim((-10,10))\n",
"\n",
"def nextframe(arg):\n",
" line1_R.set_data(x,real_psi[:,10*arg])\n",
" line1_I.set_data(x,imag_psi[:,10*arg])\n",
" line1_A.set_data(x,absl_psi[:,10*arg])\n",
" line2_A.set_data(k,abs_psik[:,10*arg])\n",
" \n",
"animate = anim.FuncAnimation(fig, nextframe, frames = int(N_t/10), interval = 50, repeat = False)\n",
"animate = display.create_animation(animate,temp = True)\n",
"display.display_animation(animate)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Probability of decay\n",
"\n",
"By integrating the probability density either side of the point of the Coulomb potential where the potential is equal to the energy of particle in the well, the probability of decay as a function of time can be visualised"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3oAAAFwCAYAAADuVwk3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmYXFW57/Hv24QYIBICMkggDGGSScYQBaUBhQAHgsJB\n4oAMF3IRhIMTigM5Kh4iehDkOqCIqCgOgKAyC42IkIBACDETgyFMiShhCpCh3/vHqqY7bSfd6XSn\nqqu/n+fZT+3aQ9VbbDbkl7X2WpGZSJIkSZLqR0O1C5AkSZIk9SyDniRJkiTVGYOeJEmSJNUZg54k\nSZIk1RmDniRJkiTVGYOeJEmSJNWZToNeRFwaEXMj4qHlHHNRRMyKiAcjYpc220dHxPSImBkRZ7XZ\nPjQibo6IGRFxU0QMWfmfIkmSJEmCrrXoXQYctKydEXEwMCIztwbGAd+rbG8ALq6cuwMwNiK2q5z2\nWeDWzNwWuA34XLd/gSRJkiRpKZ0Gvcz8M/D8cg4ZA/ykcuxEYEhEbAiMBGZl5uzMXARcWTm25ZzL\nK+uXA0d0r3xJkiRJUns98YzeMGBOm/dPVrYtazvAhpk5FyAznwU26IE6JEmSJEn0zmAs0Y1zsser\nkCRJkqR+akAPfMZTwKZt3m9S2TYQGN7BdoBnI2LDzJwbERsB85b14RFhCJQkSZLUr2XmCjWodTXo\nBctuqbsOOBX4ZUSMAuZXAtxzwFYRsRnwDHAMMLbNOccBE4CPAtcu78szzXq1bPz48YwfP77aZfR7\nzc2wcCG89lrr8vrrcMEF4/nwh8fzwgswfz5Lvb7wAjz/PPzjH63LP/8JgwfDBhvARhvB5pvDFlu0\nLiNGwMYbQ3Sn7V4d8h6qbV6f2uc1qn1eo9rm9al90Y0/eHUa9CLi50AjsF5EPAGcQ2mty8y8JDOv\nj4hDIuIR4BXgeMrOJRFxGnAzpYvopZk5rfKxE4BfRcQJwGzg6BWuXNJSGhpg0KCytLXRRrDPPl3/\nnObm1vD3zDPw+ONlueWW8vrIIyVQ7rxzWXbaqbzuuGMJiJIkSaq+ToNeZn6wC8ectoztNwLbdrD9\nX8B7ulKgpFWroQHWW68s220H++3378fMmwdTpsBDD8Hdd8Mll8Df/gZvfevS4W/nnUsL4Gqrrfrf\nIUmS1J/1xDN66ucaGxurXYKWozeuzwYbwAEHlKXF4sWlta8lAP70p2V93jzYc094xztg773Lss46\nPV5Sn+Y9VNu8PrXPa1T7vEa1zetTn6LWn3+LiKz1GiUt2/z5MHEi/OUv8Oc/w6RJpaVw9OiyjBpl\ni58kSdLyRMQKD8Zi0JO0Si1cWLp73ngj/OEPpcXv8MPhyCPhPe8x9EmSJLVn0JPU5zz2GFxzDfzy\nl/DUU3DssXD88bDNNtWuTJIkqTYY9CT1aVOnwmWXwc9+VgZxOf54OPpoWHvtalcmSZJUPQY9SXVh\n0SK44YYS+m6/HY44Aj7+cdh992pXJkmStOoZ9CTVnXnzSuD7zndg2DA4/XQ46igY4JjBkiSpnzDo\nSapbixfD734HF1wAc+bAJz4BJ5wAa61V7cokSZJ6l0FPUr9wzz1w/vlw551wyilw2mmw/vrVrkqS\nJKl3dCfoNfRWMZLUW0aNgquuKvPyPfssbLstnHoq/P3v1a5MkiSpNhj0JPVZ22wD3/8+/O1vMGRI\nGazlpJPg8cerXZkkSVJ1GfQk9XkbbQRf+xrMnAkbbgh77AEnnwxPPFHtyiRJkqrDoCepbqy3Hnz1\nqzBrFrzlLbDrrnDmmfDcc9WuTJIkadUy6EmqO+uuW1r4/va3MlrndtvBV74CCxZUuzJJkqRVw6An\nqW5tuCF8+9swcSJMnVoGbbniCmhurnZlkiRJvcvpFST1G3fdVbpyNjSU+fje8Y5qVyRJktQ5p1eQ\npOXYe+8yB9+pp8J//iccfzz861/VrkqSJKnndSnoRcToiJgeETMj4qwO9q8TEVdHxOSIuCcitm+z\n74yImFJZTm+z/e0RcXdEPBARkyJij575SZK0bA0N8JGPwPTp8OY3w447ljn5JEmS6kmnXTcjogGY\nCRwAPA3cCxyTmdPbHPN14KXM/EpEbAv8v8x8T0TsAPwC2BNYDNwIjMvMxyLiJuCbmXlzRBwMfCYz\n9+vg++26KanX3HUXnHhiCXwXX1ymapAkSaolvdV1cyQwKzNnZ+Yi4EpgTLtjtgduA8jMGcDmEbE+\n8DZgYma+nplLgDuA91fOaQaGVNbXAZ5akcIlqSfsvTc8+GAZqOXtb4fLLwf/bkmSJPV1XQl6w4A5\nbd4/WdnW1mQqAS4iRgLDgU2Ah4F3RcTQiFgTOATYtHLOmcA3IuIJ4OvA57r7IyRpZQwaBOeeCzfe\nCN/6Fhx8MMyeXe2qJEmSum9AD33OecCFEXE/MAV4AFiSmdMjYgJwC/Byy/bKOacAZ2TmbyPiKOBH\nwHs7+vDx48e/sd7Y2EhjY2MPlS1JrXbdFSZNgm98A/bYA8aPh1NOKc/1SZIkrSpNTU00NTWt1Gd0\n5Rm9UcD4zBxdef9ZIDNzwnLOeRzYKTNfbrf9XGBOZn4vIuZn5jpt9r2QmUM6+Cyf0ZO0yk2fXp7d\na2iAH/6wdO2UJEmqht56Ru9eYKuI2CwiBgLHANe1++IhEbF6Zf0k4I6WkFd5Vo+IGA68D7iictpT\nEbFvZd8BlAFfJKkmbLcd3HknHH10eY7voot8dk+SJPUdXZowPSJGAxdSguGlmXleRIyjtOxdUmn1\nu5wywMpU4MTMfKFy7p+AdYFFwJmZ2VTZ/k7gImA14DXgY5n5QAffbYuepKp65BH40Idg3XXhsssc\nmVOSJK1a3WnR61LQqyaDnqRasGgR/Pd/w6WXlrA3enS1K5IkSf2FQU+Setkdd8CHPwxjx5aROldf\nvdoVSZKketdbz+hJkir23Rfuvx+mToV3vQv+/vdqVyRJkvTvDHqStILWXx9+9zv4z/+EvfaCa6+t\ndkWSJElLs+umJK2Ee+6BY46BI4+E886zK6ckSep5dt2UpFVs1Cj4619hxgx497thzpxqVyRJkmTQ\nk6SVtt56cN11MGYM7Lkn3HJLtSuSJEn9nV03JakH3X57mXPvlFPg85+HBv86TZIkrSSnV5CkGvD0\n0+W5vbXWgp/9rLT4SZIkdZfP6ElSDdh4Y/jjH2GnnWC33WDixGpXJEmS+huDniT1gtVXh69/HS68\nEA47DL79bbBzgiRJWlXsuilJvezRR8uce1ttBT/8Iay9drUrkiRJfYldNyWpBo0YAX/5CwwdWkbl\nnDKl2hVJkqR6Z9CTpFVg0CD4/vfhi1+E/feHyy+vdkWSJKme2XVTklaxqVPhqKNg773Ls3trrFHt\niiRJUi2z66Yk9QE77AD33gsLFsA73gGzZlW7IkmSVG8MepJUBYMHwxVXwLhxpWXvqquqXZEkSaon\nXQp6ETE6IqZHxMyIOKuD/etExNURMTki7omI7dvsOyMiplSW09ud9/GImFbZd97K/xxJ6jsi4JRT\n4Prr4VOfgjPPhIULq12VJEmqB50GvYhoAC4GDgJ2AMZGxHbtDjsbeCAz3w58FLiocu4OwInAHsAu\nwGERsWVlXyNwGLBTZu4EfKMnfpAk9TV77AH33w+PPAKNjfDEE9WuSJIk9XVdadEbCczKzNmZuQi4\nEhjT7pjtgdsAMnMGsHlErA+8DZiYma9n5hLgDuD9lXNOAc7LzMWV855b6V8jSX3U0KFw7bVwxBEl\n+P3619WuSJIk9WVdCXrDgDlt3j9Z2dbWZCoBLiJGAsOBTYCHgXdFxNCIWBM4BNi0cs42wLsrXT1v\nj4g9uv8zJKnva2iAz3wG/vAHOPtsOPFEePnlalclSZL6op4ajOU8YGhE3A+cCjwALMnM6cAE4Bbg\n+pbtlXMGAEMzcxTwGeBXPVSLJPVpe+5ZunI2N8Nuu8Ff/1rtiiRJUl8zoAvHPEVpoWuxSWXbGzLz\nJeCElvcR8TjwWGXfZcBlle3n0to6+CRwdeWYeyOiOSLWy8x/ti9g/Pjxb6w3NjbS2NjYhbIlqe96\n85vhssvgl7+Egw8ug7V86lOl1U+SJNW3pqYmmpqaVuozOp0wPSJWA2YABwDPAJOAsZk5rc0xQ4AF\nmbkoIk4C9s7M4yr71s/Mf0TEcOBGYFRmvhgRJwPDMvOciNgGuCUzN+vg+50wXVK/Nns2fPjD8KY3\nweWXw7D2neclSVJd65UJ0yuDqJwG3AxMBa7MzGkRMa4S1qAMuvJwREyjjM55RpuPuCoiHgauBT6W\nmS9Wtl8GbBkRU4CfA8euSOGS1F9sthncfjvsuy/svnsZtEWSJGl5Om3RqzZb9CSp1d13w4c+BAcd\nBN/8Jqy5ZrUrkiRJva1XWvQkSbXjHe+ABx8so3Huvjvcd1+1K5IkSbXIoCdJfczaa8NPfwpf+hIc\nckh5Xbiw2lVJkqRaYtCTpD5q7NjSuvfXv8LIkWVdkiQJDHqS1KdtvDH8/vfwX/8FBx5YWvdef73a\nVUmSpGoz6ElSHxcBxx1XWvQmTy7P7k2aVO2qJElSNRn0JKlObLwx/Pa38IUvwOGHwyc/WQZtkSRJ\n/Y9BT5LqSAQccwxMmQLz5sGOO8If/lDtqiRJ0qrmPHqSVMduuQVOOQV22QW+9S3YZJNqVyRJklaU\n8+hJkpby3veW1r3tty9h7/zznYpBkqT+wBY9SeonZs2C00+Hxx6DCy4oc/BJkqTa150WPYOeJPUj\nmXD99fCJT8CIEfC//wvbbVftqiRJ0vLYdVOStFwRcOihpTvne94D73oXnHkmPP98tSuTJEk9yaAn\nSf3QwIGlVW/qVFiwoLTqfec7sGhRtSuTJEk9wa6bkiQmTy7z7s2ZA1/7Grz//aX1T5IkVZ/P6EmS\nui0Tbr4ZPvtZWH11+OpXy6idBj5JkqrLoCdJWmnNzfCb38AXvwgbbwxf/nJ5lk+SJFVHrw3GEhGj\nI2J6RMyMiLM62L9ORFwdEZMj4p6I2L7NvjMiYkplOb2Dcz8ZEc0Rse6KFC5J6h0NDXD00eX5vWOP\nhY9+tAzcctdd1a5MkiR1VadBLyIagIuBg4AdgLER0X4w7rOBBzLz7cBHgYsq5+4AnAjsAewC/EdE\nbNnmszcB3gvMXvmfIknqSQMGwPHHw4wZMHYsfPjDsN9+cOutpZunJEmqXV1p0RsJzMrM2Zm5CLgS\nGNPumO2B2wAycwaweUSsD7wNmJiZr2fmEuBPwPvbnHcB8OmV/A2SpF60+upw4ollwvUTTiiTro8a\nBb/+NSxeXO3qJElSR7oS9IYBc9q8f7Kyra3JVAJcRIwEhgObAA8D74qIoRGxJnAIsGnluMOBOZk5\nZaV+gSRplRgwAD7yEXj44TJgy0UXwVZblUnXX3ih2tVJkqS2emoevfOAoRFxP3Aq8ACwJDOnAxOA\nW4DrW7ZHxBqU7p7ntPkMx3WTpD6goQHe9z6480741a/g3nthyy3LxOuPP17t6iRJEsCALhzzFKWF\nrsUmlW1vyMyXgBNa3kfE48BjlX2XAZdVtp9LaR0cAWwOTI6IqHzmXyNiZGbOa1/A+PHj31hvbGyk\nsbGxC2VLknrbyJHwi1+U+fe+/W3Yc0945zvh5JNh9OjSCihJklZMU1MTTU1NK/UZnU6vEBGrATOA\nA4BngEnA2Myc1uaYIcCCzFwUEScBe2fmcZV962fmPyJiOHAjMCozX2z3HY8Du2Xm8x18v9MrSFIf\n8fLLpZXvkkvgySfLs30nnACbbVbtyiRJ6rt6ZXqFyiAqpwE3A1OBKzNzWkSMi4iTK4e9DXg4IqZR\nRuc8o81HXBURDwPXAh9rH/Javga7bkpSnzd4cAl299wDN9wAzz8Pu+0GBx8M11wDixZVu0JJkvoH\nJ0yXJPWqV18tE7BfcglMnw5jxsCRR8IBB8DAgdWuTpKk2tedFj2DniRplfn73+Hqq+Gqq2DaNDj0\nUDjqKDjwQFhjjWpXJ0lSbTLoSZL6jKefLt05f/MbeOABOOig0tJ3yCGlC6gkSSoMepKkPmnePLj2\n2tLSd/fdsP/+JfQddhgMGVLt6iRJqi6DniSpz3v+ebjuuhL6mppgn31K4DvooDJfnyRJ/Y1BT5JU\nV156Cf7whzKC5003lS6d73lPafHbbz9Yf/1qVyhJUu8z6EmS6lYmPPww3HIL3H473HknrLce7LVX\nmbh9r71g111h0KBqVypJUs8y6EmS+o3m5jJdw6RJMHFiWaZPh+23L6GvJQBusw00dDprrCRJtcug\nJ0nq1xYsKCN4TpzYGgCffx723LM1+I0cCRtuCLFC/7uUJKl6DHqSJLUzdy7ce29rq99995VuoNts\nA1tvXV5b1rfeGtZeu9oVS5K0NIOeJEmdyITnnoNZs2DmzLK0rM+aVaZzaBsAW0LgiBE+/ydJqg6D\nniRJK6G5uUzk3jb8taz//e/w1reW0LfZZmX9rW+FjTZqfd1oI8OgJKnnGfQkSeolixeXsDdzJsyZ\nA888A88+u/Tr3Lmw1lpLh7/2gXCDDcpooeuuC296U7V/lSSpLzDoSZJURc3NZfCXZ57pOAg+8wzM\nmwf/+ldZBg4sgW+ddWDo0PK6zjrlOcG11irLmmv++/qyXgcNcpAZSapHBj1JkvqITHj55RL45s9f\nennhhTKC6CuvlKVlfXnbFiyAhQthjTW6Fgrf9KYSNFdfvby2X2/7fvXVYcAAWG21f39d2W0NDSWc\ntiySpH9n0JMkqR9bsqQEvo5CYPuAuHBh67Jo0bLft6wvWVK6r3b02tVtHe3LLC2hLSJaw1/bENh+\n2/L2rci29uGy7fvl7evp97X0XSu7r9Y/b3n72v470n69J/cNGND6ly0DB3a+PmgQvPnNZVl77fK6\nxhr+5Uh/YtCTJEl9UubSS3Pz0q8dbVvevq5ua19DR+u9/b6Wvmtl99X653X2Xe3/HWm/3lP7Fi9u\n/cuU119f+rWj9VdfhZdeWnpZtAgGD24NgO2XlkC49trwlreU54M32KDMI7rBBuVc9R0GPUmSJKkf\nWLSodP9+8cV/D4FtlxdeKFPKzJtXlrlzy9LQsHT4e+tbYfjwMqrw8OFlGTasdN1W9fVa0IuI0cC3\ngAbg0syc0G7/OsCPgBHAq8AJmfm3yr4zgP9TOfSHmXlhZfvXgcOA14FHgeMz88UOvtugJ0mSJPWQ\nzNKFe+7c1gD41FPwxBNLL88+W0Jgy5yi227bumy2WXnOVqtGrwS9iGgAZgIHAE8D9wLHZOb0Nsd8\nHXgpM78SEdsC/y8z3xMROwC/APYEFgM3AuMy87GIeA9wW2Y2R8R5QGbm5zr4foOeJEmStIotXlym\nk5k1C2bMKMvMmeV13jwYMQJ23hl22aV12WCDalddn7oT9AZ04ZiRwKzMnF35kiuBMcD0NsdsD/wP\nQGbOiIjNI2J94G3AxMx8vXLuHcD7gW9k5q1tzr8HOHJFCpckSZLUewYMgC22KMuBBy69b8GCEvoe\neggefBD+53/K6xprwMiRsPfe8M53wu67l8FktOp1JegNA+a0ef8kJfy1NZkS4O6KiJHAcGAT4GHg\nqxExlNJF8xBKi2B7JwBXrljpkiRJkqphzTVbW/GOPbZsy4TZs2HiRLjrLjjjDJg2Dd7+9tbgt/fe\ntvqtKl0Jel1xHnBhRNwPTAEeAJZk5vSImADcArzcsr3tiRHxeWBRZv58WR8+fvz4N9YbGxtpbGzs\nobIlSZIk9YQI2HzzsnzgA2Xbyy/DvfeW4PeDH8AJJ5T9Bx8Mo0fDO97hgC8daWpqoqmpaaU+oyvP\n6I0Cxmfm6Mr7z1Kep5uwnHMeB3bKzJfbbT8XmJOZ36u8Pw44Cdi/pXtnB5/lM3qSJElSHVi8GO65\nB268EW64AR57DPbfvzX4bbJJtSusTb01GMtqwAzKYCzPAJOAsZk5rc0xQ4AFmbkoIk4C9s7M4yr7\n1s/Mf0TEcMpgLKMy88XKSJ7fBN6dmf9czvcb9CRJkqQ6NHcu3HRTCX4331ymeRg9Gg4/vHT1dGTP\norenV7iQ1ukVzouIcZSWvUsqrX6XA83AVODEzHyhcu6fgHWBRcCZmdlU2T4LGAi0hLx7MvNjHXy3\nQU+SJEmqc0uWwH33wfXXw29/W0b2POIIOPJIaGwsg8P0V06YLkmSJKkuzJoFV18Nv/41PP10GfTl\n+OPLPH79jUFPkiRJUt3529/gssvgpz+FLbcsge8DH4C11652ZauGQU+SJElS3Vq0qDzPd9llcNtt\ncNhhJfQ1NkJDQ7Wr6z0GPUmSJEn9wj/+AVdcUULfiy/CRz8Kxx1Xpm+oNwY9SZIkSf1KJjzwQAl8\nv/gF7LxzaeU78sgysXs9MOhJkiRJ6rdeew2uu66EvokTSyvf6afDFltUu7KV052gV8c9WSVJkiT1\nJ4MGwdFHl8nYH3wQVl8d9twTjjoK/vKX0vrXX9iiJ0mSJKluvfxyaeG78EJYbz0488zSrXP11atd\nWdfZdVOSJEmSOrBkCfzud3DBBfD44/Dxj8NJJ8E661S7ss7ZdVOSJEmSOrDaanDEEXDHHXDNNTB5\ncpmT7/TT4dFHq11dzzPoSZIkSepXdt8dfvYzmDIF1loLRo2C972vPMdXL+y6KUmSJKlfe+UVuPxy\nOP/8MkLnF79YJmGPFeos2Xt8Rk+SJEmSumnRojIJ+9e+BhtsUALfgQdWP/AZ9CRJkiRpJS1ZAr/8\nJZx7LgweDF/4AvzHf1Qv8Bn0JEmSJKmHNDfD1VfDV79aQt4XvlCe5WtYxSOdGPQkSZIkqYdllqkZ\nvvIVePVV+Pzny8Tsq622ar7foCdJkiRJvSQTbrqpBL7nnoOzz4YPfQgGDOjd7+21efQiYnRETI+I\nmRFxVgf714mIqyNickTcExHbt9l3RkRMqSynt9k+NCJujogZEXFTRAxZkcIlSZIkaVWKgNGj4c9/\nhu9+F378Y9hmG/jBD2DhwmpXt7ROg15ENAAXAwcBOwBjI2K7doedDTyQmW8HPgpcVDl3B+BEYA9g\nF+CwiNiycs5ngVszc1vgNuBzK/9zJEmSJKl3RcD++8Ptt8NPfgK/+Q1stRV85zvw2mvVrq7oSove\nSGBWZs7OzEXAlcCYdsdsTwlrZOYMYPOIWB94GzAxM1/PzCXAHcD7K+eMAS6vrF8OHLFSv0SSJEmS\nVrF99indOX/9a7jhBhgxAi64ABYsqG5dXQl6w4A5bd4/WdnW1mQqAS4iRgLDgU2Ah4F3Vbpprgkc\nAmxaOWfDzJwLkJnPAht090dIkiRJUjXttVcZsOX3vy9dO7fcEiZMgJdeqk49PTUw6HnA0Ii4HzgV\neABYkpnTgQnALcD1LduX8RmOuCJJkiSpT9t1V7jqKrj1VnjwwdLC95WvwPz5q7aOrowP8xSlha7F\nJpVtb8jMl4ATWt5HxOPAY5V9lwGXVbafS2vr4LMRsWFmzo2IjYB5yypg/Pjxb6w3NjbS2NjYhbIl\nSZIkqTp23BF+8QuYMQO+9rXyDN/HPgb/9V+w7rrLP7epqYmmpqaV+v5Op1eIiNWAGcABwDPAJGBs\nZk5rc8wQYEFmLoqIk4C9M/O4yr71M/MfETEcuBEYlZkvRsQE4F+ZOaEykufQzPxsB9/v9AqSJEmS\n+rRHH4XzzisTsJ90EnziE7BBFx9e65XpFSqDqJwG3AxMBa7MzGkRMS4iTq4c9jbg4YiYRhmd84w2\nH3FVRDwMXAt8LDNfrGyfALw3IlpC5HkrUrgkSZIk9RUjRpRpGB54oDy3t912Jew980zvfJ8TpkuS\nJEnSKvbUU3D++WV6hrFj4ayzYPjwjo/ttQnTJUmSJEk9Z9gw+Na3YNo0GDy4DOJy8snw2GM98/kG\nPUmSJEmqkg03LNMwzJxZ1keOhOOOK+9XhkFPkiRJkqpsvfXKNAyPPFLm4Nt7b/jgB2Hq1O59ns/o\nSZIkSVKNefFF+M534IILYN68FX9Gz6AnSZIkSTXqlVdg8GCDniRJkiTVFUfdlCRJkiQZ9CRJkiSp\n3hj0JEmSJKnOGPQkSZIkqc4Y9CRJkiSpzhj0JEmSJKnOGPQkSZIkqc4Y9CRJkiSpzhj0JEmSJKnO\nGPQkSZIkqc50KehFxOiImB4RMyPirA72rxMRV0fE5Ii4JyK2b7PvzIh4OCIeiogrImJgZfvbI+Lu\niHggIiZFxB4997MkSZIkqf/qNOhFRANwMXAQsAMwNiK2a3fY2cADmfl24KPARZVzNwY+DuyWmTsD\nA4BjKud8HTgnM3cFzgHOX/mfo2poamqqdglaDq9P7fMa1TavT+3zGtU+r1Ft8/rUp6606I0EZmXm\n7MxcBFwJjGl3zPbAbQCZOQPYPCLWr+xbDVgrIgYAawJPV7Y3A0Mq6+sAT3X7V6iq/I9DbfP61D6v\nUW3z+tQ+r1Ht8xrVNq9PfepK0BsGzGnz/snKtrYmA+8HiIiRwHBgk8x8Gvgm8AQlyM3PzFsr55wJ\nfCMinqC07n2uuz9CkiRJktSqpwZjOQ8YGhH3A6cCDwBLImIdSuvfZsDGwOCI+GDlnFOAMzJzOCX0\n/aiHapEkSZKkfi0yc/kHRIwCxmfm6Mr7zwKZmROWc85jwM7AaOCgzDypsv0jwF6ZeVpEzM/Mddqc\n80JmDungs5ZfoCRJkiTVucyMFTl+QBeOuRfYKiI2A56hDKYytu0BETEEWJCZiyLiJOBPmflypVvm\nqIgYBLwOHABMqpz2VETsm5l3RMQBwMye+EGSJEmS1N91GvQyc0lEnAbcTOnqeWlmTouIcWV3XgK8\nDbg8IpqBqcCJlXMnRcRvKF05F1Vef1D56JOAiyJiNeA14OSe/WmSJEmS1D912nVTkiRJktS39NRg\nLD2us0naVX0R8feImNwy6X216xFExKURMTciHmqzbWhE3BwRMyLipkpXa1XJMq7RORHxZETcX1lG\nV7PG/iwiNomI2yJiakRMiYjTK9u9j2pAB9fn45Xt3kM1IiLeFBETK382mBIR51S2ew/ViOVcI++j\nGhIRDZXrcF3l/QrfQzXZoleZpH0m5Zm+pynPCR6TmdOrWpiWUhl0Z/fMfL7ataiIiH2Al4GfZObO\nlW0TgH/do0SjAAAgAElEQVRm5tcrf2kyNDM/W806+7NlXKNzgJcy83+rWpyIiI2AjTLzwYgYDPyV\nMnr08XgfVd1yrs8H8B6qGRGxZmYuqDyecxdwOnAk3kM1YxnX6GC8j2pGRJwJ7A6snZmHd+fPc7Xa\noteVSdpVfUHt/jvUL2Xmn4H2wXsMcHll/XLgiFValJayjGsE5X5SlWXms5n5YGX9ZWAasAneRzVh\nGdenZW5f76EakZkLKqtvoowHkXgP1ZRlXCPwPqoJEbEJcAjwwzabV/geqtU/pHdlknZVXwK3RMS9\nldFWVZs2yMy5UP6QBGxQ5XrUsdMi4sGI+KFdmmpDRGwO7ALcA2zofVRb2lyfiZVN3kM1otLl7AHg\nWeCWzLwX76GasoxrBN5HteIC4NO0BnDoxj1Uq0FPfcPembkb5W8cTq10SVPtq73+2voOsGVm7kL5\nn67dZqqs0i3wN8AZlZaj9veN91EVdXB9vIdqSGY2Z+aulNbwkRGxA95DNaWDa7Q93kc1ISIOBeZW\nei8sr4W103uoVoPeU8DwNu83qWxTDcnMZyqv/wCuoXS5Ve2ZGxEbwhvPt8yrcj1qJzP/ka0PTP8A\n2LOa9fR3ETGAEiJ+mpnXVjZ7H9WIjq6P91BtyswXgSZgNN5DNantNfI+qhl7A4dXxsL4BbB/RPwU\neHZF76FaDXpvTNIeEQMpk7RfV+Wa1EZErFn5G1UiYi3gQODh6lalimDpvwG6Djiusv5R4Nr2J2iV\nW+oaVf6D3eL9eC9V24+Av2XmhW22eR/Vjn+7Pt5DtSMi3tLS5S8i1gDeS3mW0nuoRizjGk33PqoN\nmXl2Zg7PzC0pGei2zPwI8DtW8B6qyVE3oUyvAFxI6yTt51W5JLUREVtQWvGS8hDvFV6j6ouInwON\nwHrAXOAc4LfAr4FNgdnA0Zk5v1o19nfLuEb7UZ41agb+Doxr6YevVSsi9gb+BEyh/PctgbOBScCv\n8D6qquVcnw/iPVQTImInykARDZXll5l5bkSsi/dQTVjONfoJ3kc1JSL2BT5ZGXVzhe+hmg16kiRJ\nkqTuqdWum5IkSZKkbjLoSZIkSVKdMehJkiRJUp0x6EmSJElSnTHoSZIkSVKdMehJkiRJUp0x6EmS\nJElSnTHoSZIkSVKdMehJkiRJUp0x6EmSJElSnTHoSZIkSVKdMehJkiRJUp0x6EmSJElSnelS0IuI\n0RExPSJmRsRZyzjmooiYFREPRsSubbZfGhFzI+Khdsd/PSKmVY6/KiLWXrmfIkmSJEmCLgS9iGgA\nLgYOAnYAxkbEdu2OORgYkZlbA+OA77bZfVnl3PZuBnbIzF2AWcDnuvULJEmSJElL6UqL3khgVmbO\nzsxFwJXAmHbHjAF+ApCZE4EhEbFh5f2fgefbf2hm3pqZzZW39wCbdO8nSJIkSZLa6krQGwbMafP+\nycq25R3zVAfHLM8JwA0rcLwkSZIkaRmqPhhLRHweWJSZP692LZIkSZJUDwZ04ZingOFt3m9S2db+\nmE07OebfRMRxwCHA/ss5JrtQoyRJkiTVrcyMFTm+Ky169wJbRcRmETEQOAa4rt0x1wHHAkTEKGB+\nZs5tsz8qS+uGiNHAp4HDM/P15RWQmS41vJxzzjlVr8HF69OXF69RbS9en9pfvEa1v3iNanvx+tT+\n0h2dBr3MXAKcRhklcypwZWZOi4hxEXFy5Zjrgccj4hHg+8DH2gS6nwN/AbaJiCci4vjKrm8Dg4Fb\nIuL+iPhOt36BJEmSJGkpXem6SWbeCGzbbtv3270/bRnnfnAZ27fuYo2SJEmSpBVQ9cFY1Pc1NjZW\nuwQth9en9nmNapvXp/Z5jWqf16i2eX3qU3S3z+eqEhFZ6zVKkiRJUm+JCLIXBmORJEmSJPUhBj1J\nkiRJqjMGPUmSJEmqMwY9SZIkSaozBj1JkiRJqjMGPUmSJEmqMwY9SZIkSaozBj1JkiRJqjMGPUmS\nJEmqMwY9SZIkSaozBj1JkiRJqjMGPUmSJEmqRY8+Cief3K1TDXqSJEmSVEumT4ePfAT22gs22qhb\nH2HQkyRJkqRaMGUKfOADsO++8La3lRa9L3+5Wx/VpaAXEaMjYnpEzIyIs5ZxzEURMSsiHoyIXdts\nvzQi5kbEQ+2OHxoRN0fEjIi4KSKGdOsXSJIkSVJfdv/98L73wYEHwp57loB39tkwpPsRqdOgFxEN\nwMXAQcAOwNiI2K7dMQcDIzJza2Ac8N02uy+rnNveZ4FbM3Nb4Dbgc936BZIkSZLUF02cCP/xH3D4\n4bDffvDYY/CpT8HgwSv90V1p0RsJzMrM2Zm5CLgSGNPumDHATwAycyIwJCI2rLz/M/B8B587Bri8\nsn45cMSKly9JkiRJfcydd5bWu6OPhkMPhUcegdNPhzXW6LGv6ErQGwbMafP+ycq25R3zVAfHtLdB\nZs4FyMxngQ26UIskSZIk9T2ZcPvtpeXuuOPKs3izZsEpp8CgQT3+dQN6/BO7L6tdgCRJkiT1qEz4\n4x/hv/8b5s6FL3wBPvhBGNC7Uawrn/4UMLzN+00q29ofs2knx7Q3NyI2zMy5EbERMG9ZB44fP/6N\n9cbGRhobGzuvWpIkSZKqJRNuvbUEvOeeKwHvmGO6FPCamppoampaqa+PzOU3pEXEasAM4ADgGWAS\nMDYzp7U55hDg1Mw8NCJGAd/KzFFt9m8O/C4zd2qzbQLwr8ycUBnJc2hmfraD78/OapQkSZKkmpAJ\nN91UAt78+fDFL5Zumqut1u2PjAgyM1bonK6EqIgYDVxIeabv0sw8LyLGAZmZl1SOuRgYDbwCHJ+Z\n91e2/xxoBNYD5gLnZOZlEbEu8CtKS+Bs4OjMnN/Bdxv0JEmSJNW2TLj++jLv3SuvlIB31FErFfBa\n9FrQqyaDniRJkqSalQm//30JeK+9Bl/6Ehx5JDR0acryLulO0KulwVgkSZIkqW9oboZrroGvfrWE\nvS99CY44okcD3sow6EmSJElSVy1ZAr/5DZx7LgwcWJ7FO+wwiBVqcOt1Bj1JkiRJ6syiRfDzn8PX\nvgbrrgv/8z9wyCE1F/BaGPQkSZIkaVleew1+/GOYMAG22AK++90y6XmNBrwWBj1JkiRJam/+fLjk\nErjwQthlF7jiCnjnO6tdVZfVxpOCkiRJklQLZs+GT3wCttwSHnqojKj5hz/0qZAHBj1JkiRJgvvv\nhw9+EHbbrcx9N3ky/OxnsOuu1a6sWwx6kiRJkvqn5uYyyfn++8OYMbD77vDYY3D++bDpptWubqX4\njJ4kSZKk/uX118szd9/8Zpki4VOfgqOPhtVXr3ZlPcagJ0mSJKl/+Ne/4Hvfg29/uwywctFFpTWv\nxkfQ7A67bkqSJEmqb489BqefDlttBY88ArfcAjfcAAccUJchD2zRkyRJklSPliyBm2+G738f7roL\nTjoJHn4YNt642pWtEgY9SZIkSfXjmWfgRz+CH/wA1l8fxo0ro2cOHlztylYpg54kSZKkvq25GW69\ntbTe3XZbGVjl6qvLVAn9lEFPkiRJUt80ZUpprfvFL2DddeH//l+47DJYe+1qV1Z1Bj1JkiRJfcec\nOfDzn5fpEZ5/Hj70IfjDH2CnnapdWU2JzKx2DcsVEVnrNUqSJEnqRc8/D7/5TQl3U6bAkUeWgPeu\nd0FD/U8kEBFk5goND9qlfyoRMToipkfEzIg4axnHXBQRsyLiwYjYpbNzI+LtEXF3RDwQEZMiYo8V\nKVySJElSHVuwAK66Ct7/fth8c7jpJviv/4Knn4ZLLoF99+0XIa+7Om3Ri4gGYCZwAPA0cC9wTGZO\nb3PMwcBpmXloROwFXJiZo5Z3bkTcBHwzM2+unP+ZzNyvg++3RU+SJEnqD+bOhd//Hq69FpqaYORI\n+OAHS9hbZ51qV1c13WnR68ozeiOBWZk5u/IlVwJjgOltjhkD/AQgMydGxJCI2BDYYjnnNgNDKuev\nAzy1IoVLkiRJ6uMWLoR77ikTmN9yC0yfDgcdBMccA5dfDkOHVrvCPqsrQW8YMKfN+ycp4a+zY4Z1\ncu6ZwE0R8U0ggHd2vWxJkiRJfU4mTJvWGuzuvBO23hre+14491zYZx9405uqXWVd6K1RN7vSrHgK\ncEZm/jYijgJ+BLy3owPHjx//xnpjYyONjY09UKIkSZKkXrV4MUyeDHffDX/5C/zpTzBgQAl2xx4L\nP/4xvOUt1a6y5jQ1NdHU1LRSn9GVZ/RGAeMzc3Tl/WeBzMwJbY75HnB7Zv6y8n46sC+l62aH50bE\n/Mxcp81nvJCZQ2jHZ/QkSZKkPiAT/v53eOABuO++Eu7uuw822wze+c6y7LMPjBgBsUKPm/V7vfWM\n3r3AVhGxGfAMcAwwtt0x1wGnAr+sBMP5mTk3Ip7r4NxjKuc8FRH7ZuYdEXEAZdAWSZIkSbUsswya\nMn166YY5bRpMnVoC3hprwG67leWss2CvvXzOrko6DXqZuSQiTgNupkzHcGlmTouIcWV3XpKZ10fE\nIRHxCPAKcPxyzm0ZxOUk4KKIWA14DTi5x3+dJEmSpBW3YEGZxuCpp8rrnDkwY0ZrsGtogLe9rXU5\n5BDYdVfYcMNqV64KJ0yXJEmS+oslS0prXEuAa/vadv3VV2HjjWHYsNbXbbdtDXbrr1/tX9KvdKfr\npkFPkiRJqnXNzWUqgtdeg9dfL6+vvgovvNC6zJ+/9Pv22+fNg3/8A9Zdd+kA1/Ladn3ddX2OroYY\n9CRJUvU0N5cR9touixa1rjc3l9aE5ublr7d/3/LngMz+t74s/XVfZ+csa2n596gr21fk2JX57Obm\n1sD2+utLr3f0ungxDBwIgwaV6QcGDSrLkCFlIvEhQ5a9tOxff33YaCNYffVl/3NUTeqtwVgkSVIt\nySx/s//yy/DKK62vr7xS/kC4rKXlD4zLWhYu/Pdw1na9s/fNzWXY9NVXL68ty+qrw2qrlaWhofW1\n/XpH+yJaF+if68vSX/d1dk77paFh2dvb/zu2vGNXdnv7bQ0NS4e2jl7brg8caAubVogtepIkVcPC\nhaUbVdvlX/9q7Wb14outywsvtK6/9FIJdA0NMHgwrLVWWVrW11ij9W/62//N/7KWlmMGDmwNaR2F\ntWW9b1lv+YOsJKlH2XVTkqRqyWx9BmbevDLYwfLWX3mldKPaYIMySt3665dnYtp2wVp77aWXN7+5\nvK61ll2vJKkfMehJktQbmpvh2WfLRMCzZ7cuTzxRtreEuEGDWoPbBhssf33oUFu/JEldYtCTJGll\nPPdcmfR36lR45BF49NHy+vjjpTVt883LstlmZRk+HN761tbwNmhQtX+BJKkOGfQkSeqK+fNLmHv4\n4dZg9/DDZbCSHXYoy9Zbw4gRZdlyyxL0JEmqAoOeJEltZZbulvfcA/ffDw89VALdiy/C9tvDjju2\nBrsddyzzR9mdUpJUYwx6kqT+bfFimDwZ7roL/vzn8trcDO94B+y+O+yySwl1w4eXESIlSeoDDHqS\npP7lpZdKa11LsJs0CTbdFPbeG/bZp7xuuaWtdJKkPs2gJ0mqb//8J9x2G9x5Zwl306fDbru1Brt3\nvrNMUSBJUh0x6EmS6suSJaWV7sYbyzJtGrz73WXZZ5/SHfNNb6p2lZIk9SqDniSp71uwAG65BX77\nW/j978sAKaNHw0EHlZY7g50kqZ8x6EmS+qaXXirB7uqrS9fMPfaAI46AMWPKwCmSJPVjBj1JUt+x\neDHcfDP89Kdw/fWw775w1FFw6KGw3nrVrk6SpJrRnaDXpbGlI2J0REyPiJkRcdYyjrkoImZFxIMR\nsUtXzo2Ij0fEtIiYEhHnrUjhkqQ+KBPuuw/OOAOGDYOvfKU8a/foo3DddXDssYY8SZJ6wIDODoiI\nBuBi4ADgaeDeiLg2M6e3OeZgYERmbh0RewHfA0Yt79yIaAQOA3bKzMUR8Zae/nGSpBoxezZccUVp\nvVu4ED784TIdwtZbV7sySZLqUqdBDxgJzMrM2QARcSUwBpje5pgxwE8AMnNiRAyJiA2BLZZz7inA\neZm5uHLecz3zkyRJNWHhQrjmGvjud+Hhh+Hoo+HSS8vk5c5rJ0lSr+pK0BsGzGnz/klK+OvsmGGd\nnLsN8O6I+BrwKvDpzLyv66VLkmrS00/DJZeUZbvt4LTT4PDDYeDAalcmSVK/0ZWg1x1d+avaAcDQ\nzBwVEXsCvwK27KV6JEm9KbN0xbz44jI1wjHHlNcddqh2ZZIk9UtdCXpPAW3Htt6ksq39MZt2cMzA\n5Zz7JHA1QGbeGxHNEbFeZv6zfQHjx49/Y72xsZHGxsYulC1J6nWLFsGvfgXf+Aa8+iqcemppyRsy\npNqVSZLUZzU1NdHU1LRSn9Hp9AoRsRowgzKgyjPAJGBsZk5rc8whwKmZeWhEjAK+VWmpW+a5ETEO\n2Dgzz4mIbYBbMnOzDr7f6RUkqda89BL88IfwrW/BllvCpz9dJjVv6NJgzpIkaQV0Z3qFTlv0MnNJ\nRJwG3EyZjuHSNkEtM/OSzLw+Ig6JiEeAV4Djl3du5aN/BPwoIqYArwPHrkjhkqQqeP55uOii0kVz\nv/3gqqvK5OaSJKmmOGG6JKlz//wnXHBBGUFzzBj43OecGkGSpFWk1yZMlyT1Uy+8AF/6EmyzDcyd\nWyY7/9GPDHmSJNU4g54k6d8tWAATJpRA9+STJeD94AewxRbVrkySJHWBQU+S1GrJkjLIyjbblHD3\npz+VFjwDniRJfUpvzaMnSepr/vxnOP10GDwYrrkG9tyz2hVJkqRuMuhJUn/35JPwmc+UoHf++XD0\n0RAr9Ly3JEmqMXbdlKT+6rXX4KtfhV12ga22gmnT4AMfMORJklQHbNGTpP4mE377W/jkJ2G33cqz\neJtvXu2qJElSDzLoSVJ/MnUqnHEGPPtsGXRl//2rXZEkSeoFdt2UpP7g+efLQCv77QdHHAEPPmjI\nkySpjhn0JKmeZcLll8N228GiRfC3v8Fpp8EAO3RIklTP/D+9JNWrxx+HcePguefghhvK83iSJKlf\nsEVPkurNkiVwwQVlHrz3vAcmTTLkSZLUz9iiJ0n1ZMoUOPFEWGstuPtu2HrralckSZKqwBY9SaoH\nCxfCOeeUAVZOOgn++EdDniRJ/ZgtepLU102ZAh/9KLz1rWU0zWHDql2RJEmqMlv0JKmvWrIEJkwo\nrXinngq//70hT5IkAV0MehExOiKmR8TMiDhrGcdcFBGzIuLBiNilq+dGxCcjojki1u3+z5CkfubR\nR+Hd74abboJ77y3P5UVUuypJklQjOg16EdEAXAwcBOwAjI2I7dodczAwIjO3BsYB3+vKuRGxCfBe\nYHaP/BpJqneZcOmlMGoUHHUU3HorbL55tauSJEk1pivP6I0EZmXmbICIuBIYA0xvc8wY4CcAmTkx\nIoZExIbAFp2cewHwaeC6HvgtklTf5s8vLXePPAK33w477ljtiiRJUo3qStfNYcCcNu+frGzryjHL\nPDciDgfmZOaUFaxZkvqf++4rc+FtvHGZF8+QJ0mSlqO3Rt1c7oMiEbEGcDal22aXzpGkfikTLr4Y\nvvxl+M534D//s9oVSZKkPqArQe8pYHib95tUtrU/ZtMOjhm4jHNHAJsDkyMiKtv/GhEjM3Ne+wLG\njx//xnpjYyONjY1dKFuS+rgXXoD/839KV82774attqp2RZIkaRVoamqiqalppT4jMnP5B0SsBswA\nDgCeASYBYzNzWptjDgFOzcxDI2IU8K3MHNWVcyvnPw7slpnPd/D92VmNklR37r8fjj4aDjwQ/vd/\nYdCgalckSZKqJCLIzBXqAdlpi15mLomI04CbKc/0XZqZ0yJiXNmdl2Tm9RFxSEQ8ArwCHL+8czv6\nGuy6KUmlq+Z3vwvnnAPf/jYcc0y1K5IkSX1Qpy161WaLnqR+48UX4eSTYdo0+PWvYZttql2RJEmq\nAd1p0evShOmSpF724IOwxx4wZAjcc48hT5IkrRSDniRVUyZ8//vw3veW7prf/z6ssUa1q5IkSX1c\nb02vIEnqzEsvwbhxMGUK3HknbLddtSuSJEl1whY9SaqGhx4qXTXXXBMmTjTkSZKkHmXQk6RVKRN+\n+EM44AD4/OfL+pprVrsqSZJUZ+y6KUmryssvwymnlDny7rgDtt++2hVJkqQ6ZYueJK0KU6fCnnvC\naqvBpEmGPEmS1KsMepLU2378Y2hshLPOKutrrVXlgiRJUr2z66Yk9ZYFC+DUU8u8eLffDjvuWO2K\nJElSP2GLniT1hmnTYORIWLwY7r3XkCdJklYpg54k9bSf/hTe/W4480z4yU9g8OBqVyRJkvoZu25K\nUk9ZsABOP71Mfv7HP8LOO1e7IkmS1E/ZoidJPWH6dNhrrxL27rvPkCdJkqrKoCdJK+uKK+Bd7yqt\neVdcAW9+c7UrkiRJ/ZxdNyWpu9p21bz1Vvj/7d19jFX1mcDx74MKvrSCbitaUZdqs7i0ytqutsXG\nWc0qoBW6tKxSu9omSrK+bDYbq9U2appttWl17baVFbUR0xbforBJjS/BCdVWJdVBdJn1DSi6Mnbt\nohaL4eXZP85BhmHmzlyY4ZyZ+X6Sm3PumfO7PJcfz3Cfe37n9zv22KojkiRJAryiJ0k7p+tQTYs8\nSZJUIxZ6ktSMTLj5ZodqSpKkWutToRcRUyKiPSJeiIjLejjnhxHxYkS0RcSk3tpGxPciYkV5/r0R\nsf+uvx1JGkBvvgkzZ8JNN8GSJXD++RBRdVSSJEk76LXQi4gRwI+A04CJwNkRMaHLOVOBIzPzY8Ac\nYG4f2j4ETMzMScCLwDf65R1J0kBYvBgmTYLx4+GJJ+Doo6uOSJIkqUd9uaJ3PPBiZq7OzI3AAmB6\nl3OmA/MBMvNJYHREjG3UNjMfycwtZfsngHG7/G4kqb9t2ACXXgrnnAO33AI/+AGMGlV1VJIkSQ31\npdA7FFjT6fmr5bG+nNOXtgBfAx7oQyyStPs88wx86lOwciUsWwannVZ1RJIkSX0yUMsr9PmmlYi4\nEtiYmT/v6Zyrr776/f2WlhZaWlp2JTZJamzDBvjOd2DuXLjhBpg923vxJEnSbtPa2kpra+suvUZf\nCr3XgMM7PR9XHut6zmHdnDOyUduIOA+YBpzcKIDOhZ4kDaglS+CCC2DiRGhrg498pOqIJEnSMNP1\n4tY111zT9Gv0ZejmUuCoiDgiIkYCZwGLupyzCPgHgIj4NLAuMzsatY2IKcClwJmZ+V7TkUtSf1q3\nrijwvvxluPZauPdeizxJkjRo9VroZeZm4CKKWTKfBxZk5oqImBMRF5Tn/BJYGREvAf8B/GOjtuVL\n/zvwAeDhiHg6In7Sv29NkvogE+6+u7iCt9de8NxzMGNG1VFJkiTtksjMqmNoKCKy7jFKGqRWrSoW\nPX/5ZZg3Dz772aojkiRJ2kFEkJlNTRjQpwXTJWlIWb8evvWtYkbNE04oZte0yJMkSUOIhZ6k4WPT\nJrjtNpgwobiK19YGV14JI0dWHZkkSVK/GqjlFSSpPjJh4UK44gr48IfhrrvgM5+pOipJkqQBY6En\naWhbsgQuv7wYrvn978PUqa6JJ0mShjwLPUlDU1sbfPOb8Pzz8O1vF4uej3C0uiRJGh781CNpaHn6\n6WJ5hGnT4NRTob0dzjnHIk+SJA0rfvKRNDQ8/jiccQZ8/vNw8snFZCuXXAKjRlUdmSRJ0m7n0E1J\ng1cmPPggfPe7sGYNfP3rcM89sPfeVUcmSZJUKQs9SYPP+vVwxx1w443F0giXXQazZsGe/kqTJEkC\nCz1Jg8nq1fDjHxdr4X3uc3DTTXDSSc6iKUmS1IX36Emqt0z41a/gi1+E444rFj1/6im47z5oabHI\nkyRJ6oZX9CTVU0cHzJ9fXL3bsgUuvhh++lP44AerjkySJKn2LPQk1cfbb8MDD8CCBfDoo/CFL8C8\neTB5slfuJEmSmmChJ6lar70GixbBwoXw61/DiSfCzJnF1Tyv3kmSJO2UyMyqY2goIrLuMUpqQiYs\nX76tuHv55WJx8+nTYcoUiztJkqQuIoLMbGp4k4WepIHV0QHLlhWP3/wGHnusKOZOPx1mzChmz9xr\nr6qjlCRJqq0BK/QiYgrwbxSzdN6amdd1c84PganAeuC8zGxr1DYiDgDuBI4AVgGzMvOtbl7XQk8a\nDDZuhPb2bUXds88W2/feg2OPLR4nnFAUdocdVnW0kiRJg8bOFHq9Lq8QESOAHwGnAROBsyNiQpdz\npgJHZubHgDnA3D60vRx4JDP/AlgMfKOZwFUfra2tVYegBvq9f7Zsgd/9Dh55BK6/Hs49FyZNgv33\nhy99qRiSOWZMMUvm0qXwhz9Aa2uxuPns2RZ53TCH6s3+qT/7qP7so3qzf4amvkzGcjzwYmauBoiI\nBcB0oL3TOdOB+QCZ+WREjI6IscD4Bm2nAyeV7W8HWimKPw0yra2ttLS0VB2GetBU/2zeDG+8Aa+/\nXjzWrt22v2YNvPQSrFwJBx4IRx0FxxxTXKG78EL4+Mdh330H9L0MVeZQvdk/9Wcf1Z99VG/2z9DU\nl0LvUGBNp+evUhR/vZ1zaC9tx2ZmB0Bmro2Ig5qIWxqeMosFwzdtKoZKbty4/f6GDfDuu/CnP23b\nLl8Ot94Kb71VPNat23G7df+dd4oi7pBD4OCDi+0hh8CECXDKKUVxd+SRsN9+Vf9NSJIkqYGBWl5h\nZxa86vlGvDPOKM/ockqj582ca9tda7t2Ldx//+CKeTC03bx5x4Ju0ybYY49i8pI99yy2nff33ru4\nqrbPPtu2r7xSDKscMwZGj4bx44vtmDHbjm3djh5dvJYkSZIGtV4nY4mITwNXZ+aU8vnlQHaekCUi\n5gKPZuad5fN2imGZ43tqGxErgJbM7IiIg8v2R3fz5zsTiyRJkqRhrdnJWPry1f1S4KiIOAJ4HTgL\nOLvLOYuAC4E7y8JwXVnA/W+DtouA84DrgHOBhf3xhiRJkiRpuOu10MvMzRFxEfAQ25ZIWBERc4of\n5weh6IYAAAQ2SURBVM2Z+cuImBYRL1Esr/DVRm3Ll74OuCsivgasBmb1+7uTJEmSpGGo9gumS5Ik\nSZKa0+s6elWJiCkR0R4RL0TEZVXHox1FxKqIWBYRz0TEU1XHI4iIWyOiIyKe7XTsgIh4KCL+OyIe\njIjRVcY43PXQR1dFxKsR8XT5mFJljMNZRIyLiMUR8XxELI+IS8rj5lENdNM/F5fHzaGaiIhREfFk\n+dlgeURcVR43h2qiQR+ZRzUSESPKflhUPm86h2p5Ra9caP0F4BTgfyjuEzwrM9sbNtRuFRGvAJ/M\nzP+rOhYVIuJE4I/A/Mw8pjx2HfBmZn6v/NLkgMx0zcqK9NBHVwHvZOb1lQYnysnBDs7Mtoj4APBb\ninVfv4p5VLkG/fP3mEO1ERH7Zua7EbEH8DhwCTATc6g2euijqZhHtRER/wx8Etg/M8/cmc9zdb2i\n9/4i7Zm5Edi60LrqJajvv6FhKTMfA7oW3tOB28v924EZuzUobaeHPoKdW5ZG/Swz12ZmW7n/R2AF\nMA7zqBZ66J9Dyx+bQzWRme+Wu6Mo5oNIzKFa6aGPwDyqhYgYB0wDbul0uOkcquuH9J4WYFe9JPBw\nRCyNiPOrDkY9OigzO6D4kAQcVHE86t5FEdEWEbc4pKkeIuLPgUnAE8BY86heOvXPk+Uhc6gmyiFn\nzwBrgYczcynmUK300EdgHtXFDcClbL/OeNM5VNdCT4PD5Mw8juIbhwvLIWmqv/qN19ZPgI9m5iSK\n/3QdNlOxcljgPcA/lVeOuuaNeVShbvrHHKqRzNySmX9FcTX8+IiYiDlUK9300V9iHtVCRJwOdJSj\nFxpdYe01h+pa6L0GHN7p+bjymGokM18vt78H7qMYcqv66YiIsfD+/S1vVByPusjM3+e2G6bnAX9d\nZTzDXUTsSVFE3JGZW9d4NY9qorv+MYfqKTPfBlqBKZhDtdS5j8yj2pgMnFnOhfEL4OSIuANY22wO\n1bXQe3+R9ogYSbHQ+qKKY1InEbFv+Y0qEbEfcCrwXLVRqRRs/w3QIuC8cv9cYGHXBtrttuuj8hf2\nVn+HuVS124D/yswbOx0zj+pjh/4xh+ojIj60dchfROwD/C3FvZTmUE300Eft5lE9ZOYVmXl4Zn6U\nogZanJlfAf6TJnOolrNuQrG8AnAj2xZav7bikNRJRIynuIqXFDfx/sw+ql5E/BxoAf4M6ACuAu4H\n7gYOA1YDszJzXVUxDnc99NHfUNxrtAVYBczZOg5fu1dETAaWAMspfr8lcAXwFHAX5lGlGvTPbMyh\nWoiIT1BMFDGifNyZmf8aEQdiDtVCgz6aj3lUKxFxEvAv5aybTedQbQs9SZIkSdLOqevQTUmSJEnS\nTrLQkyRJkqQhxkJPkiRJkoYYCz1JkiRJGmIs9CRJkiRpiLHQkyRJkqQhxkJPkiRJkoYYCz1JkiRJ\nGmL+H5CTDfg1lJ76AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"b_x = np.where(V(x) > E_particle)[0][-1]\n",
"P_undecayed = [np.trapz(absl_psi[:b_x,n]**2,x[:b_x]) for n in np.arange(N_t)] \n",
"P_decayed = [np.trapz(absl_psi[b_x:,n]**2,x[b_x:]) for n in np.arange(N_t)] \n",
"\n",
"t = np.arange(N_t)/200\n",
"%matplotlib inline\n",
"fig1 = plt.figure(figsize = (15,6));\n",
"ax = fig1.add_subplot(211);\n",
"ax.plot(t,P_undecayed,'b');\n",
"ax = fig1.add_subplot(212);\n",
"ax.plot(t,P_decayed,'r');\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [default]",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}