{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Hydrogen Wavefunction

\n", "

\n", "After some variable changes the radial equation for the Hydrogen atom is:\n", "\n", "$$\n", "\\frac{1}{\\rho^2}\\frac{d}{d\\rho}\\Bigl (\\rho^2\\frac{dR}{d\\rho} \\Bigr)+\\Bigl [\\frac{\\lambda}{\\rho} -\\frac{1}{4}-\\frac{\\ell (\\ell +1)}{\\rho^2}\\Bigr ]R=0\n", "$$\n", "\n", " A form suitable for behavior at large $\\rho$ is:\n", "\n", "$$\n", "R(\\rho)=F(\\rho) e^{-\\rho/2},\n", "$$\n", "\n", "When substituted, we obtain the equation solved here using Runge-Kutta:\n", "$$\n", "\\frac{d^2F}{d\\rho^2}+\\Bigl ( \\frac{2}{\\rho}-1\\Bigr)\\frac{dF}{d\\rho}\n", " +\\Bigl [ \\frac{\\lambda-1}{\\rho}-\\frac{\\ell (\\ell +1)}{\\rho^2}\\Bigr ]F=0 \n", "$$ " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hc1bXw4d8a9S6ruKlYknvDtiw3TAu9JaYGE7ghCQlwAyk34SaQSnr5ciG5N4UQSEJJ6BBMqKZjcC+yJVe5yCq2JKtbtvr+/pgzMIhRs+fMmRmt93n0eHTarOPRzJqz19l7izEGpZRSqi+X0wEopZQKTpoglFJK+aQJQimllE+aIJRSSvmkCUIppZRPmiCUUkr5pAlCnTQRuUtEHnE6jkATkc+JyCqv34+KSMEQ9ssTESMikfZGODwicrqI7HI6DhU8NEGMQCJyQETO7bPsIx924cg6xx7rg7xFRIpF5FJ/Hd8Yk2iM2Xeyx7Fen+Mi0ioiTSLyvojcIiK2vl+NMe8aY6b2iePcgfZxmoj8XUQ6rdfU8xPhdFzhQhOEslWwfUsGVhtjEoFU4I/AYyKS6nBMvnzSGJMETAB+CXwbeMDZkILWr63k7PnpcTqgcKEJQn2MiPy3iDzdZ9n/ichvrcf5IvK29Q13JZDhtZ2n+eRGETkIvGEt/5SIlFrfiN8Skele+xSKyGbreE+KyOMi8lOv9ZeKyBavb9OneK07ICK3i8hWEWm29o0d7ByNMb3Aw0ACMNnreE+KyGHrWO+IyEyvdekissK6+lgHTOzzf2REZJL1+BLrnFpEpEJE7hospn7ibDbGrACuAW4QkVnW8WNE5DciclBEakTkXhGJs9adJSKVIvJNEakVkUMi8nmvOC8Wke3W/3eViNzuvZ/1+GEgF3je+lb+LRF5QUS+0uect4rIZX3j9vo7uMGK8YiIfPdE/g+Ug4wx+jPCfoADwLl9ln0OWGU9Hge0AanW75FALTDf+n01cDcQA5wBtAKPWOvyAAM8hPvDNw6YYh3vPCAK+BZQBkRbP+XA16x1VwCdwE+t4xVaz70IiABusOKP8TqXdcB4IA3YAdzSz3l7n2MEcKv1XKO9tvkCkGSd22+BLV7rHgOesM5rFlDlOZ613gCTrMdnAbNxfwk7BagBLuvzfxQ51NfHWn4Q+E/r8W+BFdY5JwHPA7/weu5u4MfW/+nFwDFglLX+EHC69XgUUOi1X2V/cQCfBtZ6/T4HqAeifcTqOce/WH8Dc4AOYHo/53wH0NTfzwB/y38HGqyfjcCVTr+/wunH8QD0x4EX3f3GP9rnTXisz4fdS8CXrMeXAtutx7nWh0+C17b/5OMJosBr/feBJ7x+d1kfrmfhTjBVgHitX8WHCeJPwE/6xL8LONPrXK73Wvdr4N5+zvtzVuxNQBdwHPj0AP9Pqda5pOBOKF3ANK/1P6efBOHjWL8F7unzfzTcBLEG+C4guBPuRK91S4D91uOzrHOL9FpfCyy2Hh8EbgaS+xz/LAZOEDHWB/Fk6/ffAH/s5xw855jttWwdsNzPf8uFQDruLzEX4/6ystTp91i4/GgT08h1mTEm1fMDfLnP+geB663H1+NujgH3N/VGY0yb17blPo5f4fV4vPc2xt28UwFkWeuqjPVu97HvBOCbVvNSk4g0ATnWfh6HvR4fAxJ9xOOxxjrfUbi/gZ/uWSEiESLySxHZKyItuD8gwd2Elon7Q8g7Nl/n7TnWIhF5U0TqRKQZuAWvprgTlIX7AzoTiAc2ev2fvGwt96g3xnR7/e79/3Il7g/TcqupcMlQntwY04H7Cup6q2B+LR/+XfRnOK/NsBljNhlj6o0x3caYF4F/4L4KVX6gCUL151/AKVab96W433jgbp4YJSIJXtvm+tjf+wO/GvcHPQAiIrg/5Kus42VZyzxyvB5XAD/zTmbGmHhjzKMnemIAxpijuJPif4jIPGvxZ4BlwLm4rxryPCEDdbivPrxj83XeHv/EnYByjDEpwL3WcU6IiCzAnSBWAUdwXyHM9Po/STHu4vugjDHrjTHLgNG4X+cn+tvUx7IHgeuAc4BjxpjVwzwVn0TkO33uRPrIzzAOZTiJ/2f1UZoglE/GmHbgKdwfdOuMMQet5eXABuBHIhItIqcBnxzkcE8Al4jIOSISBXwTd3v0+7jrGT3AbSISKSLLgIVe+/4FuMX6Ri4ikmAVgJP8cI71wP3AD6xFSVZc9bi/of/ca9se4BngLhGJF5EZuOsh/UkCGowx7SKyEHfyGTYRSRb3rbiP4W7G22Zdgf0FuEdERlvbZYnIBUM4XrSIXCciKcaYLqAF9/+/LzXAR/p1WAmhF/gfBr96GDJjzM/NR+9E+shPf/uJyFUikigiLhE5H/fV7gp/xTXSaYJQA3kQd6G17wfBZ3AXjRuAH+IuSPfLGLML9xv3/3B/+/0k7ts4O40xnbibBG7EXRu4Hvg37g9qjDEbgC8BvwcacRe3P3fyp/aB3wIXW3dGPYS72agK2I67zd/bbbibSA7jLo7+bYDjfhn4sYi04k5A/X1L78/z1r4VuOsOdwOf91r/bdz/F2us5rDXgKkfO4pv/wEcsPa7hQ+bEvv6BfA9qxnrdq/lD+H+uwiGzpFfw/16NQH/D3fd7C1HIwoj8tGmX6U+JCK5wE5grDGmJYDPuxZ3oXmgD2DlEBH5LHCTMeY0p2NR9tIrCOWTVYT8BvCY3clBRM4UkbFWE9MNuG8LfdnO51QnRkTicV8d3ed0LMp+wdbLVQUBqwBdg7u55cIAPOVU3E0wicBe4CpjzKEAPK8aBqvG8Qzu5qx/OhyOCgBtYlJKKeWTNjEppZTyKWyamDIyMkxeXp7TYSilVEjZuHHjEWNMpq91YZMg8vLy2LBhg9NhKKVUSBGRfkcE0CYmpZRSPmmCUEop5ZMmCKWUUj5pglBKKeWTJgillFI+2ZogRORCEdklImUicoeP9THiniKyTETWikietTxP3JO2b7F+7rUzTqWUUh9n222uIhIB/AH3NJOVwHoRWWGM2e612Y24J5+ZJCLLgV/hnnsXYK8xZq5d8SmllBqYnf0gFgJlxph9ACLyGO7JWLwTxDLgLuvxU8Dv+0wco5QaRGd3L2v21bPrcCudPb2MT43l1IkZjEmOdTo0FeLsTBBZfHR6xkrccwj43MYY021NzZhurcsXkc24JzT5njHm3b5PICI3ATcB5OYONLmXUuGnq6eXh1aX84c3y2ho6/zIOhE4f8YYvnXhNCZm+nWWTzWC2JkgfF0J9B0ZsL9tDgG5xph6EZkP/EtEZvYddtoYcx/WsMNFRUU66qAaMWpb2rn1n5tYf6CR0yZl8LlT81iQl0ZMlIu9dUd5cdshHlpdzkW/e5dvXTCVG0/LRy/O1XDZmSAq+ej8vdm45yb2tU2liETinge4wZrA3jOj2EYR2QtMwT3VpVIj2uHmdpbft5qalg5+t3wuy+ZmfWT9zPEpzByfwg2n5vGdZ0r46Qs72Ft3lJ8sm0VkhN64qIbOzr+W9cBkEckXkWhgOR+fK3YFH87rexXwhjHGiEimVeRGRAqAycA+G2NVKiS0tndx3f1rOHK0k0e+uPBjycHb6KRY7vuP+dz6iYk8uq6Cbz+9jd5evdBWQ2fbFYRVU7gNeAWIAP5qjCkVkR8DG4wxK4AHgIdFpAz3/MbLrd3PwD2fbzfuCdVvMcY02BWrUqGgt9fwjSeKOVB/jEduXMT8CWmD7uNyCf99wTRiIiO4e+VuRsVH8b1LZwQgWhUObB3N1RjzIvBin2U/8HrcDlztY7+ngaftjE2pUPPg6gOs3F7DDy6dwZKJ6YNu7+0rZ0+ioa2T+1ftZ8b4ZK4ozLYnSBVWtEFSqRBQ0XCMX7+8izOnZPL5pXnD3l9E+O4l01lckMadz2xjxyFbpxlXYUIThFIh4PvPleAS+Nnls074bqSoCBe//0whSbFR/NfjW+jo7vFzlCrcaIJQKsit2nOEt3bV8bVzJ5M9Kv6kjpWRGMMvr5jNzsOt/O61PX6KUIUrTRBKBbHeXsMvXtpBVmocn12S55djnjtjDNcU5XDv23spqWr2yzFVeNIEoVQQe7HkEKXVLdx+wRRioyL8dtzvXDKdUfHR/OC5Er31VfVLE4RSQcoYw5/e2ktBZgLL5vTf3+FEpMRFcefF09l0sImnNlX69dgqfGiCUCpIrSo7Qml1CzefUYDL5f9hMq6Yl0XRhFH88qWdtLZ3+f34KvRpglAqSN379l7GJMdw2Tz/Xj14uFzCDz85k4a2Tv7yjg5UoD5OE4RSQWjn4RbeK6vn80vziYn0X+2hr9nZKVxyyjj+8u5+alvbbXseFZo0QSgVhB5de5DoCBfXFOUMvvFJuv38qXT19PJ/r5fZ/lwqtGiCUCrIHO/s4ZnNVVw0eyyjEqJtf778jASuWZDDo+sOcrD+mO3Pp0KHJgilgsy/t1bT2t7NtQsDNwnWV8+ZjEuEP729N2DPqYKfJgilgsyj6w5SkJnAovzBR2v1lzHJsXx6QTZPbazgUPPxgD2vCm6aIJQKIgeOtLHpYBOfLsoJ+AxwN58xEWPgPr2jSVk0QSgVRJ4vdk+6+Kk54wP+3Dlp8Vw2L4tH1x3kyNGOgD+/Cj6aIJQKEsYYniuuZmF+GuNT4xyJ4ctnTaSju5cHVu135PlVcNEEoVSQ2HGolbLao45cPXgUZCZy0ayx/GNNOW0d3Y7FoYKDJgilgsRzxVVEuoSLZ49zNI4bTyugpb2bp3WMphFPE4RSQcAYw7+LD3H65AzSAtD3YSDzJ4xibk4qf3vvgI70OsJpglAqCGw/1EJV03EumuXs1YPHjafls/9IG2/srHU6FOUgTRBKBYGV22sQgbOnj3Y6FAAumjWW8Smx3L9Kb3kdyTRBKBUEVm6vYX7uKDISY5wOBYDICBefW5rHmn0NOuvcCKYJQimHVTUdp7S6hfNmjHE6lI+4ZkEu8dERPLy63OlQlEM0QSjlsJWlhwGCLkGkxEWxbO54niuuovm4Tig0EmmCUMphK3fUMGl0IgWZiU6H8jHXLZpAe1cvz+gtryOSJgilHHS0o5u1+xo4J0iK033NykphTk4q/1h7EGP0lteRRhOEUg5avbee7l7DmVMynQ6lX9cvyqWs9ihr9jU4HYoKME0QSjno3T11xEdHUDQhcEN7D9cn54wnJS6KR9ZqsXqk0QShlIPe2V3HkoJ0oiOD960YGxXBVfOzeaXksM5bPcIE71+lUmHuYP0xDtQf4/TJGU6HMqjrFuXS3Wt4Yn2F06GoALI1QYjIhSKyS0TKROQOH+tjRORxa/1aEcnrsz5XRI6KyO12xqmUE97eUwfAGUFcf/AoyEzk1InpPLa+QsdnGkFsSxAiEgH8AbgImAFcKyIz+mx2I9BojJkE3AP8qs/6e4CX7IpRKSe9u7uO7FFx5GckOB3KkFyzIIfKxuOs3lfvdCgqQOy8glgIlBlj9hljOoHHgGV9tlkGPGg9fgo4R6x5FkXkMmAfUGpjjEo5oqunl/f31nP65MyATy16oi6YOZbk2Eie2KDNTCOFnQkiC/D+S6q0lvncxhjTDTQD6SKSAHwb+NFATyAiN4nIBhHZUFdX57fAlbLbloomjnZ0c0YI1B88YqMiWDY3i5dKDtN8THtWjwR2JghfX4v6Nl72t82PgHuMMUcHegJjzH3GmCJjTFFmZvC34yrlsWZvPSKwZGK606EMyzULcujs7mVFcZXToagAsDNBVAI5Xr9nA9X9bSMikUAK0AAsAn4tIgeArwPfEZHbbIxVqYBas7+eaWOTSY13dnKg4ZqVlcKMcck8rs1MI4KdCWI9MFlE8kUkGlgOrOizzQrgBuvxVcAbxu10Y0yeMSYP+C3wc2PM722MVamA6ezuZWN5I4sLgrdz3EA+XZRNSVULpdU6DHi4sy1BWDWF24BXgB3AE8aYUhH5sYh8ytrsAdw1hzLgG8DHboVVKtxsrWyivauXRfmh1bzkcdm8LKIjXDy5QQfwC3eRdh7cGPMi8GKfZT/wetwOXD3IMe6yJTilHLLGuk10UX5oXkGkxkdz/swxPLu5ijsumkZsVITTISmbaE9qpQJszb4Gpo1NYlRCaNUfvF2zIIfm412s3F7jdCjKRpoglAqgzu5eNpQ3sLggNJuXPJZOzCArNU77RIQ5TRBKBdC2Knf9IVQL1B4ul3BlYRbvlR3hcLMO4BeuNEEoFUCeORUWhmiB2tvlhdn0GvjXFu0TEa40QSgVQGv3NzBlTCJpIVx/8MjPSKAwN5WnN1bqbHNhShOEUgHS02vYXN5IUV5oNy95u3J+Nntqj1JS1eJ0KMoGmiCUCpA9ta20dnQzP3eU06H4zaWzxxMd6eLpTdonIhxpglAqQDaWNwIwf0L4JIiU+CjOmz6GFcXVdHb3Oh2O8jNNEEoFyMbyRjISo5mQHu90KH51RWEWDW2dvL1bR1QON5oglAqQTeWNFOaOCpn5H4bqjCmZpCdE84w2M4UdTRBKBcCRox0cqD8WVs1LHlERLpbNzeL1HbU0Het0OhzlR5oglAqATWFYf/B2RWEWnT29PL/1kNOhKD/SBKFUAGwsbyQqQpiVleJ0KLaYOT6ZaWOTtJkpzGiCUCoANpY3MisrJWxHPhURrijMYvPBJvbWDTgRpAohmiCUsllHdw9bq5rDqv+DL5fNzcIl8OwmHXojXGiCUMpmpdUtdHb3hm39wWN0ciynT87k2c1V9Pbq0BvhQBOEUjbzFKgLwzxBgLtYXdV0nDX7650ORfmBJgilbLalooms1DjGJMc6HYrtLpg5lqSYSJ7eqM1M4UAThFI2K65s4pTs8Lx7qa/YqAgunj2Ol0sOcayz2+lw1EnSBKGUjRraOqloOM6cnFSnQwmYK+dn09bZwyulh50ORZ0kTRBK2ai4sglgxFxBABRNGEVOWpw2M4UBTRBK2WhrRTMiMDtMO8j54nIJV8zL5r29RzjUfNzpcNRJ0AShlI2KK5uYlJlIUmyU06EE1BWFWRgDz27Wq4hQpglCKZsYY9ha2cQp2SOn/uAxIT2BBXmjeGZTlU5HGsI0QShlk6qm4xw52sncnJHTvOTtisJsymqPsrWy2elQ1AnSBKGUTYor3B+MI+kOJm+XnDKO6EiXDuAXwjRBKGWTrZVNREe4mDY22elQHJEcG8X5M3Q60lCmCUIpm2ypaGL6+GSiI0fu2+zK+dk0HuvizV21ToeiTsDI/ctVykY9vYaSqmbmjqD+D76cPimDzKQYnt6ozUyhSBOEUjbYW3eUts6eEXkHk7fICBeXzR3Pm7tqaWjT6UhDjSYIpWxQXOHuQT1SC9TerijMpqvH8HxxtdOhqGGyNUGIyIUisktEykTkDh/rY0TkcWv9WhHJs5YvFJEt1k+xiFxuZ5xK+VtxZRNJMZEUZCQ4HYrjpo9LZsa4ZJ7Wu5lCjm0JQkQigD8AFwEzgGtFZEafzW4EGo0xk4B7gF9Zy0uAImPMXOBC4M8iEmlXrEr527bKZmZlpeByidOhBIUrCrPYWtnMnppWp0NRw2DnFcRCoMwYs88Y0wk8Bizrs80y4EHr8VPAOSIixphjxhjPWMGxgHbFVCGjq6eXHYdbmT3CC9Tels3NIsIlPK3TkYYUOxNEFlDh9XultcznNlZCaAbSAURkkYiUAtuAW7wSxgdE5CYR2SAiG+rq6mw4BaWGr6z2KJ3dvcwcPzL7P/iSmRTDmVMy+dfmKnp0OtKQYWeC8HVt3fcvo99tjDFrjTEzgQXAnSLysem4jDH3GWOKjDFFmZmZJx2wUv5QUuXuQT1rBI3gOhRXFmZzuKWd9/cecToUNUR2JohKIMfr92yg720MH2xj1RhSgAbvDYwxO4A2YJZtkSrlR6XVLcRHR5CfrgVqb+dMH01ybCTPaDNTyLAzQawHJotIvohEA8uBFX22WQHcYD2+CnjDGGOsfSIBRGQCMBU4YGOsSvlNaXUzM8Yla4G6j9ioCC6dM56XSw5ztEOnIw0FQ0oQIvK0iFwiIkNOKFbN4DbgFWAH8IQxplREfiwin7I2ewBIF5Ey4BuA51bY04BiEdkCPAt82Rij16Uq6PX2GrZXt2jzUj+uLMzmeFcPL2075HQoagiGeuvon4DPA/8rIk8CfzfG7BxsJ2PMi8CLfZb9wOtxO3C1j/0eBh4eYmxKBY0D9W20dfYwQwvUPhXmppKfkcDTmyq5uihn8B2Uo4Z0RWCMec0Ycx1QiLupZ6WIvC8inxeRkTVVllIDKKluAWDWeL2C8EVEuGJeFmv2NVDZeMzpcNQghtxkJCLpwOeALwKbgd/hThgrbYlMqRBUWtVMdISLyWMSnQ4laF02z323+7NarA56Q61BPAO8C8QDnzTGfMoY87gx5iuAvhOUspRWtzB1bBJRETrMWX9y0uJZXJDGM5t1OtJgN9S/4vuNMTOMMb8wxhwC9zhKAMaYItuiUyqEGGMoqW5mVpbWHwZzRWE2+4+0selgo9OhqAEMNUH81Mey1f4MRKlQV9V0nKZjXczU+sOgLp49jvjoCJ5YrwP4BbMBE4SIjBWR+UCciMwTkULr5yzczU1KKUupVaDWITYGlxgTyaWnjOP5rdXaJyKIDXab6wW4C9PZwN1ey1uB79gUk1IhqbSqmQiXMH2cJoihuGZBLk9sqOSFrdVcsyDX6XCUDwMmCGPMg8CDInKlMebpAMWkVEgqqW5hUmYisVERTocSEgpzU5k0OpHH11dogghSAyYIEbneGPMIkCci3+i73hhzt4/dlBqRSqqaOW1ShtNhhAwRYfmCHH76wg5217QyZUyS0yGpPgYrUntGG0sEknz8KKWA2tZ2als7mKlDbAzL5fOyiIoQHl9fMfjGKuAGa2L6s/XvjwITjlKhqfSDHtRafxiO9MQYzpsxhmc3V/GtC6cSE6nNc8FkqB3lfi0iySISJSKvi8gREbne7uCUChWl1hwQOgbT8F2zIJeGtk5e217rdCiqj6H2gzjfGNMCXIp7DocpwH/bFpVSIaa0uoW89HiSYnVosuE6bVIG41NieXyDNjMFm6EmCM9f/cXAo8aYhoE2VmqkKalu1vrDCYpwCVcX5fDunjodwC/IDDVBPC8iO4Ei4HURyQTa7QtLqdDRfKyLiobj2kHuJFxdlA3Akxu0Z3UwGepw33cAS4AiY0wX7ilAl9kZmFKhovSQNQe1DrFxwrJHxXPapAye2lhJT68O4BcshjPk5HTgGhH5LO7pQc+3JySlQktplQ6x4Q/LF+RS1XScd/bUOR2KsgxpRjkReRiYCGwBeqzFBnjIpriUChkl1c2MS4klPTHG6VBC2nkzxpCRGMM/1pTziamjnQ5HMfQpR4uAGUYHb1fqY0qrW3QEVz+IjnSxfEEOf3yrjMrGY2SP0vFAnTbUJqYSYKydgajQV9FwjJKqZjq7e50OJWCOdXazt+6ozgHhJ9cuco/J9Oi6gw5HomDoVxAZwHYRWQd0eBYaYz5lS1QqpNS2tnP7k1t5Z7e77TgzKYafXTaL82eG/3eKHYdaMAa9gvCTrNQ4zp42hsfXV/C1c6YQHakz8zlpqAniLjuDUKGrsa2T5X9eQ3Xzcb514VSyR8Vz3zt7ufmRjfz2mrksm5vldIi2+mCIDb2C8JvrF+fy2o4aXik9zCfnjHc6nBFtSAnCGPO2iEwAJhtjXhOReEAHTRnhjDF86+mtVDYd55EbF7EwPw2A86aP4Ya/reNbT21l2thkpo4N33EdS6qaSU+IZmxyrNOhhI0zJmeSmxbPw2vKNUE4bKhjMX0JeAr4s7UoC/iXXUGp0PByyWFWbq/h9vOnfJAcAOKiI/jjdYUkxERy5zNb6Q3j+9pLqlqYMT4ZEXE6lLDhcgnXLcpl3f4Gdte0Oh3OiDbUBr5bgaVAC4AxZg+g96GNYD29hv/36i6mjkniC0vzP7Y+IzGG71w8nU0Hm3h+a7UDEdqvo7uH3TWtzNIhNvzu6qIcoiNd/GNNudOhjGhDTRAdxphOzy8iEom7H4Qaof69tZp9dW187dzJREb4/jO6Yl4WU8ck8bvX94Rl79g9NUfp7jXag9oGaQnRXDJ7HE9vqqJN56x2zFATxNsi8h0gTkTOA54EnrcvLBXsHli1n0mjE7lwgDuVXC7ha+dOZl9dGy+VHApgdIFRYg3xrT2o7XH94lyOdnTzry1VTocyYg01QdwB1AHbgJuBF4Hv2RWUCm4lVc1srWzm+kW5uFwDt71fOHMsuWnxPPR++DUVlFQ3kxQTSW6aduiyQ2HuKGZlJfO39w6gfXSdMdTB+npxF6W/bIy5yhjzF+1VPXL9c91BYiJdXF6YPei2Lpdw/eJc1h1oYOfhlgBEFzil1e4C9WBJUp0YEeELS/Mpqz3Ku3uOOB3OiDRgghC3u0TkCLAT2CUidSLyg8CEp4JNR3cPzxdXc8nscaTEDW1ynKvn5xAT6eKfa8Ond2x3Ty87DukQG3a75JRxZCbF8Nf39jsdyog02BXE13HfvbTAGJNujEkDFgFLReS/Bju4iFwoIrtEpExE7vCxPkZEHrfWrxWRPGv5eSKyUUS2Wf+ePewzU7Z4r+wIre3dw7o/fVRCNOfOGMO/tx6iqyc8huHYd6SN9q5erT/YLCYygv9YPIG3dtVRVnvU6XBGnMESxGeBa40xH6RvY8w+4HprXb9EJAL4A3ARMAO4VkRm9NnsRqDRGDMJuAf4lbX8CPBJY8xs4Abg4aGdjrLbi9sOkxQbydJJGcPab9mc8TS0dbIqTJoKSqvdBerZ2XoFYbfPLMolOtLF39/Xq4hAGyxBRBljPvaONsbU8eE0pP1ZCJQZY/ZZt8g+xscnGVoGPGg9fgo4R0TEGLPZGOO5eb4UiBURHUvZYZ3dvbxaepjzZowZ9hg5Z00dTUpcFM+FyR0pJVUtxEa5KMhIcDqUsJeRGMNlc8fz9MYqmo51Dr6D8pvB3uUDvRqDvVJZgPcs5JXWMp/bGGO6gWYgvc82VwKbjTEdKEet2VdPS3s3F88aN+x9oyNdXDx7LK+U1nC8s2fwHYJcSVUz08Ym99sHRPnX55fmc7yrh8fWVwy+sfKbwf6654hIi4+fVmD2IPUgozAAACAASURBVPv6urWj751PA24jIjNxNzvd7PMJRG4SkQ0isqGuTmehsttbu+qIiXRx2uThNS95XDx7HMe7enivLLSbmXp7DdurW3SAvgCaPi6ZUyem8+D7B8KmjhUKBkwQxpgIY0yyj58kY8xgTUyVQI7X79lA3zEXPtjG6p2dAjRYv2cDzwKfNcbs7Se++4wxRcaYoszMzEHCUSfr7d21LCpIJzbqxMZpXJSfTlJMJCu31/g5ssCqaDxGa0e39qAOsC8szedQczsvlRx2OpQRw87r4/XAZBHJF5FoYDmwos82K3AXocE9z/UbxhgjIqnAC8Cdxpj3bIxRDVFFwzH21rVx5pQTT8TRkS7OnJrJ6ztrQnrojZIqzxDfmiAC6expo5mYmcC9b+3VjnMBYluCsGoKtwGvADuAJ4wxpSLyYxHxTDT0AJAuImXAN3D32MbabxLwfRHZYv3o4IAO8kwkf+aUE2te8jhvxhiOHO1kS0WTP8JyREl1M5EuYfKYRKdDGVFcLuHmMyey/VAL74TJ3XDBbqgTBp0QY8yLuIfl8F72A6/H7cDVPvb7KfBTO2NTw/PO7jqyUuOYmHlyH4pnTR1NpEtYub2G+RNG+Sm6wCqtbmHKmCRiInVKlEC7bG4Wd7+6m3vf2ntSV7NqaPQWDDWonl7D+3vrOX1yxknPe5ASF8XC/DTe3Fnrp+gCyxhDaVWzFqgdEh3p4oun57N6Xz2bDzY6HU7Y0wShBrXjUAut7d0sLuh7B/KJOWNKJrtqWqltaffL8QLpcEs79W2dWn9w0PKFuaTERXHv2z7vXVF+pAlCDWrd/gaAj8wadzJOt26TDcUB2DwFah1iwzmJMZHcsGQCr26v0eE3bKYJQg1q3f4GctLiGJ8a55fjTR+bTHpCNO/uCb2+K6XVzYi478tXzrnh1DxiIl3c945eRdhJE4QakDGGdQcaWJjnn+YlcN+NctrkDFaVHQm5+apLqlqYmJlIfLSt93eoQaQnxrB8QS7Pbq6isvGY0+GELU0QakBltUdpaOtkkZ+alzxOn5zJkaOd7AixOSJKq5uZpc1LQeGmMwoQhD+8qVcRdtEEoQa01qo/LCrwd4IIvTpE/dEODjW36xwQQWJ8ahzXLMjhyQ0VehVhE00QakDr9jcwJjnG79NqjkmOZcqYxJCqQ5RWWwVqvcU1aHz5ExNxifCHN8ucDiUsaYJQA9pY3siCvLST7v/gy9JJGWwsb6SjOzRGdy2x5oDQK4jgMS4ljuULc3hyQyUVDXoV4W+aIFS/alvbqWo6ztycVFuOv7ggnfauXrZWNttyfH8rrWohNy1+yFOtqsD4z7P0KsIumiBUv7ZWuD+47UoQC/PcdY01e+ttOb6/lVY3a/+HIDQuJY5rF+bw1Ea9ivA3TRCqX8WVTUS4xLYmlVEJ0Uwbm8Sa/cGfIFrauzhQf0x7UAep/zxrEi6X8NvX9jgdSljRBOGQju4entpYyZ/f3kt5fZvT4fhUXNnMlDFJxEXbNyjd4oL0kKhDbK/WHtTBbGxKLDcsmcAzmyvZGWK3TgczTRAOaGnv4tP3rub2J4v5xUs7Of+ed3h9R3BNomOMobiiibk59n5jDpU6xDYrPr2CCF5fPmsSiTGR/PrlXU6HEjY0QQSYMYbbnyimtLqFP3ymkPfvOJspY5L4yqOb2VcXPOPKlNcfo/l4F3Oy7ak/eHg64K3dF9zNTFurmslKjSMjMcbpUFQ/RiVE8+WzJvHGzlrWBPnfU6jQBBFgb+2q49XtNdx+wVQuOWUc41Pj+Mtni4gQ4Sf/3u50eB8ornRP6HOKzQnigzrEvgZbn+dkba1s4pRsvXoIdp9fmsfY5Fh++dJOnXXODzRBBJAxhl+9vJO89Hi+sDT/g+VjU2L56jmTeXNXHesPBMcHZXFFM7FRLqYEYNa0xQXpbChvoLM7OCejbzrWSXn9MWZrggh6sVERfOO8KWypaOJlnbv6pGmCCKDVe+vZebiV286eTHTkR//rr188gdT4KP7yzj6Hovuo4somZmelEBlh/5/Ih3WI4JyGdFuVu/5gd3Ob8o8rCrOYPDqRX7+yK2i/dIQKTRAB9NDqckbFR3HpKeM+ti4uOoLrF01g5Y4ax+/l7urppaSq2fbmJQ9PHWJ1kPaH2KoF6pASGeHiOxdPZ/+RNh58/4DT4YQ0TRABUn+0g5U7avh0UQ6xUb5vG71mQQ7GwHNbqgIc3Uftrmmlo7uXOTZ1kOvLU4fwDAwYbIormsjPSNAe1CHkE9NGc/a00fzu9T3UtobezIXBQhNEgLxcepieXsOyuVn9bpOTFs/CvDSe3VzlaIGt2NODOoBNKovy09hY3khXT/A1CWyratYCdQj6/qUz6Oju0dteT4ImiAB5cdsh8jMSmD4uacDtLi/MYm9d2wcjhzqhuKKJUfFR5KT5Zwa5oVhUkM7xrp6g6w9R29rOoeb2gDW3Kf/Jz0jgC6fl89TGSrZUBGd9K9hpggiA+qMdrN5bzyWzxw06Kur5M8bgEli53bmOc8WVTZySnWrLCK798cx3vTbIht3wjEelVxCh6StnTyYzKYYfrigNudkLg4EmiAB4fWctvQYumj120G3TE2MozB3F6zudSRDHOrvZXdMasPqDR0ZiDJNGJ7I2yPpDbK1qxiU6xEaoSoyJ5M6LplFc0cRj6yucDifkaIIIgDV760lPiGbGECe6P2f6GEqqWjjUfNzmyD6upKqFXgNzHPjGvCg/jQ0HGugOojrE1sompoxJ0jmoQ9jl87JYXJDGL17aQW2LFqyHQxOEzYwxrNlXz+KC9CE32Zw3YzQAr+2otTM0n4orAtOD2pdFBem0dfY4Wn/xZoxha6UWqEOdiPDzy2fT0d3Lj4JotIJQoAnCZhUNx6lubmfxMOZ0npiZSE5aHO/sDvx0nMWVTWSlxpGZFPgxhxYHWR2isvE4DW2dzNYCdcgryEzkK5+YxAtbDwXdwJjBTBOEzTyDhi0uSB/yPiLC0okZrNlXH/DmluLKJubYPIJrf0Ynx5KfkRA0dQjPeFRONLcp/7v5zIlMHp3I9/9VQltHt9PhhARNEDZbs89df5g0enhjGp06KYPW9m5KAtjcUn+0g4qG444OKbEoP411BxroCYI7TjaVNxEb5WL6EGtHKrhFR7r45ZWzqW5u59cv73Q6nJCgCcJGJ1J/8Dh1ovuK472yI3aE5tNWz5hDAb6DyduigjRa27vZccj5OsSmg42ckpVKVADGo1KBMX9CGp87NY8HV5fz7p7AN+GGGv3Lt9GJ1B88MhJjmDY2KaAJoriiCZfAbAfHHFqU706MTg+70d7VQ2l1M/MmaP0h3Nxx0TQmZibw309upflYl9PhBDVbE4SIXCgiu0SkTETu8LE+RkQet9avFZE8a3m6iLwpIkdF5Pd2xminE6k/eFs6KYMN5Y20dwVmOs7iiiYmjU4kIca5WzrHp8aRkxbn+ARCpdXNdPUYCnNHORqH8r/YqAju/vRc6o52cNfzpU6HE9RsSxAiEgH8AbgImAFcKyIz+mx2I9BojJkE3AP8ylreDnwfuN2u+ALhROsPHksnpdPZ3cum8kY/R/ZxxhiKK5uDYkjrRfnprDvQ4GjP103l7gK1JojwNCcnla+cPYlnN1fx4rZDTocTtOy8glgIlBlj9hljOoHHgGV9tlkGPGg9fgo4R0TEGNNmjFmFO1GEpJOpP3jMn5CGCKwLwCRCnls6naw/eCzKT6PpWBe7a1sdi2HTwUZy0py53VcFxq2fmMSc7BTufGYblY3ODrEfrOxMEFmAd9/2SmuZz22MMd1AMzDk9hgRuUlENojIhrq64Co4nUz9wSMlLoppY5MDMsvch7d0Op8gPE1yTt3uaoxh08FGvXoIc1ERLv732nn09hpu++dmnVzIBzsThK+vzX3bDIayTb+MMfcZY4qMMUWZmZnDCs5uJ1t/8FiYN4pN5U22D4NdXNFEdKSLqWMHHm02ELJHxZGVGudYh7nq5nZqWjo0QYwAE9IT+NVVp7Cloolf6a2vH2NngqgEcrx+zwaq+9tGRCKBFCA4ekmdpJOtP3gsyE/jeJf9w08UVzYzY1zyx6ZCdYKIuPtD7G9wZF4MT81HE8TIcPHscdywZAIPrNrPq6U6j7U3Oz8N1gOTRSRfRKKB5cCKPtusAG6wHl8FvGGcnCnHT/xRf/BYmOduotpgYzNTd08v2yqbmRsE9QePRQVpHDnayd66owF/7k0HG4mNcjFtkLk7VPj4ziXTOSU7hW8+Wcw+B/7mgpVtCcKqKdwGvALsAJ4wxpSKyI9F5FPWZg8A6SJSBnwD+OBWWBE5ANwNfE5EKn3cARW0/FF/8BidHMuE9HjW2dgvYE/tUY539QRXgrD6Q6x2oA6xsbyRU7K1g9xIEhMZwR8+U0hUhIsvPrSB5uPaPwJs7gdhjHnRGDPFGDPRGPMza9kPjDErrMftxpirjTGTjDELjTH7vPbNM8akGWMSjTHZxpiQGYbRX/UHj6IJaWwob7StucUz21Yw3MHkMSE9njHJMQHvD9Ha3kVJVfMHAweqkSMnLZ4/XVfIwfpjfPXRzUEx3IvT9CuSDdbsqycj8eTrDx4L80fR0GZfc0txRRMpcVHkpcfbcvwT4a5DpLM2wHWIjeWN9BpYmO+f5K5Cy6KCdH5y2Sze3l3HL17c4XQ4jtME4Wee+sMiP9QfPBZYdYh1++3pMLelook5OYGdYnQoFhWkUdfawf4jbQF7zrX7G4h0CYU6xMaIde3CXG5YMoH7V+3nodUHnA7HUZog/OzD+oP/voHmZySQkRhtS6G6rcM9xWgw1R88nBiXad3+Bk7JTtEZ5Ea47186g3Onj+aHK0pHdE9rTRB+5qk/LPFDgdpDRFiQl2ZLj+ptVc30Gpjr0BwQA5mYmUBGYuDqEMc7e9ha2aTNS4rICBf/d20hhbmj+PpjW1i9NzgmsQo0TRB+5qk/TMz0T/3BY2F+GpWNx6lu8u881Z4pRoOhB3Vfnv4QgapDbDrYSFePYZEfk7sKXXHRETxwQxET0uO56aENlFjD4Y8kmiD8yI76g8fCfE8dwr9XEVsqmshJiyM9MTjHHFpUkMah5nYONtg/Vs7a/Q24BIomaAc55ZYaH82DX1hIclwU192/dsQlCU0QfmRH/cFj2thkkmIj/d4eX1zRxNyc4P1AXDopA4B399g/L8aaffXMGJ9MUmyU7c+lQsf41Dgeu2kxiTGRXP/AWkqrR06S0AThR3bUHzwiXFYdwo/jE9W2tFPd3B7Ucy4XZCSQlRrH27vtHYyxraObzQcbP0hISnnLSYvn0S8tJj4qYkRdSWiC8KPVNtUfPBbmp7G3ro0jRzv8crxNB931h3m5wVd/8BARzpyayftlR2wdbXPt/nq6egxnTA6uQR9V8MhNj+exm5YQHxXB8vvWjIjCtSYIPzHG8P7eI7bUHzw8dYj1fmpm2nCggehIF7McnGJ0KM6ckklbZw+bDto3cdI7u48QE+livtYf1ABy0+N56j9PZVxKLDf8dR0vhfktsJog/GTfkTZqWjpYOtG+JopZ41OIi4rwWx1ifXkjc3NSiYmM8Mvx7HLqxHQiXWJrM9OqMndyj40K7v8L5bzxqXE8ecsSZmUl8+V/buLh1QecDsk2miD85H3rcvPUifbdQx8d6aJwQqpfEsSxzm5Kq5pZkBf835iTYqMonDCKt3fZkyCqm45TVnuUMyZr/UENTWp8NP/44mLOnjqa7z9Xyvf+tc32OVucoAnCT1bvPcL4FPfIq3ZalJ/OzsMtNB87udEmtxxsorvXUJQXGvf8nzklk+2HWqht9f8stKusO6RO0wShhiEuOoL7PlvEzWcU8Miag1x//1rq/VQfDBaaIPygt9ewem89p07KsH08o4X5aRgDG8pP7ipi/YFGREJnUpwzp7iLx3ZcRby1u5bRSTFMHaPzP6jhiXAJd148nXuumcPmiiY+9fv3bK2VBZomCD/YcbiFxmNdtjYveczNSSU6wnXSHeY2lDcwbWwyKXGhcc//zPHJjEuJ5dXtNX49bkd3D2/vquOc6WOCbrBCFToun5fNkzcvQQQ+fe9q/vTWXnrDYLhwTRB+8H6Zp/5gfxNFbFQEc3JSTqoO0d3Ty6byxpCoP3iICBfMHMs7u+s41tntt+O+v7eets4ezp8xxm/HVCPTnJxUXvjq6Vwwcyy/enknn/3rOg43+79JNJA0QfjBW7trmTw6kbEpsQF5voX5aWyraqal/cTqEFurmmnr7PngttlQcf7MMXR09/q1mem17TXER0ewJABXfyr8pcRF8fvPzOMXV8xmQ3kD593zNo+tO+jI3Or+oAniJLW2d7F2XwPnTA/cN9DTJ2fS02s+uHIZrlV7jiCCrbfk2mFhXhqj4qN4xU8Ty/f2Gl7bUcOZUzL19lblNyLCtQtzeflrZzBjXDJ3PLON6x9Yy8F6+8cT8zdNECfp3T1H6O41nDN9dMCeszB3FIkxkbyz58S+Sa/ac4RZ41MYlRDt58jsFRnh4pzpY3h9Z61felVvq2qmpqWD87R5SdkgLyOBR7+0mJ9eNostB5s49563+Z9Xd/m1idRumiBO0us7akmNj2JeACfciY50sWRiOm/vqhv2pWtbRzebQnjMoQtnjqW1vZtVZSffzLSiuJroCBfnTNMEoezhcgnXL57Aa988k4tmjeX/3ijj7N+8zbObK0OiiK0J4iR09/Ty1q5azpqSSWREYP8rz5ySSVXTcfYNczrOtfvr6e41nB6i9/yfMSWTUfFRPLOp6qSO09NreL64mrOmZpISHxp3cqnQNS4ljt8tn8dTtywhMymG/3q8mIv/911eKT0c1PUJTRAnYfW+eurbOrlw1tiAP/eJ9gt4e1ddSI85FB3p4lNzxvPq9hqaj594Z8HVe+upbe3gsnlZfoxOqYEV5aXx3K1L+e01c+no7uXmhzfyyd+vYuX2mqC8otAEcRJWbKkmKSaSs6YGrv7gkZMWT0FGAm/uqh3yPsYYVm6v4fTJoV2Uvbwwm87u3pMaKO1fW6pIionk7GmBf+3UyOZyCZfNy2Llf53Bb66eQ8vxbr700AbOvfttHl59IKhqFJogTlB7Vw8vlxzmglljHfuwPX/mWN7fW09jW+eQti+tbqG6uZ3zZ4Z2m/uc7BQKMhN4cmPlCe3ffLyLf2+t5pJTxoV0olShLTLCxVXzs3n9m2fyu+VzSYyN5PvPlbLkF2/wsxe2U1bb6nSImiBO1Culh2nt6GbZ3PGOxXDJ7HH09LqvCobi1e01uATOCfFvzSLCZxbmsrG88YQmbnlmUyXtXb1cv3iCDdEpNTxRES6Wzc3iuVuX8uQtS1hSkM5f3zvAuXe/w+V/fI9H1x084T5PJ0sTxAl6eHU5eenxjvYlmJWVTE5aHC8MoanFGMMLW6spmpAWtPNPD8fVRTnER0fw1/f2D2s/YwyPrClnbk5q0M+DoUYWEfeskff+x3zW3HkO3714Okfbu7nzmW0U/eQ1bvz7ep7cUEHTsaG1GPhDZMCeKYyUVDWzobyR710yHZfLufF7RISLZ4/jgXf3c+RoBxkDfPAXVzazt66NL55eEMAI7ZMSF8VV87N5bF0Fd1w0jdFJQ+vF/vqOWvbWtXH3p+fYHKFSJy4zKYYvnVHAF0/PZ0tFEy9sPcRLJYd5fWctkS5hcUE6Z03N5IwpmUwenWjbOGJ6BXECfvf6HpJiI7m6KMfpULh6fjbdvYYnNwzcHv/0xkpiIl1ccsq4AEVmv88vzafHGP745t4hbW+M4X/f2ENuWjyfnONc06BSQyUizMsdxfcuncGqb3+CFbct5UtnFHC4pZ2fvrCD8+95h1N/+Qa/e22PLc+vCWKYNh9sZOX2Gm46vSAoRkKdNDqJRflp/HNdeb+3yR3t6OZfW6q4YOZYkmOdj9lf8jMSuHp+Nv9ce5DKxsGHMVi5vYatlc3c+omJRAW434pSJ0tEOCU7lW9fOI3XvnEm799xNr+8YjbzclMx2HOLrL5LhqGrp5fvPFvC6KQYPn9avtPhfOC6xROoaDjOyh2+i9WPrTtIa3s3NwZRzP7ytXMng8DPX9wx4HbtXT385IXtTBqdyBWF2QGKTin7jE+NY/nCXP543Xy+fu4UW55DE8Qw/PzFHew41MJPLptFYkzwlG8unjWWgowE7n51Nz19riKOdnTz53f2sSg/jTkBHA4kUMalxPG1cybz4rbD/Htrdb/b/eaVXVQ0HOfHn5qpVw9KDZGt7xQRuVBEdolImYjc4WN9jIg8bq1fKyJ5XuvutJbvEpEL7IxzML29ht+8sou/vXeALyzN54KZge85PZDICBffPH8qu2pa+Vufu3rufnU3da0d3HHRNIeis9/NZxQwJyeVO57eRmn1x297fW5LFfev2s9nl0zg1BAdg0opJ9iWIEQkAvgDcBEwA7hWRGb02exGoNEYMwm4B/iVte8MYDkwE7gQ+KN1vIBqbOvkha2HuPLe9/n9m2UsX5DDdy+ZHugwhuTi2WM5d/oYfv3yLl7bXoMxhsfWHeSv7+3n+sW5zAuRqUVPRGSEi3uvLyQ5NpLr7l/Lq9b4Np3dvdz/7j7+6/EtLMpPC9rXTqlgJXYNFCUiS4C7jDEXWL/fCWCM+YXXNq9Y26wWkUjgMJAJ3OG9rfd2/T1fUVGR2bBhw7Dj3HGoha88upnunl66egzdvb109xg6e3ppbXd3eR+XEsvt50/lisKsoJ6WsulYJ9fdv5bS6hYyEmM4crSD0ydncP8NRcREhn+P4QNH2rjlkY3sPNxKRmI0HV29tHZ0c96MMfz2mrkkBFGzoFLBQkQ2GmOKfK2z8x2TBVR4/V4JLOpvG2NMt4g0A+nW8jV99v3YqGoichNwE0Bubu4JBRkfHcGUMYlEulxERghRLhdRkUKky0VWahyzslJYmJ9GhIP9HYYqNT6ap245lUfWlLPjUAuFE0ZxzYKcEdPmnpeRwHO3LeX54kOs3VdPTJSLc6eP4cwpmUGd2JUKVnYmCF/vyL6XK/1tM5R9McbcB9wH7iuI4QYIMCE9gT9eN/9Edg1KcdERfOmM8OgMdyJiIiO4an42V83XO5WUOll2frWsBLx7kmUDfW8z+WAbq4kpBWgY4r5KKaVsZGeCWA9MFpF8EYnGXXRe0WebFcAN1uOrgDeMuyiyAlhu3eWUD0wG1tkYq1JKqT5sa2Kyagq3Aa8AEcBfjTGlIvJjYIMxZgXwAPCwiJThvnJYbu1bKiJPANuBbuBWY0yPXbEqpZT6ONvuYgq0E72LSSmlRrKB7mIaGbe3KKWUGjZNEEoppXzSBKGUUsonTRBKKaV8CpsitYjUAeUncYgM4IifwgkFI+18Qc95pNBzHp4JxphMXyvCJkGcLBHZ0F8lPxyNtPMFPeeRQs/Zf7SJSSmllE+aIJRSSvmkCeJD9zkdQICNtPMFPeeRQs/ZT7QGoZRSyie9glBKKeWTJgillFI+jfgEISIXisguESkTkTucjicQROSAiGwTkS0iEpYjHIrIX0WkVkRKvJalichKEdlj/RtWE3X3c853iUiV9VpvEZGLnYzR30QkR0TeFJEdIlIqIl+zloflaz3A+dryOo/oGoSIRAC7gfNwT1K0HrjWGLPd0cBsJiIHgCJjTNh2JhKRM4CjwEPGmFnWsl8DDcaYX1pfBkYZY77tZJz+1M853wUcNcb8xsnY7CIi44BxxphNIpIEbAQuAz5HGL7WA5zvp7HhdR7pVxALgTJjzD5jTCfwGLDM4ZiUHxhj3sE9x4i3ZcCD1uMHcb+xwkY/5xzWjDGHjDGbrMetwA7c89eH5Ws9wPnaYqQniCygwuv3Smz8zw4iBnhVRDaKyE1OBxNAY4wxh8D9RgNGOxxPoNwmIlutJqiwaGrxRUTygHnAWkbAa93nfMGG13mkJwjxsWwktLktNcYUAhcBt1pNEyo8/QmYCMwFDgH/42w49hCRROBp4OvGmBan47Gbj/O15XUe6QmiEsjx+j0bqHYoloAxxlRb/9YCz+JuahsJaqw2XE9bbq3D8djOGFNjjOkxxvQCfyEMX2sRicL9YfkPY8wz1uKwfa19na9dr/NITxDrgckiki8i0bjnxF7hcEy2EpEEq7iFiCQA5wMlA+8VNlYAN1iPbwCeczCWgPB8SFouJ8xeaxER3HPb7zDG3O21Kixf6/7O167XeUTfxQRg3Q72WyAC+Ksx5mcOh2QrESnAfdUAEAn8MxzPWUQeBc7CPQxyDfBD4F/AE0AucBC42hgTNkXdfs75LNzNDgY4ANzsaZsPByJyGvAusA3otRZ/B3e7fNi91gOc77XY8DqP+AShlFLKt5HexKSUUqofmiCUUkr5pAlCKaWUT5oglFJK+aQJQimllE+aIJTyQUR6rFExS0TkeRFJHeb+d4nI7dbjH4vIufZEqpR9NEEo5dtxY8xca1TUBuDWEz2QMeYHxpjX/BWYuOl7V9lO/8iUGtxqrEEcRSRRRF4XkU3WnBofjP4rIt+15hZ5DZjqtfzvInKV9fiAiGRYj4tE5C3r8ZleY/lv9vR29zpGnjUHwB+BTXx0iBilbBHpdABKBTNrzpBzcA9vANAOXG6MabE+6NeIyAqgEPdQLfNwv6824R6rf6huB241xrxnDcTW7mObqcDnjTFfPrGzUWp49ApCKd/iRGQLUA+kASut5QL8XES2Aq/hvrIYA5wOPGuMOWaNrjncMb3eA+4Wka8CqcaYbh/blBtj1pzAuSh1QjRBKOXbcWPMXGACEM2HNYjrgExgvrW+Boi11g1l3JpuPnzfefbDGPNL4ItAHO6rkmk+9m0b7kkodTI0QSg1AGNMM/BV4HZrmOUUoNYY0yUin8CdQADeAS4XkTirfvDJfg55AJhvPb7Ss1BEJhpjthljfgVsAHwlCKUCShOEUoMwxmwGinHXGP4BFInIJDhCGQAAAGdJREFUBtxXEzutbTYBjwNbcI/V/24/h/sR8DsReRfo8Vr+deuW2mLgOPCSHeei1HDoaK5KKaV80isIpZRSPmmCUEop5ZMmCKWUUj5pglBKKeWTJgillFI+aYJQSinlkyYIpZRSPv1/3ifOaRFuWkcAAAAASUVORK5CYII=\n", "text/plain": [ "

" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "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, 2017. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", "\n", "# Hdensity.py: Hydrogen Radial density calling rk4Algor\n", "\n", "import numpy as np\n", "import matplotlib.pylab as plt\n", "\n", "#from rk4Algor import rk4Algor if rk4algor in same directory as this file\n", "\n", "%matplotlib inline\n", "n = 5; el = 2; dr = 0.01 # n = npr+el+1 \n", "rVec = np.zeros((2500),float) # array for plot\n", "RhoVec = np.zeros((2500),float) # Density array\n", "fvector = np.zeros(2) \n", "y = np.zeros(2); y[0] = 1e-8; y[1] = 0 \n", " \n", "def f(r,y): # RHS of ODE \n", " fvector[0] = y[1]\n", " fvector[1] = -(2/r-1)*y[1]-((n-1)/r-el*(el+1)/r**2)*y[0]\n", " return fvector\n", "\n", "f(0.001,y) # Integration, f(t= 0)\n", "i = 0\n", "\n", "def rk4Algor(t, h, N, y, f):\n", " k1=np.zeros(N); k2=np.zeros(N); k3=np.zeros(N); k4=np.zeros(N);\n", " k1 = h*f(t,y) \n", " k2 = h*f(t+h/2.,y+k1/2.)\n", " k3 = h*f(t+h/2.,y+k2/2.)\n", " k4 = h*f(t+h,y+k3)\n", " y = y+(k1+2*(k2+k3)+k4)/6.\n", " return y \n", "\n", "for r in np.arange(0.001,25,dr):\n", " rVec[i] = r\n", " y = rk4Algor(r, dr, 2, y, f) # call rk4 algorithm \n", " RhoVec[i] = 4*3.141593*(y[0]*np.exp(-0.5*r) )**2 *r**2 \n", " i = i+1\n", "plt.figure()\n", "plt.plot(rVec,RhoVec)\n", "plt.title('Hydrogen Radial Density n = 5')\n", "plt.xlabel('Radius r')\n", "plt.ylabel('Density')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [] } ], "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 }