Skip to content

Instantly share code, notes, and snippets.

@berquist
Created March 10, 2016 06:04
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save berquist/7e3855c7b7d4acafe1b5 to your computer and use it in GitHub Desktop.
Save berquist/7e3855c7b7d4acafe1b5 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Explaining the Relative Energies of Various Vibrational Modes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"http://chemistry.stackexchange.com/questions/45138/explaining-the-relative-energies-of-various-vibrational-modes"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import sympy as sy"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Vendor: Continuum Analytics, Inc.\n",
"Package: mkl\n",
"Message: trial mode expires in 27 days\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"IPython console for SymPy 1.0 (Python 3.5.1-64-bit) (ground types: python)\n",
"\n",
"These commands were executed:\n",
">>> from __future__ import division\n",
">>> from sympy import *\n",
">>> x, y, z, t = symbols('x y z t')\n",
">>> k, m, n = symbols('k m n', integer=True)\n",
">>> f, g, h = symbols('f g h', cls=Function)\n",
">>> init_printing()\n",
"\n",
"Documentation can be found at http://docs.sympy.org/1.0/\n"
]
}
],
"source": [
"sy.init_session()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# help(symbols)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"alpha, beta, hbar = symbols('alpha, beta, hbar', real=True)\n",
"omega_1, omega_2, omega_a, omega_s = symbols('omega_1, omega_2, omega_a, omega_s', real=True)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAADsAAAAyBAMAAAAQBL5IAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhB2mUQi3bvN\nZqsoIwvDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABdElEQVQ4Ec3Ur0/DQBQH8O+W/pooNICZYjSh\nlmnUEoZAkExhWYJAY1AkTOAww5LgkRMkKJJq1P6BJfsLCCFBd7fee+8O1isExZlr3+cu7aXfV2wX\nH6gezaJIsNk/rFb4/YMEW9XYmo0Bz8ln/lMNxz3c1vADsF/DQ+CtlqO24ex89OUA0dHl60D4fuS3\nW11rQeCnFz3hDeA5UKtlNIEgZw4/gStVMSMAonfmxgTY2TUI3ABxwnzaAfa6NqvC2px5farYVuTA\n3Zg5nCO+9gbWghd4j+qWPkk6m4aZpdFxli5XO74Yn+JvfEIPcuxWxy6Hg0ldr/avuEw9v5HMcrAy\n9VLmC2adeq7KzKxTL2W+YB6WqeeqzIaXqV8ZxJR6i3VfEFPqDVNfEFPqDVNf8G6demHuC2JKvXBj\novuCuKNTL8x9QZzr1AtzXxBT6oVD6gvNnHphUF9o5tQbpqtfMafesZtT7+CVMhfoYHz7fVb8ww+7\n/ne/ADnVTA2yNxppAAAAAElFTkSuQmCC\n",
"text/latex": [
"$$\\left[\\begin{matrix}\\alpha & \\beta\\\\\\beta & \\alpha\\end{matrix}\\right]$$"
],
"text/plain": [
"⎡α β⎤\n",
"⎢ ⎥\n",
"⎣β α⎦"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# hamiltonian = Matrix([[hbar * omega_1, beta],\n",
"# [beta, hbar * omega_2]])\n",
"hamiltonian = Matrix([[alpha, beta],\n",
" [beta, alpha]])\n",
"hamiltonian"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [],
"source": [
"normal_mode_eig = hamiltonian.eigenvects()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"normal_mode_eigenvals = [x[0] for x in normal_mode_eig]\n",
"normal_mode_eigenvects = [x[2][0] for x in normal_mode_eig]"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAI0AAAAUBAMAAABCGaMdAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdt3NMolEmSIQ71S7\nZqvC3fgDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABp0lEQVQ4EaWTMUvDQBiG32CbJtZoaREHF6ki\nuJXUHxCcXIR2qKKLGTqJIOLmYnARhEId3AQFFToWf0EHwUnIoCBCob+gCE5ufte7s3chaQt+wyXf\ne+/7cHe5IL+O+Gp8xOtR1eoGjutjO6qL3mjfeAlTurxoPyOdzHmD2dED8Z3j4XAExwqRXYpPDtRs\nTkw2gKcRHAOY7U3CIV5/NOcgmIxjtiSnuuBHI/WLwlpUA4ZGuS9z89htCs65b7estp7KXe886Ap1\nilFyDLt26gnOI/BiNPVUBSjoCnWKUXLSgFHinMwPcEI9bLdIteyxPB3fEXuq9WdMl8vufbnco0n6\nHuYX50x3gL1dNUBzHvCpS4BqlOupA84d59zSHvbbeogtr69LgGqUHMrS/Rjc57mQOJHMFGC1IhpU\no+SUgMtAnE8PzlmqqaVmQmz4sL81MaMYJecVqVVwDmrdMFPVEqi/z7MTK+p0xSg45la1Rp7E/7TC\nsdlAxw87wWHnSDWOI2zcrI12MGjHcOwcD11p2ZhGGJLWw64Pq5A/kke6PqyIk1/hr/8ZjaL/Cxe/\nWe2emOs9AAAAAElFTkSuQmCC\n",
"text/latex": [
"$$\\left [ \\alpha - \\beta, \\quad \\alpha + \\beta\\right ]$$"
],
"text/plain": [
"[α - β, α + β]"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"normal_mode_eigenvals"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAH4AAAAyBAMAAABol3KsAAAAKlBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADmU0mKAAAADXRSTlMA74lUMhDN3auZdmYiZQh4sgAA\nAAlwSFlzAAAOxAAADsQBlSsOGwAAASVJREFUSA1jkL17iYGB8e5dBQYkwHT3rgASF8EEqwZyIQp0\n715kEHZxBeo3cQlAKGJgYHVxxq4frBqoEqIgxMWQQQSkj9EBRCIDFuz6warB6iAKHGmgn1kUaAE+\n+w1BfsVpf8SMq3BpsEuRCaD7WTt78epn4MKvn4Fh7ah+LOEfVg4CKQyj4UdZ+uFsvNYNC17kpAtm\ng/LfDNmMDTAFtMh/YIvw5T8kBaP2D77yEyl6wExkgh7lNxXtZ72G7Hggm5D7wXUHSA80/ZaDimok\nANaPu/6A1B0g9VD9zAeQNAOZQP146w9w3QHSAtXPhKod4n489Qe6/igK9QMLChQA9j8J9qNoBnJG\n9QMDYfCGH6TuAEUanvJ7tP5ABNDA11+Utt8p7D8AALdkfj2ULWJSAAAAAElFTkSuQmCC\n",
"text/latex": [
"$$\\left [ \\left[\\begin{matrix}-1\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\1\\end{matrix}\\right]\\right ]$$"
],
"text/plain": [
"⎡⎡-1⎤, ⎡1⎤⎤\n",
"⎢⎢ ⎥ ⎢ ⎥⎥\n",
"⎣⎣1 ⎦ ⎣1⎦⎦"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"normal_mode_eigenvects"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"# from IPython.display import set_matplotlib_formats\n",
"# set_matplotlib_formats('pdf', 'svg')"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"beta_plot = -np.linspace(start=0, stop=500)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"eigval_a = normal_mode_eigenvals[0]\n",
"eigval_s = normal_mode_eigenvals[1]"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"local_mode_frequency = 1750.0"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"eigval_a_plot = [eigval_a.subs(alpha, local_mode_frequency).subs(beta, x) for x in beta_plot]\n",
"eigval_s_plot = [eigval_s.subs(alpha, local_mode_frequency).subs(beta, x) for x in beta_plot]"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEVCAYAAADU/lMpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VfWd//HXJwiKAmWRRdkCQQu4gIholda4QILtKAxt\nxVLQoVWrVdHitFj9jdSpVZl2VNpHO/VXBXEpRZ1xqTQBq9FW2rIICgJWoqwKyk9ri4Msyef3xzlJ\nbmISksu599zl/Xw88vCe7z33nu89N+HtOd/N3B0REZFUK4i7AiIikh8UOCIikhYKHBERSQsFjoiI\npIUCR0RE0kKBIyIiaZH2wDGzPmb2vJm9bmZrzOy6Bs/PMLNqM+uaUHaTmb1pZuvNbGxC+Qgze83M\n/mpm96Tzc4iISOvEcYVzAPiOu58AfA74tpkNhiCMgDHA5pqdzWwI8FVgCDAO+LmZWfj0L4BvuPvx\nwPFmVpK+jyEiIq2R9sBx9x3uvjp8vBtYD/QOn74b+NcGL7kIWODuB9x9E/AmMMrMegEd3X15uN98\nYHyq6y8iIsmJtQ3HzAqB4cBfzOxCYKu7r2mwW29ga8L29rCsN7AtoXwbdcElIiIZ5rC4DmxmHYDH\ngelAFfB9gttpIiKSg2IJHDM7jCBsHnL3p8zsRKAQeDVsn+kDvGJmowiuaPolvLxPWLYd6NtIeWPH\n04RxIiJJcHc7+F4tE9cttQeAde5+L4C7r3X3Xu4+0N0HENweO8Xd3wOeBi42s3ZmNgAYBCxz9x3A\nR2Y2KgypqcBTTR3Q3fXjzq233hp7HTLlR+dC50LnovmfqKX9CsfMzgImA2vMbBXgwPfdvSxhNwcM\nwN3XmdlCYB2wH7ja687Et4F5wBHAogbvISIiGSTtgePuLwNtDrLPwAbbdwB3NLLfSuCkSCsoIiIp\noZkG8kxxcXHcVcgYOhd1dC7q6FykjqXiPl2mMTPPh88pIhIlM8NzoNNAxliwYAFFRUW0adMGM9NP\nFv20adOGoqIiFixYEPevkYi0QGzjcDLBggULuOGGG3jiiScYOXIk7dq1i7tK0gr79u1jxYoVTJw4\nEYBJkybFXCMRaU5e31IrKirioYce4swzz4yhVhKVpUuXMmXKFCorK+OuikhOifqWWl4HTps2bdiz\nZ4+ubLLcvn37aN++PVVVVXFXRSSnKHCS0FTghCczhhpJ1PRdikRPnQYkZ/3xj39kyJAhcVdDRFJE\nVzh58PkzQUFBARs3bmTgwIEH3zkJ+i5FoqcrHMlKZs3/zqr9RST3KXAy2F133cWgQYPo1KkTJ554\nIk8++SQAlZWVFBcX07lzZ3r06MEll1wCwDXXXMONN95Y7z0uuugi7r33XgAGDBjAj3/8Y4YNG0bH\njh25/PLLee+997jgggvo1KkTY8eO5aOPPgJg8+bNFBQUMG/ePPr160e3bt345S9/yYoVKxg2bBhd\nu3bl2muvrXesBx54gKFDh9KtWzfGjRvH1q3BMkZnn3027s7JJ59Mp06deOyxx3jxxRfp27cvs2fP\n5phjjmHatGm1ZTW2bdvGxIkT6dGjB927d+e66+qtRi4i2Sbu2UjTNOOpN6ap8kzx+OOP+44dO9zd\nfeHChd6hQwffsWOHX3LJJf6jH/3I3d337t3rL7/8sru7L1u2zHv37l37+l27dvlRRx3l77//vru7\nFxYW+uc+9zl///33/Z133vEePXr4qaee6q+++qrv3bvXzz33XL/tttvc3X3Tpk1uZn7VVVf53r17\nfcmSJX7EEUf4hAkTfNeuXb59+3bv0aOHv/TSS+7u/uSTT/pxxx3nb7zxhldVVfntt9/uZ555Zm1d\nzMzfeuut2u2Kigo/7LDD/KabbvJ9+/b5J5984hUVFd63b193d6+qqvJhw4b5jBkzfM+ePfU+Z2My\n/bsUyUbh31V0/xZH+WaZ+pOtgdPQ8OHD/emnn/ZLL73Ur7zySt+2bdun9hk6dKg/99xz7u7+s5/9\nzL/4xS/WPldYWOiPPvpo7fbEiRP96quvrt3+6U9/6hMmTHD3IHAKCgr83XffrX2+W7duvnDhwnqv\nv/fee93dfdy4cf7AAw/UPldVVeVHHnmkb9myxd2DwKmsrKx9vqKiwg8//HDft29fvbKawFm6dKn3\n6NHDq6qqWnRusu27FMkGUQeObqkdhNmh/yRr/vz5nHLKKXTp0oUuXbrw+uuvs2vXLmbPnk11dTWj\nRo3ipJNOYu7cubWvmTp1Kg8//DAADz/8MFOmTKn3nj179qx93L59+09t7969u97+PXr0aNH+mzdv\nZvr06XTt2pWuXbvSrVs3zIzt2xtdEw+A7t2707Zt20af27ZtG/3796egQL+iIrkir6e2aQmPqePT\nli1buOKKK3jhhRf43Oc+B8App5yCu9OjRw/uu+8+AF5++WXOP/98zj77bAYOHMjXv/51TjrpJF57\n7TU2bNjA+PHj01Lfvn37csstt9S2J7VEcx0J+vbty5YtW6iurlboiOQI/SVnqI8//piCggKOPvpo\nqqurmTt3LmvXrgXg8ccfr71y6Ny5MwUFBbX/KPfu3ZuRI0cyZcoUJk6cyOGHH550HbwVafutb32L\nH/3oR6xbtw6Ajz76iMcff7z2+V69evHWW2+1+P1GjRrFMcccw8yZM/nf//1f9u7dy9KlS1teeRHJ\nOAqcDDVkyBBmzJjBGWecQa9evXj99dcZPXo0AMuXL+f000+nU6dOjB8/njlz5lBYWFj72ksvvZS1\na9cyderUeu/Z8IriYF2VD7Z/4vb48eOZOXMmkyZNonPnzpx88smUldUtwDpr1iymTp1K165d6wVR\nUwoKCnjmmWd488036devH3379mXhwoUHfZ2IZC4N/MzBz/+HP/yBKVOmsGnTprirkja5+l2KxEkD\nP6VZ+/fv59577+Xyyy+PuyoiIvUocHLIhg0b6NKlCzt37mT69OlxV0dEpB7dUsuDz58P9F2KRE+3\n1EREJCspcEREJC008FNEROr58EP4/e+jf18FjohInquqguXLobw8+Fm7FsJhf5FSp4E8+Pz5QN+l\nSOts21YXML//PfTuDSUlwc/o0XDEEdF3GlDg5MHnzyQXXHABl1xyyacmFT1U+i5FmrdnD7z0Ul3I\n7NwJ559fFzLHHvvp1yhwkqDASb0f/OAHVFZWMn/+/FiOr+9SpD53WL++LmBefhmGDasLmFNPhTZt\nmn+PqAMn7W04ZtYHmA/0BKqB+9z9p2Y2G/gnYC9QCfyLu/89fM1NwDTgADDd3ReH5SOAecARwCJ3\nvz7NH0dawd0POn+biCTvww/huefqQqagIAiXyy+HBQugc+eYKxjl4jot+QF6AcPDxx2AN4DBwPlA\nQVh+J3BH+HgosIogHAuBjdRdmf0FOC18vAgoaeKYzS0ulLHuvPNO7927t3fs2NEHDx7sjzzyiB95\n5JH+wQcf1O6zcuVK7969ux84cMDnzZvnZ511lt9www3euXNnLyoq8qVLl/q8efO8b9++3rNnT3/w\nwQdrX3vZZZf51Vdf7ePGjfMOHTr46NGjfceOHX799dd7ly5dfMiQIb569era/d955x2fOHGid+/e\n3QcOHOhz5sxxd/eysjJv166dt2vXzjt06ODDhw93d/fi4mK/+eab/ayzzvIjjzzSKysrvbi42O+/\n//7a97zvvvt8yJAh3rFjRz/hhBN81apVSZ2rTP8uRVLhwAH3P/3J/dZb3c84w71DB/dx49zvucd9\n/Xr36upDe39ybcVP4EngvAZl44GHwsczge8lPPc74PQwuNYllE8CftHEMZo7mRnpjTfe8L59+9Yu\nMb1582avrKz0L37xi/5f//VftfvdcMMNft1117m7+7x587xt27b+4IMPenV1td9yyy3er18/v+aa\na3zfvn2+ePFi79ixo3/88cfuHgRO9+7dfdWqVbVLTA8YMMAffvjh2tefc8457u5eXV3tp556qv/w\nhz/0AwcO+Ntvv+1FRUW+ePFid3efNWuWT5kypd5nKC4u9v79+/v69eu9qqrK9+/fXy9wFi5c6H36\n9PGVK1e6u3tlZWXtCqGtlcnfpUiUtm51/9Wv3L/yFfeuXd1PPNF9xgz3xYvd9+yJ9lhRB06sAz/N\nrBAYTnClkmgawRULQG9ga8Jz28Oy3sC2hPJtYVlOaNOmDfv27WPt2rUcOHCAfv36MXDgQKZOncpD\nDz0EQHV1Nb/+9a/rLUMwYMAApk6diplx8cUXs23bNm699Vbatm3LmDFjaNeuHRs3bqzdf8KECQwf\nPpx27doxYcIE2rdvz+TJk2tfv3r1agCWLVvGrl27uPnmm2nTpg2FhYV885vfZMGCBc1+jssuu4zB\ngwdTUFDAYYfVv4N7//33893vfpcRI0YAMHDgQPr27RvJ+RPJFXv2BLfHvvMdOOEEGD4cliyBcePg\ntddgzRr48Y9hzJigZ1kmi20cjpl1AB4naJPZnVB+M7Df3X8dV90S2Q8Ovc3Bb219Y3ZRURH33HMP\ns2bNYt26dZSUlPCf//mfXHTRRVx11VVs3ryZ9evX07lzZ0499dTa1zVcAhrg6KOPrleWuIx0S5ec\n3rJlC9u3b6dr167BZ3KnurqaL3zhC81+juYCZOvWrRQVFTX7epF844009p98MpSWwty5LWvsz1Sx\nBI6ZHUYQNg+5+1MJ5ZcBFwDnJuy+HUj8V6tPWNZUeaNmzZpV+7i4uJji4uIW1TWZsIjKpEmTmDRp\nErt37+aKK67ge9/7Hg8++CBf/epXeeihh9iwYUPk3Yub0rdvXwYOHMgbb7zR6PNNdQY42DLSlZWV\nkdRPJJslNvYvXgxm8TT2V1RUUFFRkbL3j+sK5wGC9pd7awrMrBT4V+AL7r43Yd+ngUfM7G6CW2aD\ngGXu7mb2kZmNApYDU4E5TR0wMXCywV//+le2b9/OWWedRbt27Wjfvj3V1dUATJkyhalTp/L+++9z\nxx13NPs+wW3Y5NW8ftSoUXTs2JHZs2dz3XXX0bZtWzZs2MCePXsYOXIkPXv25LnnnmtVT7RvfvOb\nzJgxg7POOosRI0ZQWVlJ27Zt6dev3yHVWSTTNTWyv6QEbrwRPvvZIHTSreH/jP/gBz+I9P3T3oZj\nZmcBk4FzzWyVmb1iZuOAnxL0WlsSlv0cwN3XAQuBdQTtOld73b+i3wbuB/4KvOnuZeSIvXv3MnPm\nTLp3786xxx5bL1zOPPNMCgoKGDFixEHbPFq7rHRTry8oKOC3v/0tq1evZsCAAfTo0YPLL7+cv//9\n7wB85Stfwd3p1q0bI0eObPJYiWVf/vKXufnmm/na175Gp06dmDBhAh9++GGr6ieSLbZtg/vvh69+\nFXr0gCuugN274bbb4L33YNEimD4dBg+OJ2zSQQM/s/Tzn3feeUyePJlp06bFXZWMkM3fpeSmTz6p\nP7J/x46gYb+kBMaObXxkf6bRTANJyLXAWb58OSUlJWzdupWjjjoq7upkhGz9LiV3NGzsX7o0aOwv\nKQka/EeMyL7G/qyfaUAOzWWXXcZTTz3FnDlzFDYiMauZxr/hyP4rrsiQkf0ZRlc4efD584G+S0mH\n5hr7S0vh+ONzq/1Ft9SSoMDJffouJVVaMo1/rlLgJKGpwGnTpg179uyhXbt2MdRKorJv3z7at29P\nVVVV3FWRHNDUNP6lpdnT2B+VjGjDMbOjgE/cPav/wgsLC1mxYgVnnnlm3FWRQ7BixQoKCwvjroZk\nqcTG/rKyoLG/Zhr/bB/Zn2laFDhmVkAwOeZk4DSCJQQON7NdwLPAL919YzNvkZFuv/12Jk6cyBNP\nPMHIkSN1pZNl9u3bx4oVK5g4cSJ333133NWRLNLUNP5XXgm/+Y0a+1OlRbfUzOxF4DngKWCtu1eH\n5V2Bc4CvAf/j7g+nsK5Ja+qWGsCCBQu4+eab2bRpU+1IfskOBQUFFBYWcvvttzNp0qS4qyMZrKax\nv6ysrrH/85+va4uJa2R/poulDcfM2rr7/kPdJy7NBY6I5KaGjf3HHlvXmyzXG/ujok4DSVDgiOS+\nPXvgD3+ou4qpaeyvGdnfO2cWL0mftAeOmR0DJB7w3Ey9ddYUBY5I7mluGv+SEjX2RyGOwLkQuAxY\nTRA8x7v75KgqkA4KHJHc0FRjf0kJnHeeGvujFlcbTk933xk+7uHu70VVgXRQ4Ihkp+ZG9quxP/Vi\nbcMxszPc/c9RHTxdFDgi2SOfR/ZnmrgHfn4mqgOLiEBdY3/NwMuaxv5x4+Duu9XYn0taGzi6TBCR\nQ9LcNP4a2Z/bWhs4ulsqIq2mafwFWt+Gc6y7v5PC+qSE2nBE0ivfpvHPVRr4mQQFjkjqqbE/98Td\nS20kcDPQn+B2nAHu7idHVaFUUOCIRE/T+Oe+uAPnDeBfgTVA7UyX7r45qgqlggJH5NA1N42/Rvbn\nprgD54/uPjqqg6eLAkckOU2N7C8thXPPVWN/ros7cM4DLgF+T7AmDgDu/t9RVSgVFDgiLZPY2F9W\nBq+/rpH9+SzuwHkYGAy8Tt0tNXf3aVFVKBUUOCJNU2O/NCXuwHnD3T8b1cHTRYEjUqepkf2axl8a\nintqm6VmNtTd10VVARFJrcam8a9p7NfIfkmn1l7hrAeKgLcJ2nDULVokA6mxX6IQ9y21/o2Vq1u0\nSLw0jb+kgmYaSIICR3KRGvsl1aIOnIJWHvxBM+ucsN3FzB5o5Xv0MbPnzex1M1tjZtclvNdiM3vD\nzMrN7DMJr7nJzN40s/VmNjahfISZvWZmfzWze1pTD5Fss2cPLF4MM2bACSfA8OHBbbMLLoA1a+C1\n1+A//iPoAKCwkUzU2ltqq9z9lIOVHeQ9egG93H21mXUAVgIXAf8C/D93n21m3wO6uPtMMxsKPAKc\nBvQBngOOc3c3s78A17j7cjNbBNzr7uWNHFNXOJJ1mpvGXyP7JR3i7qVWYGZd3P3DsDJdW/se7r4D\n2BE+3h12ROhDEDpnh7s9CFQAM4ELgQXufgDYZGZvAqPMbDPQ0d2Xh6+ZD4wHPhU4ItlC0/hLLmtt\n4PwE+JOZPRZufwW4PdmDm1khMBz4M9DT3XdCEEpm1iPcrTfwp4SXbQ/LDgDbEsq3heUiWaO5xv4b\nb9Q0/pJbWnt1Mt/MVgDnhkX/nOyYnPB22uPA9PBKp+E9L90Dk5zUVGP/bbepsV9yW4sCxxIaQcKA\n+VTIWCsaSszsMIKwecjdnwqLd5pZT3ffGbbzvBeWbwf6Jry8T1jWVHmjZs2aVfu4uLiY4uLillRV\n5JA1No3/mDFBY/8992gaf8kcFRUVVFRUpOz9W9RpwMwqgCeAp9x9S0J5O2A0cCnwgrvPa9FBzeYD\nu9z9OwlldwEfuPtdTXQaOJ3gltkS6joN/Bm4DlgOPAvMcfeyRo6nTgOSNprGX3JFLONwzOwIYBow\nGRgA/A1oT9CtejHwc3df1aIDmp0FvESwpo6HP98HlgELCa5aNgNfdfe/ha+5CfgGsJ/gFtzisPxU\nYB5wBLDI3ac3cUwFjqRUYyP7S0uDgNHIfslWsQ/8NLO2wNHAnppAyHQKHImapvGXfBB74GQjBY5E\noWFj/7HH1l3FqLFfcpECJwkKHElGY439msZf8okCJwkKHGmJ5qbxV2O/5KO4Z4u+Fni4ZqaBbKHA\nkaZoGn+RpsU9tU1PYLmZvQI8AJTrX3LJJs019t94oxr7RVIpmV5qBowlmGxzJEFX5vvdvTL66kVD\nVzj5TdP4iyQn7iscwgGXNRNwHgC6AI+b2RJ3/25UFRNJ1p498Ic/1F3FaGS/SGZobRvOdGAqsAv4\nFfCku+83swLgTXcvSk01D42ucHJbU9P413RZHjFCjf0iyYj7CqcrwYSd9ZaUdvdqM/tSVJUSORhN\n4y+SfVobOIXARzUbZtYF+Im7T3P39VFWTCRRY9P4f/7zmsZfJJukfcXPOOiWWnZSY79IvOK+pXbI\nK36KNEXT+IvktlhX/JT85g7r1tVv7K8Z2T93rkb2i+SaZMbhDKVuxc/nk13xM510Sy1zaBp/keyh\nudSSoMCJT1UVLFtWFzCaxl8ke8Q9l9rhwESC3mq1t+Pc/baoKpQKCpz0UmO/SG6Iu9PAUwTdolcC\ne6OqhGS3pqbxHzcO7r5b0/iLSKC1Vzhr3f3EFNYnJXSFE63GpvFPHNmvxn6R3BD3Fc5SMzvJ3ddE\nVQHJDk1N43/55RrZLyIt09ornHXAccBbBLfUjGA+z5NTU71o6Aqn9RpO4792bdD+UnMVo8Z+kdwX\nd6eB/o2VN5xbLdMocFqmYWP/scfWNfZ//vNq7BfJN3EHjgGTgYHufpuZ9QN6ufuyqCqUCgqcxjXV\n2F9SAmPHqrFfJN/FHTi/AKqBc919SDh552J3Py2qCqWCAifQWGN/zch+NfaLSENxdxo43d1HmNkq\nAHf/0MzaRVUZiZ4a+0UkU7Q2cPabWRvAAcysO8EVj2SIxqbxrxnZf+ONauwXkfi09pbaZOBi4FRg\nHvBl4BZ3f6y518Ut12+paWS/iKRC7HOpmdlg4Lxw8/lsWHgt1wKnqWn8axr7NY2/iEQh7k4D/9ZY\nueZSS63mpvFXY7+IpErcnQY+Tnh8BPAloNVXOGZ2f/janTWDRs1sGPBf4fvuB6529xXhczcB04AD\nwHR3XxyWjyC4tXcEsMjdr29tXTJVU9P4X3kl/OY3auwXkexzSMsThLNHl7t7cStfNxrYDcxPCJxy\n4CfuvtjMxgHfdfdzwvV3HgFOA/oAzwHHubub2V+Aa9x9uZktAu519/JGjpfxVzgNp/FfuzYYbKlp\n/EUkLnFf4TR0JEEItIq7/7GRWQuqgc+EjzsD28PHFwIL3P0AsMnM3gRGmdlmoKO7Lw/3mw+MBz4V\nOJlq69ZPN/aXlsK//7sa+0Uk97QqcMxsDWGXaKAN0B2Iqv3mBqDczH5CMEfbmWF5b+BPCfttD8sO\nANsSyreF5RmrqZH9F1wA99yjkf0ikttae4XzpYTHBwjaYA5EVJerCNpnnjSzLwMPAGMieu9YNDeN\n/9y5auwXkfzSqsBJ8SSdl7r79PA4j5vZr8Ly7UDfhP36hGVNlTdq1qxZtY+Li4spLi6OpNINaWS/\niGSriooKKioqUvb+re0W/Z3mnnf3/2zFexUCz7j7SeH26wQ90140s/OAO939tIROA6cT3DJbQl2n\ngT8D1wHLgWeBOe5e1sixUtZpoGZkf1lZEDCvv143sl+N/SKSzeLuNDCSoLfY0+H2PwHLgDdb8yZm\n9ihQDHQzsy3ArcDlwJxw6pxPgCsA3H2dmS0E1lHXXbomPb5N/W7RnwqbVGhqZL8a+0VEmtbaK5yX\ngC+6+z/C7Y7As+7+hRTVLxKHeoXTVGN/aalG9otI7or7CqcnsC9he19YllOam8Zfjf0iIslpbeDM\nB5aZ2f+E2+OBB6OtUjyaGtl/xRVq7BcRiUIyk3eOAD4fbr7k7qsir1XEGrul1tw0/qWlcPzxauwX\nkfwW9+SdWb3EdMPG/mOPrQsYNfaLiNQXd+Bk7RLTxw/fxf/b2o3zz6+bxl8j+0VEmhZ3p4GsXWJ6\n28QBDO05mOOKSjhuUAk9jzmDQ59KTkREWipvlpj+8KZdvLzlZcory7n2d9ey6W+bOKfwHEoHlVJS\nVEL/zg3nEhURkSglu8T0CILeaVm7xPSO3TtYUrmE8spyFlcupmv7rpQUlVAyqISz+5/NUe2Oiqm2\nIiKZIbY2nLDDQB/gKIIlpg34fS4sMV3t1azesZqyjWWUV5bzyruvcHrv02sD6KQeJ2HqsiYieSbu\nTgNrauY+yyatnWng73v/zgtvv0B5ZTnlleXs2b+HsUVjKR1UypiBY+h2ZLcU1lZEJDPEHTgPAj9L\nWPQsKxzq1DYbP9hI+cYgfCo2VTD46MG1Vz9n9DmDwwrU+UBEck/cgbMBGARsBj4muK3mNctEZ6oo\nZ4veV7WvtvNBeWV5beeDmgAq7FwYyXFEROIWS+CY2UPuPsXMrgf+p+HzKV4n55ClcnkCdT4QkVwV\nV+CsA84HfkewrEC9Crj7B1FVKBVSGTiJqr2aVe+uqr36Sex8UDqolBN7nKjOByKSNeIKnOsIloAe\nSLCqZmIF3N0HRlWhVEhX4DSU2PmgbGMZe6v2MrZoLCVFJep8ICIZL+42nF+4+1VRHTxd4gqcRO7O\nxg82srhyMWWVZby46UV1PhCRjBZr4GSrTAichmo6H9QEkGY+EJFMo8BJQiYGTkM7d+9kceXiRjsf\nFBcWc2TbI+OuoojkGQVOErIhcBLVzHxQM/Zn5bsrNfOBiKSdAicJ2RY4DTU180FJUQljisZw9JFH\nx11FEclBcXca+E4jxR8BK919dVSVilq2B05DiTMfvLj5RT7b7bPqfCAikYs7cB4FRgLPhEVfAl4D\nCoHH3H12VBWLUq4FTqJ9VftYunVp7cSjiTMflA4qVecDEUla3IHzEnCBu+8OtzsAzwKlBFc5Q6Oq\nWJRyOXAa0swHIhKVuANnA3CSu+8Ptw8HXnX3wWa2yt1PiapiUcqnwEnUXOcDzXwgIgcTd+D8H2AC\n8BTBbAP/FD7+CXCfu0+OqmJRytfAaahh54NPDnyimQ9EpEmx91Izs5HAWQTLTC919xVRVSZVFDiN\n07ILItKcuK9wDgcmEnQSqP3XyN1vi6pCqaDAOTgtuyAiDcUdOGWE3aCBqppyd/9JVBVKBQVO6zXs\nfNClfRdKi0rV+UAkj8QdOGvd/cRDPqjZ/QRdqncmLt5mZtcCVwMHgGfdfWZYfhMwLSyf7u6Lw/IR\nwDzgCGCRu1/fxPEUOIeguWUXNPOBSO6KO3DuA37q7msO6aBmo4HdwPyawDGzYuD7BN2uD5jZ0e6+\ny8yGAI8CpwF9gOeA49zdzewvwDXuvtzMFgH3unt5I8dT4ESo4bILnxz4hJJBJep8IJJj4g6cdQRL\nTL8N7OUQlpg2s/7AMwmB8xvgl+7+fIP9ZobHuCvc/h0wi2CZ6+drxv6Y2STg7MaWT1DgpE7Nsgs1\nVz9adkEkd0QdOK39l2BcVAduxPHAF8zsR8Ae4EZ3Xwn0Bv6UsN/2sOwAsC2hfFtYLmlkZhzX7TiO\n63Yc14zhOl6BAAAO/UlEQVS6pl7ng2t/d62WXRCRWq0KHHffnKqKENSli7ufYWanAY8RrDAqWaRd\nm3acM+AczhlwDneef2e9ZRduef4WLbsgksdaFDhm9kd3H21m/yAYf1P7FMHtrk4R1GUr8N8Eb7jc\nzKrMrBvBFU2/hP36hGXbgb6NlDdq1qxZtY+Li4spLi6OoMpyMD079GTKsClMGTal3swHs1+ezcWP\nX8zpvU+vvfrRzAci8aqoqKCioiJl7x/b8gRmVkjQhnNSuH0F0NvdbzWz44El7t7fzIYCjwCnE9wy\nW0Jdp4E/A9cBywnmdJvj7mWNHEttOBnoH3v/wQubXqideFQzH4hklrg7Ddzl7t87WFkL3udRoBjo\nBuwEbgUeAuYCwwk6JMxw9xfD/W8CvgHsp3636FOp3y16ehPHU+BkAS27IJJZ4g6cV9x9RIOy15Lp\npZZOCpzsU7PsQvnGcsoqyzTzgUgMYgkcM7uKYEDmQKAy4amOwMvu/vWoKpQKCpzsl9j5QDMfiKRH\nXIHzGaALcAcwM+Gpf7j7B1FVJlUUOLklsfNBWWWZZj4QSZHYZ4vORgqc3NZw2YU9+/cwtmgspYNK\n1flA5BDEdYWT2B264cGj6hadMgqc/KJlF0SioSucJChw8peWXRBJXty91P6tsXKthyPZouGyC4kz\nH6jzgUh9cQfOjITNIwiWGFjv7tOiqlAqKHCkMVp2QaR5GXVLLVwBtNzdi6OqUCoocKQltOyCSH2Z\nFjhdgOXuPiiqCqWCAkdaS8suiMR/S20Ndb3V2gDdgdvc/WdRVSgVFDhyqBrrfHDugHODANKyC5Kj\n4g6cxL+qAwRLRB+IqjKposCRqO3cvZMlbwWdD8o3lmvZBclJcQfOV4Ayd/+Hmd0CjAB+6O6vRFWh\nVFDgSColznxQXlnOyndXatkFyQlxB85r7n6ymY0Gfgj8B/Bv7n56VBVKBQWOpFPNsgs1U+9o2QXJ\nVnEHzip3P8XM7gDWuPujNWVRVSgVFDgSJy27INkq7sD5LcGqmmOBU4A9wDJ3HxZVhVJBgSOZInHZ\nhfLKct7+29ua+UAyVtyBcyRQSnB186aZHQOcVLMgWqZS4EimamzZhZqeb8WFxZr5QGIVd+AY8HVg\ngLvfZmb9gF7uviyqCqWCAkeygZZdkEwTd+D8AqgGznX3IeHAz8XuflpUFUoFBY5kIy27IHGLO3Be\ncfcRiR0FzOxVteGIpJ6WXZB0iztw/gKcSTCdzQgz605whaNeaiJppGUXJB3iDpzJwMUEAz4fBL4M\n3OLuj0VVoVRQ4Eiu07ILkgqxBU7YYaAPcBRwHsHKn7939/VRVSZVFDiST2o6H5RtLNOyC3JI4r7C\nWePuJ0V18HRR4Eg+a6zzgZZdkJaIO3AeBH7m7sujqkA6KHBE6tR0PiirLNOyC9KsuANnAzAI2Ax8\nTHBbzd395KgqlAoKHJHGadkFaU7cgdPob5+7b46qQqmgwBFpmZqZDxa/tVjLLkhmrfiZLRQ4Iq3X\n1LILJUUllA4q1bILeUCBkwQFjsih07IL+ScnAsfM7ge+RLBi6MkNnptBsM7O0e7+QVh2EzCNYJXR\n6TWThZrZCGAecASwyN2vb+J4ChyRiGnZhdyXK4EzGtgNzE8MHDPrA/wK+Cxwqrt/YGZDgEeB0wjG\nAT0HHOfuHs58cI27LzezRcC97l7eyPEUOCIplLjsQlllmWY+yBE5EThQ2wHhmQaB8xhwG/A0dYEz\nk6An3F3hPr8DZhH0lHve3YeG5ZOAs939qkaOpcARSaOGMx9o2YXsFHXgZMw1r5ldCGx19zUNGiJ7\nA39K2N4elh0AtiWUbwvLRSRmvTr0YsqwKUwZNqVe54PZS2cz6YlJmvkgT2VE4JhZe+D7wJi46yIi\n0SqwAkYcM4IRx4zgps/fVG/mgwm/mVC77EJJUQljisZw9JFHx11lSZGMCBygCCgEXk2Ys+0VMxtF\ncEXTL2HfPmHZdqBvI+WNmjVrVu3j4uJiiouLo6m5iLRKp8M7cdHgi7ho8EVAXeeDX6/9Nd969lvq\nfBCjiooKKioqUvb+cbbhFBK04XxqbjYzexsY4e4fmtlQ4BHgdIJbZkuo6zTwZ+A6YDnwLDDH3csa\neT+14YhkgZrOBzUTj6rzQbxyotOAmT0KFAPdgJ3Are4+N+H5t4CRDbpFfwPYT/1u0adSv1v09CaO\np8ARyULNLbugmQ9SLycCJ90UOCLZT8supJ8CJwkKHJHc09iyC2OLxlI6qFQzH0REgZMEBY5I7tOy\nC9FT4CRBgSOSX/Ye2BvMfJCw7II6H7SeAicJChyR/LZj9w4WVy6mvLKcJZVL1PmghRQ4SVDgiEiN\naq9m1buraq9+EjsfaNmF+hQ4SVDgiEhTGnY+0LILdRQ4SVDgiEhLadmFOgqcJChwRCQZ+T7zgQIn\nCQocEYnCzt07azsf5MOyCwqcJChwRCRqicsulFeWs/LdlTk384ECJwkKHBFJtaZmPsjmZRcUOElQ\n4IhIuuVC5wMFThIUOCISp31V+3h5y8tZN/OBAicJChwRySTNLbtwdv+zM6bzgQInCQocEclUzc18\nEHfnAwVOEhQ4IpItMmnZBQVOEhQ4IpKt4lx2QYGTBAWOiOSCppZdKB1USklRCf0794/0eAqcJChw\nRCQXpXrZBQVOEhQ4IpLrUrHsggInCQocEck3/9j7D17Y9EJt+08yyy4ocJKgwBGRfJfMzAcKnCQo\ncERE6tQsu1ATQG//7e1GZz5Q4CRBgSMi0rTGll0oLSplzgVzFDitpcAREWmZxGUXvv+F7ytwWkuB\nIyLSelHfUiuI6o1ERESao8AREZG0iCVwzOx+M9tpZq8llM02s/VmttrMnjCzTgnP3WRmb4bPj00o\nH2Fmr5nZX83snnR/DhERabm4rnDmAiUNyhYDJ7j7cOBN4CYAMxsKfBUYAowDfm51w2V/AXzD3Y8H\njjezhu8pDVRUVMRdhYyhc1FH56KOzkXqxBI47v5H4MMGZc+5e3W4+WegT/j4QmCBux9w900EYTTK\nzHoBHd19ebjffGB8yiuf5fTHVEfnoo7ORR2di9TJ1DacacCi8HFvYGvCc9vDst7AtoTybWGZiIhk\noIwLHDO7Gdjv7r+Ouy4iIhKd2MbhmFl/4Bl3Pzmh7DLgcuBcd98bls0E3N3vCrfLgFuBzcAL7j4k\nLJ8EnO3uVzVyLA3CERFJQpTjcFK3VNzBWfgTbJiVAv8KfKEmbEJPA4+Y2d0Et8wGAcvc3c3sIzMb\nBSwHpgJzGjtQlCdMRESSE0vgmNmjQDHQzcy2EFyxfB9oBywJO6H92d2vdvd1ZrYQWAfsB65OmDbg\n28A84AhgkbuXpfWDiIhIi+XF1DYiIhK/jOs0EDUzKzWzDeHg0O/FXZ90MbMvm9laM6sysxENnsur\ngbQaVFzHzG4zs1fNbJWZlYXDC2qey6tzUcPMZphZtZl1TSjLq3NhZrea2TYzeyX8KU14Lrpz4e45\n+0MQqBuB/kBbYDUwOO56pemzfxY4DngeGJFQPgRYRXA7tTA8PzVXun8BTgsfLwJK4v4cEZ2L84GC\n8PGdwB3h46F5eC46JDy+FvhFvp6L8PP0AcqAt4GuYVk+/o3cCnynkfJIz0WuX+GMAt50983uvh9Y\nAFwUc53Swt3fcPc3SeiYEbqIPBtI6xpUXMvddydsHgXUnJe8Oxehuwk6KyXKu7+RUGOdqyI9F7ke\nOA0HjWpwqAbS5v2gYjP7YdhZ52vAv4XFeXcuzOxCYKu7r2nwVN6di9A14W3nX5nZZ8KySM9FnN2i\n5RCZ2RKgZ2IR4MDN7v5MPLWKR0vORb4MKj7YuXD3W4BbwjbNa4FZ6a9lejRzLm4h6Bk7Jo56xaG5\n3wvg58Bt7u5m9kPgJ8A3o65DrgfOdqBfwnafsCwnuHsyfyzbgb4J2zXnpKnyrHCwcxEOKr4AODeh\nOC/PRYJHgWcJAievzoWZnUjQJvFqOBlwH+CVcFxfU/9u5OS5aMT/BWr+hzXS34tcv6W2HBhkZv3N\nrB0wiWAgab5JvDf7NDDJzNqZ2QDqBtLuAD4ys1HhH+BU4KkY6hq5hEHFF/qnBxXn27kYlLA5HtgQ\nPs6rc+Hua929l7sPdPcBBLeETnH39wjOxcX5ci4AEnsrAv8MrA0fR/p7kdNXOO5eZWbXECx9UADc\n7+7rY65WWpjZeOCnwNHAb81stbuP8/wcSPtTNKi4xp1mdjxBZ4HNwLcA8vRcJHLC/zHL03Mx28yG\nE/xebAKuhOjPhQZ+iohIWuT6LTUREckQChwREUkLBY6IiKSFAkdERNJCgSMiImmhwBERkbRQ4IiI\nSFoocEREJC0UOCJZyMwuNLNj4q6HSGsocESSYGZHmNl4MxtnZlc28lxFOMdUKo7dE7iMcCoWM2tr\nZi+amf6eJaPpF1QkOf8EPOXuvwPOaPDcNOAJT9G8Ue6+k2D12prt/cBzBJPTimQsBY5IK5lZH2Bj\nuHbIQGBLg10mkzBzrplNNbNXzWyVmT0YlvUP14ifa2ZvmNnDZnaemf0x3B4Z7nesmZWY2djwvzXh\n1vDq6anwuCIZK6dnixZJkZPdfVG4gNkpwI01T5hZW2CAu28Jt4cSLPT1OXf/0Mw6J7xPETAxnJF3\nBXCJu48OV6K8GZjg7u8A7yQe3Mx6AMcTrO3zcFi8FjgtFR9WJCq6whFJkrvfBcyl/q2so4G/JWyf\nCzzm7h+Gr0l87m13Xxc+fh34ffh4DdC/meO+5+6T3f3hhLJqYK+ZHZXs5xFJNQWOSOsl3hkoAj5I\n2N4DtG/h+yQuBledsF1NcncfDgc+SeJ1ImmhwBFpBTPrRrCcbo1S4L9rNsIrmIJwhVmA54GvmFnX\n8PVdEt+uuUO1sl5dgV3uXtWa14mkk9pwRFpnOLDezP6ZYE33WQ1uk0Gwwuxo4PmwfeZ24EUzOwCs\nIujFBsEqkzTyuLHtgzkHeLaVrxFJK634KdIKZnaBuy86yD6nANe7+6VpqhZm9gTwPXffmK5jirSW\nbqmJtM5Bb1m5+yrghVQN/Gwo7Bn3PwobyXS6whERkbTQFY6IiKSFAkdERNJCgSMiImmhwBERkbRQ\n4IiISFoocEREJC0UOCIikhYKHBERSQsFjoiIpMX/B+CKhCQBi6wtAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10dc01da0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(beta_plot, eigval_a_plot, label='asymmetric')\n",
"ax.plot(beta_plot, eigval_s_plot, label='symmetric')\n",
"ax.legend(fancybox=True, loc='best')\n",
"ax.set_xlabel(r'$\\beta$ (cm$^{-1}$)')\n",
"ax.set_ylabel(r'resulting frequency (cm$^{-1}$)')\n",
"ax.set_xlim(ax.get_xlim()[::-1])\n",
"fig.savefig('mpl.pdf', bbox_inches='tight')"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAABIAAAATBAMAAABvvEDBAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIokQdkQymVRmzbur\n3e+SS/cOAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAfUlEQVQIHWNgYGBgBGIIYIIxGBAsE7hYGJyl\nAGWxBRY2QpgVjBO4BBiyGBi4J3AKMB1IP8vAMBdoCPsFhl4Ghu1AFtcBEOslgySDfAOItYRhIlAc\nxIplUOAuALPYKreUAY0DijEwTgCZC2JxCoBYfUAMdpP4/SNAMZAQEAAAim4WL+0cFKgAAAAASUVO\nRK5CYII=\n",
"text/latex": [
"$$\\phi_{1}$$"
],
"text/plain": [
"φ₁"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phi1, phi2 = symbols('phi_1, phi_2')\n",
"phi1"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAOgAAAAZBAMAAADJSX/3AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAu90iEM0yiXZEmVRm\nq++H7A6pAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACvElEQVRIDc2WP2gUQRTG3+3e7XK5NYVYyGEx\nbBS1ULeOkRzY2CmRBUEC16QwjQtiExC2tgoRtD2j2KgQtRECuprOahEESxs7QUEbETnfzLzZ2dvd\ndxcrfZDMm9nvfb95++8WDoYRtMIF2HdshUmTVvrMDFnbDY+DQKUzU10WuHl5ZnJhkqmjqj3XBGW2\n7JLdNKjRVMnGkocm1RI99/cBNZqqg7H8v6DttLTPnp1coOVpp9dolPRl4dNOKWU77Q5JIYcS9C4t\nT4MajZJaaGHJQieuSwk6cV1oA3YQKjUaNbHQwpKDXls8G1mvAurdefBCL/OdWs0k1Foy0LmdJbjc\nAH3sZh29GR5qNRPQkiUDPQJ7cAu8S8Q1nbayIPJz//5DABaqNbBmtmxOr7bc6A+otng5tFdCjGPD\nQEAC295TISv9OF75HMc5pm/BhzlxAz5MgWrN+ntZCjfj+HQcr2KmLf3oADbS3GlrhNBV6AlZiGE6\n7SO0ky/A9YjvVGvgua4EoE61pRN53znofOZm7q869AsswfLOSVje5KFaU4NqS2fT/clBO2kQ4dZq\nnV7BS93HFl4PeChpqp2SJfR+cNCu8KV5DboFCe4F4AQVmlNYjAKANFUoWYKTclBYf3N72AD1nhx+\nhAB/ZKAXzf2NqzIEAGmqULKEDVQ130iA1w6j1im4mVx/h3/6kXGOyrkNganW1KDaMhhJhXwSikcG\nc4o9ORbQdkTLgUy6SXdorim9n+iw+l1WGnsjnTfHlOU9wBczB02ktPfJFJhRvT7Xdp8NDFTpzFF1\nevGpkvHKLlImpd7V3UM8NEOFu/07xaEcgZx8HY+pENopzkshMFeaxTMfS8sqzfD//Hj8jYdWK2pz\nfU0VoHRMlHI+5U4vX0FHNLQqE9WFxvm/hIb4ver83Xdv3tSE9JkZW2GO372n/gBqoLuB4N3QrQAA\nAABJRU5ErkJggg==\n",
"text/latex": [
"$$\\left [ \\left[\\begin{matrix}- \\phi_{1} + \\phi_{2}\\end{matrix}\\right], \\quad \\left[\\begin{matrix}\\phi_{1} + \\phi_{2}\\end{matrix}\\right]\\right ]$$"
],
"text/plain": [
"[[-φ₁ + φ₂], [φ₁ + φ₂]]"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"normal_modes = [nm.T * Matrix([phi1, phi2]) for nm in normal_mode_eigenvects]\n",
"normal_modes"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2250.00000000000\n",
"1250.00000000000\n"
]
}
],
"source": [
"print(eigval_a_plot[-1])\n",
"print(eigval_s_plot[-1])"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.1"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment