Skip to content

Instantly share code, notes, and snippets.

@goerz
Created April 13, 2019 17:38
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 goerz/ae7db4b2f1bc0219c6cd5f1787d54bb4 to your computer and use it in GitHub Desktop.
Save goerz/ae7db4b2f1bc0219c6cd5f1787d54bb4 to your computer and use it in GitHub Desktop.
Frequencies vs Angular Frequencies in Quantum Dynamics
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Frequencies vs Angular Frequencies in Quantum Dynamics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Author: Michael Goerz <<http://michaelgoerz.net>>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook checks that for an energy separation of $\\omega$ in the Hamiltonian ($\\hbar = 1$, as assumed by `qutip.mesolve` and any other propagator for quantum dynamics), a resonant driving field should oscillate with $\\cos(\\omega t)$, not $\\cos(2 \\pi \\omega t)$.\n",
"\n",
"Furthermore, when plotting the spectrum of a control, it is most intuitive to use angular frequency on the x-axis, so that matching peaks in the spectrum to transitions in the Hamiltonian is direct."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import qutip\n",
"import matplotlib.pylab as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"w = 10.0\n",
"mu = 1.0\n",
"H0 = qutip.Qobj([[0, 0], [0, w]])\n",
"H1 = qutip.Qobj([[0, mu], [mu, 0]])\n",
"def eps(t, args):\n",
" return np.cos(w*t) # NOT 2π ω t !!!\n",
"H = [H0, [H1, eps]]"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"tlist = np.linspace(0, 10, 2000)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"ket0 = qutip.ket('0')\n",
"ket1 = qutip.ket('1')\n",
"\n",
"proj0 = qutip.ket2dm(ket0)\n",
"proj1 = qutip.ket2dm(ket1)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"res = qutip.mesolve(H, ket1, tlist, e_ops=[proj0, proj1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For a resonant pulse, we expect full Rabi-cycling (with small osciallations on top of that, if $\\omega$ is small enough for the rotating wave approximation not to hold; as $\\omega \\rightarrow \\infty$, these small oscillations will average out."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def plot_population(result):\n",
" fig, ax = plt.subplots()\n",
" ax.plot(result.times, result.expect[0], label='0')\n",
" ax.plot(result.times, result.expect[1], label='1')\n",
" ax.legend()\n",
" ax.set_xlabel('time')\n",
" ax.set_ylabel('population')\n",
" plt.show(fig)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXd4XNW1t9896mUkWb3LsqptucvGxoBNM2BaAoQWklASchOS8N0QLqTc3JvkJiHJTSFAbmgpQAKhJTRT3TFuMrblpi7L6tXqHrXZ3x9bso2xrZHmzJzRzH6fZ57jOTqzz0/yObPOXmvttYSUEo1Go9FoACxmC9BoNBqN56CNgkaj0WiOo42CRqPRaI6jjYJGo9FojqONgkaj0WiOo42CRqPRaI6jjYJGo9FojqONgkaj0WiOo42CRqPRaI7jb7aAiRIbGyunT59utgyNRqOZUuzatatNShk33nFTzihMnz6doqIis2VoNBrNlEIIUePIcdp9pNFoNJrjaKOg0Wg0muNoo6DRaDSa42ijoNFoNJrjaKOg0Wg0muO4zCgIIf4khGgRQuw/w8+FEOL3QogKIUSxEGKhq7RoNBqNxjFcOVP4C3D5WX5+BZAz+rob+D8XatFoNBqNA7hsnYKUcpMQYvpZDrkWeEaqfqDbhBBRQogkKWWjSwQd2Q7VmyBjGWQsByFcchrNp2npsVF0+CgWIViUMY04a5DZknwHux0q3oemfRA9A/JWQ0Cw2ap8BtvQCFur2umxDTMrKYKsuDCEh3/3mLl4LQWoPel93ei+TxkFIcTdqNkE6enpkzvbka2w/n/Uv5Pmw4Xfg5xV2jgYTHvvAM3dA/QPDlPV2seb+xr5sLwV+2grcCFgYfo0LpudwHnZcSRHBRMZEuDxN8qUQ0ooeQvW/wxaDpzYHxYP534TFt8FgWHm6fNCbEMjNHfb6LEN09xt4/2DzazZ10i3bfj4MZmxYVw2O5Hl2TGEBvoRHRZERnQoFovnXP9CPai7aHA1U3hTSllwmp+9CTwkpfxw9P1a4AEp5VmXKxcWFspJr2ge6IX9r8DmX0NnDViTISgcpB1ic+HC70Pip6RqHKCl28a3X9zLhxVtn9ifEhXCZxYkc9nsRIbtkg/L23j3QBMHGrqPHxPkb+GSWQk8dN0crMEB7pbuHdR8BJt/A90NyvL2t0NPI8Rkw4oHIX811O6ALb+Dqg0QEAZR6cowBIXDgi/AnBvM/i2mJMMjdh56u4RnttUwOGw/vj800I/LZydyzfxkEiODKTp8lHcPNLG1sp1h+4nv3dyEcH5/ywLyEyNcqlMIsUtKWTjucSYahceBDVLK50fflwIrx3MfOWUUxhgZgr0vKHfSyKC6iao3wZAN7npPG4YJMjRiZ/XDm6nvPMbXVmSRHR9OeLA/cdYg8hKsp50F1Hb0s6e2k+ZuG9VtffxjZy3LsmL46x1LPOqpaUpQ8QH87UawJkLyArUvIETNhGdfB36nOARqd8C+l6G7Hob64WgNdFTCVb+Fwjvdr3+K88PX9vPM1hpuWJTKOZnRRIQEEBkSwPy0KIID/D51fGf/IAcbuxkctlN39Bi/X1sOwBvfPI+ECNe59hw1CkgpXfYCpgP7z/CzK4G3AQEsBXY4MuaiRYukS+iql/JXuVI+tkzK4SHXnMNLeeajapnxwJvy3f2Nkx7juW2HZcYDb8pntx42UJkPMNAr5S9mqOv2WNfkxhgekvLZ66T8SbyUbRXG6vNyypq6ZeaDb8of/mvfpMcoaeyWeT9YI7/+t10GKvs0QJF04DvWlSmpzwNbgTwhRJ0Q4i4hxL8JIf5t9JA1QBVQATwJfN1VWhwiIhmu+IXyv+5+1lQpU4m+gWEeXlvOksxoLp2VMOlxbl2SzjmZ0fzug3KODY4YqNDL2fkU9LfB1b+D4Em6H/z84ZpHQfipGITGYX71bilhgf7ce0nupMfIS7TybyuyeKu4kR3VHQaqmxwuMwpSyluklElSygApZaqU8mkp5R+llH8c/bmUUt4jpcySUs6R48QS3MKsayHtHNj4S+VK0ozLU5uraesd5MEr8p0KFgsh+M5lebT1DvDcNoeKOWoG+2DLw5B1MaQtcW6siCRY8mUVc2stM0afl7Or5ijvHWzm7gtmEB0W6NRYX70gi4SIIH7xTsmYJ8U09IrmkxECLvpP6GmAoqfNVuPxtPUO8MSmSi6fncjC9GlOj7d4ejTn58Tyx42V9A0Mj/8BX2fnUyqgvPJBY8Y791sqFrHpl8aM58VIKfnF2yXEhgdx1/mZTo8XEujHvRfnsqvmKGsPtRigcPJoo3AqmedD5gqVyTHQa7Yaj+bRdRXYhu3cf3meYWP++6W5tPcN8pePDhs2plcy0Ds6S7jI+VnCGGGxsHh0ttBWYcyYXsr60hZ2HO7g3ouzCQ00JrP/c4WpZMaG8at3Sxmxmzdb0EbhdFz8Q+Wn3f5Hs5V4LBUtvTy3rYYbC9PIigs3bNyF6dO4OD+exzdW0tU/ZNi4Xsf2/xudJXzX2HHP/Sb4BcHm/zV2XC9iaMTOz9eUkBETys1LJrlu6jQE+Fm4b1Uupc09vLan3rBxJ4o2CqcjtRByr4CPfg/HjpqtxuOwDY3wwCvFhAb68Z1Vkw+wnYn7VuXRP6jOYeYTk8fS3QBbfq9WJxs1SxgjPF4tbCv+B9RsNXZsL+EP6yspb+nlB1fOIsDP2K/Q1QVJFKRE8It3Sqg72m/o2I7i0nUKrsCQdQqO0LQPHr8A0pbCsnsgPAHi8iaf4THFqe3o5+87jlDW1MOBhm6aum08eusCrpqb7JLzPbmpip+uOcSMuDDyEqykRYdyx/LpJEWGuOR8Hs9gP3RUQV8LfPDf0FYOX90MsdnGn2ugB/54Hti64fxvgzUJ4mdCwmzjzzUFGBy28/reBvbVdVLW3MvWqnY+Mz+Z39403yUr8ffXd3HLE9sYtktmJllJiAjm5iXprMgdt73yWfGIxWuuwG1GAaD4RXjz2zDYo94HR8Gt/4D0pe45v4dQ2drLZx7bwrHBEbLjw8mICeXGwjQunjn5FFRHeH1vAy/urKWp20ZNex8xYUG89a3ziAn3sdpJFR/Ay3eBrVO9DwyHG/4EuZe57pztlfDS7dBUfGLfyu/Bygdcd04PZGjEzm1PbWd7dQfWIH8yYkO5ICeO/3dJLoH+rnO0VLb28uct1VS09FLZ2kdb7wBPf6mQi/Inf89po2AUAz3QVgY9TfDeD9Rq6Hu2+1TdmJuf2EpJUw+v3bOcjBhzfu/99V189g9buG5BKr+4Ya4pGkyhpxkeWQjTpqun9pBoSJ4PIc5ne42LlNDXCv0dKsaw7yW4633jXVYezBObKvnZmhIeum4ONy1OM6VG17HBET7z2Bb6BodZd9/KSRsjR42CjimMR5AVUhZB/pVqgU9XLez5u9mq3EZxXSfbqjr4xoXZphkEgIKUSG5dks4/d9fT2jNgmg638+FvYdgGNz4DBddD1oXuMQigUrTD4yE+H676HYTFwSbfCUAPDtt5+sNqzs+J5eYl6aYVbQwJ9OPB1fnUHT3GW/saXH4+bRQmwvTlkLwQiv5kthK38fKuOkIC/LhpcZrZUvjCsukMjij/rk8wPADFL8DMayAmy1wtQeGw6A4ofw+6zMuMcScbSlto7h7gzuXOr0NwlpW5cXz5vEzyElwf09RGYaLMuxlaDqpAnw+wsayVc7NiPKJ6aXZ8OPmJVt7Z75qWGx5HxQcq+23+581Woph7IzBaktsH2FDWSniQP8uzY82WghCCH1w1i1nJ2ih4HvlXqu2h183V4QYOt/VR097Pijznsh6M5PKCRIpqjtLW6wMupPL3IdAKM1aYrUQRmwNx+T5x7Usp2ViqHohcGVD2RHzrtzWCyFRVnrjsXbOVuJwNpWq5vbOpcEZyYV48UsKWU/o2eB1SQuVayLwA/MyfpR0nb7VqWGXrMluJS6ls7aO+85hHPRC5C20UJkPOZVC3U2VleDEby1rJjA0zNcB8KgUpkUSGBPBhuZcbhY4q6DyiAsueRM4qsA9D5TqzlbiUjWWtAFyQo42CxhFyVqlubRVrzVbiMsZ6y3rSLAHAzyI4NyuGLRVtpleTdCljX7pZF5mr41RSF6v1OuXvm63EpWwsayUrLoy06FCzpbgdbRQmQ/ICCI1VmRheys7DHdiG7B5nFACWZ8fS0KU6tnktFWvV2gSzs45Oxc8fsi9W177dPv7xUxDb0Ajbq9pZkRtvthRT0EZhMlgskHOpyg6xe2dDmI2lrQT6W1g6I8ZsKZ/i/ByVDXJqP2ivYWQIDm+GGR7mOhoj5zK1qK1xt9lKXMK2qnYGhu0+GU8AbRQmT86lcKwD6neZrcQlbChr5ZzMaEICP91j1mwyYsJIiw5hs7fGFep2wmCv57mOxsi+BBBe60LaWNZKkL+FczKjzZZiCtooTJasi1T7Qi90IdW091HR0svKPM+dPq/IjWNLRRu2IS+cqVWuA2FRmUeeSFiMqiTshRl4UkrWlbSwdEYMwQGe90DkDrRRmCwh01TrTi80Ch+Mdn661MUF75zhkpkJ9A+OsLWy3WwpxlO5HlIKISTKbCVnJucyaPgYes3tEmY0la291LT3O9VvfKqjjYIz5FwKjXtVsTwv4v2DTeQmhJMe47mZF8uyYggL9OPdA971t6e7Qbkkcy41W8nZGdNX8YG5Ogzm/YPKyF0803Nnya5GGwVnGCtd7EW+1SPt/Wyr6uDKOa7pk2AUQf5+XF6QxBt7G+ixeVGHtv2vAhJmX2e2krOTNA/CE73KhSSl5NWP65iXFuW7fTvQRsE54mdBRKrqUuUlOfN/2FCBv0V4RAG88fjSuRn0DY7wpw8Pmy3FGOwjsOvPKuXZFc1zjEQI9VBU/p7XzJQ3lLZS3tLL5w1ssTkVMabjtK8ihOpp+84D8OFvYMEXIXzqpbGVNfewrqSFipZeXt5Vx1fOzyQxMthsWeMyNzWKq+cl8/DaMspbepgeE8bV85LJS7SaLW1iSAm9zar6bnsFfO6vZityjOX3qjLyL90BK+5XD0nWRLNVTYiBYRWXqu3o57H1lcyIC+MzC1LMlmUqusmOs4wMwYtfhNI16n14guq7kLvKXF0Osr2qndue3s7QiCQ8yJ/VcxL5yWcKCPKfGpkXfQPD/HTNIdaXtNDcbSPQ38KrX1vulmqShnDoDXjrO9A7+rQ961plFEyq3T9hil+E178Fw8fU+/m3wbWPTgn9fQPD3Pj4Vg40dAOQEhXCU18qZGbSFLl2JojuvOZO7Hao3Q4Nu2HXX9RT37d2Q6hn5znb7ZJLf7uREbvkxa8uIz7C82cHZ6Ol28YVD28mP8nK3748BVqmHj0Mjy1V1UcXfEGtXp6xEixTwyAfx9YFDXvgwKvq+r/pOZh5tdmqxuXnaw7x5OYq/vdz81iWFUO8NRg/i+cbs8miO6+5E4sFMpbBsq+r3rm2TtjzN7NVjcvmijYqW/v490tzp7xBAIiPCObO8zLZUtFOTfsUKIHx0aOqhtYtL8A5d6vyEVPNIAAER6ry3lf+RpXm2P642YrGpW9gmL9vP8KVc5O5bmEqSZEhXm0QJoI2CkaTWKDWL0yBlp3vHmgiLNCPywumlh/4bHx21B/8ZrGHN+IZHoR9L8KsayDSS3zYFj9YcJsq0dFVZ7aas/JhRRs9A8PcMgUSKtyNNgquIP9K1Z2t23PbRkopWXeohfNz4qZM/MARkqNCyE+08lGlh5fAqNup3C6zrjVbibHkjTah8vAKwusOtWAN8mexj5ayOBvaKLiCsUJmVRtMlXE2Klv7aOq2sdILi34tz46l6PBRzy6BUbVelUnx1FIWkyV+JliT1e/nwWwub+X83FgC/PRX4Knov4grSCiA0Bio3my2kjOy+8hRABZlTDNZifEsmxHDwLCdvbWdZks5M1UbIWWh8sd7E0Ko+FrtTrOVnJHmbhsNXTYWZehZwulwqVEQQlwuhCgVQlQIIR48zc/ThRDrhRC7hRDFQojVrtTjNiwWVbum4WOzlZyRPbWdWIP8yYoLN1uK4cxPVzWD9tV7aMvIkWFoKobUJWYrcQ0phdBdB92eGdfZM/qwMD/Ng2tLmYjLjIIQwg94DLgCmAXcIoSYdcphPwBelFIuAG4G/uAqPW4nZSG0lsJAj9lKTsvuI53MS4vC4oUZF7HhQaREhbC3zkONQmsJDNvUymVvJHWx2tZ7WOr4KHtrO/G3CGZPlbUsbsaVM4UlQIWUskpKOQi8AJwaVZPA2P9MJOC5kdmJkrwQkKpgnodxbHCE0uYeFqR775PSnJRI9tV5qPuocY/aJs83V4erSJwDlgAVTPdA9tR2MjMpwmdLY4+HK41CClB70vu60X0n89/AbUKIOmAN8E0X6nEvKQvVtt7zXEgHGroYsUvmpXqxUUiN5HB7P139Hlgsr2E3BFoh2sNabRpFQLBKzW7wvM5sdrukuK6LeWleFssxELMDzbcAf5FSpgKrgWeFEJ/SJIS4WwhRJIQoam1tdbvISREWC5HpHhlX2D/qa5+T6r03xtzR380j4woNe1SVUYvZt58LSV6ofk8P6+Nc3d5H78Awc734gchZXHlV1gMnrwxJHd13MncBLwJIKbcCwUDsqQNJKZ+QUhZKKQvj4qZQCmXyfHVjeBj7G7qJDQ8i3hpkthSXMTdF3fTF9R7mQhoZgqZ93us6GiN5AQx0Q0el2Uo+wfEHohTvfSByFlcahZ1AjhAiUwgRiAokv37KMUeAiwGEEDNRRmGKTAUcIHk+HK2GY0fNVvIJ9td3UZASgZgCRcsmS2RoABkxoezztGBzawmMDECSlxuFMfeph7mQDjR0E+hvITve+7LujMJlRkFKOQx8A3gXOITKMjoghPixEOKa0cPuA74ihNgLPA/cLqdahb6zMXbje1Cw2TY0QnlLLwXJ3v+kNCclkmJPMwoNXh5kHiM2DwJCPS6mtr++i5mJVr1o7Sy4tJ+ClHINKoB88r4fnvTvg8ByV2owlbGUw4Y9qvqlB1Da1MOIXVKQ4v3peHNTI3mzuJH23gFiwj3EVebtQeYx/Pwhca5HzRSklOyv7+KqeZ7dVdBstLl0JaHREJV+IgXRA9jfoJ6cZ/vETGEsruBBs4XGPWqW4M1B5jFSFqpZ8siw2UoAqO04Rrdt2Cdmyc7gA1emySQv8Kinpf313USGBJA6zft70Kq4CZ4TVxgZgqb93u86GiN5gWq+01ZqthLgxAORL8ySnUEbBVeTNF81U/GQYPPuI0eZmxrp1UHmMazBAcyIDfOcuELLId8IMo+R7FlrdfbUdhLoZyE3YYq1a3Uz2ii4mrGnQg+YLXT1D1Ha3MOS6b5TCGxuahR76zrxiPyFsRW+3lre4lSiZ0BQhEdc+wA7qjuYmxqpVzKPgzYKriZlkSqRfHiL2UooqulASnyqhvzi6dG09gxQ3tJrthTVfMaarL4sfQGLRT0U1e8yWwn9g8Psr+9iiQ9d+5NFGwVXExypCoRVmt90ZGNZK8EBFp+qDjnWL2JDaYu5Qux2VUo984Ip0dTeMNLPVcHmXnP//turOxi2S20UHEAbBXeQdZFKS+1pNk3CiF3y7oEmVubG+9T0eawT25p9TeYKadgN/W2ql7EvMfNqQELpmnEPdSXv7GsiPMifZVkxpuqYCmij4A4KrgMk7H7GNAnv7G+iuXuAa+f7Xo72TYvT2FPbyfaqdvNE7P07+AWpVq2+RMJsiMmGoj+DSXGdrmNDvLWvkVWzEryq9ayrcOniNc0osTmQfSl8+LBa6RmbC3F5LncjDA7beeXjOorrOnmruJG8BCurZie69JyeyI2FaTy1uZqvPFPEeTmxFKREctd5me75gujvgJotsPs5KLje+zqtjYcQcP598K+vwevfgLSlauYceWrBZOM53NbHtqp23ihuoH9wmC+f7yOxHCcRHpGVMQEKCwtlUZFnNu84K51H4E+XQ/doTcC8K+HGZ9TKTxcgpeRrz33MOweamBYawMykCH722TlMjw1zyfk8ncrWXn71TiklTd0cbu/nliVp/Py6ua474fAAvHo3HPyXem9Nhq+shQjfm6lht8Pb98POpwEJQZFw17uqn7OLeP9gM197bhfDdkmgv4UfXDmTLy6b7rLzTQWEELuklIXjHqeNghsZ7FdBt8p1sOmXcPXDsOh2l5xqe1U7Nz2xje+syuWeC7N9Yl2Co/zkzYP8aUs1a7+9ghmuake68Zew/qew/F71ZJyyCIJ8PD9+yAZtZfDMtSot9wuvuuQ0tqERVvxqPdFhQTx26wLSo0Px17WOHDYK+i/lTgJDVVPzC7+nFjDtfMplp3phZy3WYH/uOm+GNgin8NULlBvh9b0uavQ3Mgw7noTcy+HSH6u6V75uEEA130maC0u/prLxuupccpr1JS00dw/w4BX5zIgL1wZhgui/lhkIAXNvUnX12yoMH95ul6wvbeGy2YmEBOrA2qnERwRTmDGNdw+4KBusegP0tcDCL7pm/KlOwfVqe/DUSvrG8Oa+RmLDgzgv+1OtWTQOoI2CWeRdobZV6w0f+lBTN539QyzP1ul3Z2JFbhyHGrvpOuaCdp1VG8EvEGZcaPzY3kBMlspIqt5k+NBSSrZVtnNBTix+Fj1DngzaKJjFtOkq+FjzkeFDb61UqZfLZugnpTOxMH0aoGpBGU7NFhVDCAw1fmxvIX0ZHNlqeLvOqrY+2vsGfWrVvtFoo2AWQqj4Qu12w4feW9dFSlQIiZHBho/tLcxLi8Ii4OMjBrfrHOxTCxUzvLdNiCGkLwNbJ7SXGzps0eEOQJU30UwObRTMJGm+SlHtM3ZR1aHGbmYl6/LAZyMsyJ/cBCvFdQYbheaDIEdOtKPUnJ6k0XTgpn2GDnuwoZvwIH+y4nwz9doItFEwk8Q5atts3I1hGxqhqrWXmUnaKIzHzKQISpt6jB20eb/aJsw2dlxvIzYPLAGGG4XS5h5yE8J1xp0TaKNgJmNGoWm/YUOWNvVglzArSadAjkdeopXGLpuxwebmA6rdZmS6cWN6I/6BEJd/wogagJSS0qYe8hL1te8M2iiYSVgshCca+rR0sLEbgFlJPlZOYRLkjTZbKWs2cLbQfAASZvlGu01nSZxj6LXf2jvA0f4h3UTHSfSVazaJcwx9WjrUqHyqvtBu01nGnihLjHIhSTlqFAqMGc/bSSyA3mbobTVkuLIm1TMjTxsFp9BGwWwSZqml/yPGuDDKmnvISQjHonO0xyUpMhhrsD9lRhmFrjoY6NLxBEeJn6W2LQcNGa50dMan3UfOoY2C2cTPgpFBaK80ZLjy5l5y4/VN4QhCCPISrMYFm8e+3LRRcAyDjUJZUw+x4YHEhAcZMp6voo2C2Ry/MQ44PVR77wDtfYPkJLioyJsXkpdopbS5x5gezi2H1DYu3/mxfIHweAiNUS43Ayhp7tHxBAPQRsFs4vJUD+exLxQnGOtDnKNvDIfJS7TSdWyI5u4B5wdrLQFrEoT4TrtTpxBCPRQZcO3b7ZJybRQMQRsFs/EPUnVgmp2fQpeP+lRz9UzBYca+REqNyEBqOaRnCRMlYbb6uzlZ7qK+8xj9gyM6nmAA2ih4AvEzDXEflbf0Yg3yJzFCl7dwlDGj4HSw2W5XCQMubBzjlcTPhKE+6KxxapixDDJtFJxHGwVPIGE2HD0MA71ODVPW3EO2Xs05IaLDAomzBjmfltpZA0P9yh2ocZz40aC8ky6k0ia1Pke7j5xHGwVPYCzY3Frq1DAVLb3kxGvX0UTJT7Q6v4CttURt4/RMYULEj7rbnJwplzb3kjothPAg3XbeWbRR8ATGXA5O3BgdfYO09Q7qJ6VJkJtgpbylhxG7ExlIx42CnilMiCArRGU4HVMrberWi9YMwmGjIITwE0IkCyHSx16uFOZTTMuEgFCnboz99V0AuhDeJMhLsGIbsnOko3/yg7QcUv0xdObRxEmY7dRahcFhO1WtfTqeYBAOGQUhxDeBZuB94K3R15sOfO5yIUSpEKJCCPHgGY65UQhxUAhxQAjx9wlo9x4sFpW14sSNsW/UKBQk65pHEyV39MtkzC89KRr2nCgHrZkY8TOhvQKGJ5cWXNnay7BdaqNgEI7OFO4F8qSUs6WUc0ZfZ70DhBB+wGPAFcAs4BYhxKxTjskBvgssl1LOBv7fhH8DbyF+lnNGoa6L6TGhRIYGGCjKN8hLsOJvEeyp7ZrcAAM9KvMoWfdQmBTxs8A+DG2Ta7ize7RR0rxUPUszAkeNQi0w0TtmCVAhpaySUg4CLwDXnnLMV4DHpJRHAaSULRM8h/eQMBv6WqG7ccIflVKyt66TghQ9S5gMIYF+zE2NZOdo164J07gXkJC8wFBdPsNYCfmGjyf18V01R4kJCyQjRrc/NQJHjUIVsEEI8V0hxLfHXuN8JgVlTMaoG913MrlArhBiixBimxDi8tMNJIS4WwhRJIQoam01pqKix5FxrtpWb5zwR6va+mjssrF0RozBonyHJZkxFNd1cmxwZOIfPrJNbXW3tckRmwvhCVC1YVIf//jIURZmTNOp2AbhqFE4goonBALWk17O4g/kACuBW4AnhRCfmgNKKZ+QUhZKKQvj4uIMOK0HkjhX1YGpXD/hj24qU4byghwv/du4gXOzYhgakWwqn8RDR+U69f8XFmu8MF9ACJixUhmFCa5srm7ro7qtTz8QGYhDSb1Syh8BCCHCR987ssqqHkg76X3q6L6TqQO2SymHgGohRBnKSOx0RJdXYbFA1sVQ/i4MHYMAx/ohSCn55+568hKspOvp86Q5NyuG2PBAXiqq47LZiY5/sLcVarfDud9ynThfIPcyKP4HVK2D7Esc/thbxQ0AXF4wgf8zzVlxNPuoQAixGzgAHBBC7BJCjFcfeCeQI4TIFEIEAjcDr59yzL9QswSEELEod1LVBPR7F4u+BMeOws6nxj10eMTOrpoOfvVuKcV1XXxhWYYbBHov/n4WPn9OBh8cauY375XyUUUbtiEHXEk7nlBB0nk3u16kN5N/tXIhrfupKiM/PHjWwweGR9he1c5fPjrM8uwYUqJ0UymjcHT53xPAt6WU6wGEECuBJ4Fzz/QrnsIhAAAgAElEQVQBKeWwEOIbwLuAH/AnKeUBIcSPgSIp5eujP1slhDgIjAD3SynbJ/3bTHUylkPOKnjvB7DlYUhdDNc88im3RGf/IDc/se14aYYVuXHcvDjtdCNqJsDXVmZxoKGb36+r4PfrKshPtPLq188lNPCU26SvHdb+COqK1ILDguv1ojVn8Q+EK34JL98BjyyEgDD4zB9g9mc+deiTm6r4zftlHBsaITTQj++tdt8q8qGhIerq6rDZbG4750QJDg4mNTWVgIDJZSIKR+rICyH2SinnjbfPHRQWFsqioiJ3n9Z9DB2DXX9VvWv3vQTZF8Mtz3/ikP9+/QDPbqvhF9fP5ZzMaFKnheggm4G09NhYX9LCA6/s4/7L8rjnwuwTP5QS/rwa6ouUHzxxLpz/bQgMM0uud9FaCnU71Wy5tQz+X/EnHoqKDndwwx+3clF+PDcWplI4PZpYNzbVqa6uxmq1EhMT45H3nJSS9vZ2enp6yMzM/MTPhBC7pJSF443h6EyhSgjxn8Czo+9vw5fdPK4kIASW/pv697TpsP5/oKXkeI0Y29AIr+yq49p5ydywKNU8nV5MvDWYmxan8/reBl7YeYSvr8w68QVweDMc+Qiu+i0U3mmuUG8kLk+9UpfAY4vh42eU0R3lT1uqiQ4L5LFbFxIS6Od2eTabjenTp3ukQQDVTTAmJgZnsjQdzT66E4gDXh19xY3u07iShV9U20NvHN+1qayVnoFhPrPg1OxejdFcPTeZ2o5jx5sXAbD3BQiKgHm3mCfMF4jLhZRFULrm+K6+gWHWHmrhqrlJphiEMTzVIIzhrD6HjIKU8qiU8ltSyoWjr3vHFpxpXIg1Qd0YZe8c37WtqoMgf4tOwXMDy7OV22JLRZvaISVUrFVxHwezwzROkHu5itv0qTBjUc1RBobtXDorwWRh5vLOO++Ql5dHdnY2Dz30kOHjn9UoCCF+N7p9Qwjx+qkvw9VoPk3mCmjco2INwK6aDuanRRHorwvcupq06FDSokPYUT260rmjCnqbTiw01LiWjOWAhPpdAOw63IFFwIL0aebqMpGRkRHuuece3n77bQ4ePMjzzz/PwYPOd208mfFiCmMxhP819Kwax0lZpFIem/bRn7CQ/Q3dfG1FltmqfIa5KVEU16vaOhzZqrYZy80T5EskzQNhUUYhdxU7Dx9lVnKET/dM2LFjB9nZ2cyYMQOAm2++mddee41Zs2aN80nHOetfV0q5a/Sf86WUD5/8MyHEvcDEazJoJkbKIrWt30WpzGHELpmbqmscuYvZKRG8ta+Rrv4hIus/hqBInX7qLoLCVdOi+l1IKdlf38W1C5LNVnWcH71xgIMNTlTWPQ2zkiP4r6vPvASsvr6etLQT6eepqals377dUA2O+iC+dJp9txuoQ3MmIpJUnf6G3ccDnrqRjvuYPVqK/EBDl+qZED9TlWXQuIfkBdC4l6ZuGz0Dw7qRjhs460xBCHELcCuQeUoMwQpMsqSkZsIkzIKWQ1QE9xLkbyEtWpezcBezk1XTooMNXZzbchBmf9ZkRT5GfD7seY7qI6q2Zo4HGYWzPdG7ipSUFGprT9QZraurIyXF2EzE8ZxzHwGNQCzw65P29wDFhirRnJm4fDj8IeWBnWTFheNn0U+q7iImLJCo0ABaGo+ArfNEP22Nexjted1eXQyE+HwP8sWLF1NeXk51dTUpKSm88MIL/P3vxvYmGy+mUAPUAMsMPatmYsTlwbCN3qZqcmboxvDuRAhBdlw49qbRWv/x+u/vVkbjN0NNh4gOW0qMG1cveyL+/v48+uijXHbZZYyMjHDnnXcye7axMxaHwvhCiKXAI8BMVPlsP6BPSqkbAruDOLWa2dpbSU78YpPF+B7Z8eGE7C9Tb7RRcC+RqRAYTtDRcnLiHa+e6s2sXr2a1atXu2x8RwPNj6L6HZQDIcCXUa02Ne4gNheAXFFHdrzn+FR9hay4cNKGDmMPjdM9E9yNEMjYXKKPVZGT4NuuI3fh8AooKWUF4CelHJFS/hk4bZc0jQsIiaI/OIEcSx25+sZwO9nx4eRa6uiNzDFbik9ii8phhqwlRz8QuQVHV4H0j/ZE2COE+CUq+KyX1LqR5sAMcvsbSNeZR24nOy6MKNFAY+BitL/U/TQFZpApOsmPmkSrVM2EcfSL/QuoOMI3gD5UR7XrXSVK82kqZCrZlgb8deKR20n268QqjlEpdRFCMyhHVQPOtTSYrMQ3cLQgXo2U8piUsltK+SMp5bdH3UkaN7F7IIEQbNBVO/7BGkPx6ygHoNjm24XYzGKPTbXajOqtNFmJbzDe4rV9wBm78Egp5xquSPMp+geH2dEbr/K+Wkthmm696VZaVebRR126Mq0ZFB0NwyaCCG4rNVuKTzBeTOEqt6jQnJWq1j7K7aOui9YSyF1lriBfo62UAb8w9nYF02Mbwho8uTaHmokjpaS0pZ+24OmktpaYLccjuPPOO3nzzTeJj49n//79ho9/VvfRqNvojC/D1WhOS3lLD12EMxwar2YKGvfSWootMgsQlDX3jnu4xjjaegfpOjaELSpHdSDUcPvtt/POO++Mf+AkcSimIIToEUJ0j75sQogRIYSx5QE1Z6S8uZcAP4ElYSa0HjJbju/RVo5/glpAWNbcY7IY36J89O/tnzATehrA1mWyIvO54IILiI6Odtn4DqWkSimPJwgL1evtWmCpq0RpPklZcy+ZsWFY4vJhz99UBzBdqdM9HOuE3iZCkmYSesCP0iZtFNzJWGXgqIw5sAcV30nzkFX9bz8ITfuMHTNxDlxhfDe1iTDhtQZS8S/gMhfo0ZyGsuYetXAnLg8Ge6G73mxJvkOz8tlaEgvISbDqmYKbKWvuISLYn8j0OWqHnim7HEdrH1130lsLUAjYXKJI8wk6+wc50tHPLUvST9TdaSlRNWE0rqdxtBhw0jzyEppYe6gFKaXHN2/3FvY3dDM7ORIxLQP8gz0rrmDyE72rcHSmcPVJr8tQpbOvdZUozQn21ikf6rzUyBNlmxt3m6jIx2jcC+GJYE2gICWS9r5B6o4eM1uVTzAwPMKhhm7mpkWCxU+5Vup3jf9BjVM4unjtjpNeX5FS/lRK2eJqcRrYW9uJEFCQGgkhUaq+/BFj2+9pzkLjHkhSy3EWT1fBvR3Vur+UOyhp7GFwxM781Ci1I+0caPgYhnzbSXHLLbewbNkySktLSU1N5emnnzZ0fEezj2YIId4QQrQKIVqEEK8JIWYYqkRzWrZVtZMbbyViLDc+fSnU7gC7rgPjcnpb1bqQdJVTkZdgJSLYXxsFN7HzsPo7z08fNQrpy2BkUBlqH+b555+nsbGRoaEh6urquOuuuwwd31H30d+BF4EkIBl4CXjeUCWaT9FtG2JHdQcX5sef2JmxHAa6oP5j84T5CtUb1TZzJQAWi2B5dixrS1oYHrGbp8tHWFfSQm5COEmRIWpHxrkg/KDsXXOFeTmOGoVQKeWzUsrh0ddzQLArhWng9T0NDNslq2afVHMndxX4BcL+l80T5iscfA1CYyFp3vFdn1mQQlvvAB8c0t5TV9LUZWNHdQeXzDzp2g+Nhhkr1LU/MmyeOC/HUaPwthDiQSHEdCFEhhDiP4A1QohoIYTrVlH4IFJKNpe38vjGSn79Xinz06JYkBZ14oDgSJh5DXz8DNTtArt+YjUcKVXcpuQtmHcz+J1I0rsoP56suDC+9899/PStg/xrdz12+xnLg2kmSO/AMLtqjvKfr+3HLqXKujuZwrug8wi8932ViaSvf8NxtJ/CjaPbr56y/2ZUwTwdXzCIH752gGe3qQoiM2LD+PWN8z6d/njxD5Vr46mLwD8ECu+Ey36qF7Q5i90Oa38Eu/6sVs6GJ8B5//6JQwL8LDz+hUK++2oxf/2ohsERO4cau/nuat2m01l21Rzljj/voNumZgHfvSKftFP7h+RfCQu/CNv/qF5R6XD7GohKc5tOT09JltK5hxTh7ABnHVyIy4GHUb0YnpJSnjaxVwhxPfAysFhKWXS2MQsLC2VR0VkPmbIcauzmioc388VlGdx3aR4RIf5nvvj62qD0bah4X7k5bnwGZuksYafY+RS8dR/M+gxMPw/yVkPkmXso2O2SB14p5tXd9Xz4wIUnfN+aCWO3Sy75zUaG7HZ+eNVschPCyYgJO/MH2iuhehO8+33lUv3cX9yis7q6GqvVSkxMjEcaBikl7e3t9PT0kJmZ+YmfCSF2SSkLxxvD0cVrAcDXgAtGd20AHpdSDp3lM36oPs6XAnXATiHE61LKg6ccZwXuBXw+z/K5bTWEBPhx36V5RIaOU4kzLBYWfgHm3wqPLoZt/6eNgjNICR89AmlL1ReMAze8xSL41sU5vLSrjlc/rueeC7Ndr9NL2VTeSlVbHw/fPJ9LZznQtyImS706KmHrH6CvHcJcX9o8NTWVuro6WltbXX6uyRIcHExq6uQXtzrqPvo/IAD4w+j7L4zu+/JZPrMEqJBSVgEIIV5ALXg7eMpxPwF+AdzvoBavRErJhtJWVuTGjW8QTsbiB/NvgXX/Az1NYE10nUhvpm4nHD0MKx6YkBsuLTqUgpQINpa2aqPgBO8dbCY8yJ8rCpIm9sGCG5QxL31LuZVcTEBAwKeewL0NRwPNi6WUX5JSrht93QGMV5UqBTi5TVjd6L7jCCEWAmlSyrccVuylVLf1Ud95jPNyYif+4exL1bZqg6GafIrKdYCAvCsm/NGVufHsOnKUHtsZJ86acfiwvI2lM6IJ9J9gObakeWBN1te+gTj6PzAihMgaezO6cM2p1VNCCAvwG+A+B469WwhRJIQo8uRpmzN8VNkOwHnZkzAKiXMhJFr5WDWTo+YjSCyAkGkT/mjh9GmM2CX76nVZ58lQ29HPkY7+yV37QkD6OWpBp8YQHDUK9wPrhRAbhBAbgHWM/2VeD5ycEpA6um8MK1AAbBBCHEaV4n5dCPGpQIiU8gkpZaGUsjAuLs5ByVOLkqZurMH+ZMSEjn/wqVgskFoIDbom0qQYHlTuo4zlk/r43NEyDMV12ihMhgMNqjXLgvSJG2QAUpeo3uXdDQaq8l0cNQpbgMcBO9Ax+u+t43xmJ5AjhMgUQgSi0ldfH/uhlLJLShkrpZwupZwObAOuGS/7yFspbeohL8E6+YyGpPmqJMNgn7HCfIG2Uhjqh9TJ1emPDgskLTqE4rpOg4X5BmXNPQgBOQnhkxsgbYna1vnkV4fhOGoUngEyUUHhR1DrEp492weklMPAN4B3gUPAi1LKA0KIHwshrpm8ZO9DSklJUw95idbxDz4TyQtA2qHJ+J6tXk/LaI3+sSq0k2BOSqR2H02S0qYe0qNDCQ10NO/lFOJnAQJaTs1h0UwGR/8XCqSUJ98x64UQ4/4PSCnXAGtO2ffDMxy70kEtXkdjl40e2zD5zhiFsVIMTcXKx6pxnJaDYAmAmMlnD+UlRPD2/ib6B4cn/+Xmo5Q0dZOX4MS1HxgK0ZnaKBiEozOFj4UQx9tvCiHOAfRczSBKR7t55SVGTH6QiGQIijzx1KtxnJZDEJsD/oGTHiI3IRwpobJFu+8mgm1ohMPt/c49EIGaLTRro2AEjhqFRcBHQojDo0HhrcBiIcQ+IUSxy9T5CGN9f516WhJCtetsLTVIlQ/RfPBEV7tJMuYP1+06J0ZFSy8jdkmuEUahoxKGdAMkZ3F0nnu5S1X4OKVNPSRGBE9s0drpiM9XRdw0jmPrhq4jsOhLTg2TERNGgJ843mhe4xhjRtTpmULCLBVTayv7RFVbzcRxyChIKWtcLcSXKXU2yDxG3ExVPbW3FcK9M3XXcMZmVk4EmUEVypsRG065nilMiNKmHgL9LUw/W50jR4ifrbbNB7VRcJIJLh/UGM3wiJ2K1l6DjEKe2rZ6UHNzT2csOOmk+wiUC6msRRuFiVDS1EN2XDj+fk5+FUXPUH1GdLDZabRRMJnD7X0MDtudiyeMMfbFpo2C47QcgoAwiMpweqjcBCu1HcfoH9QNYBzFsFmynz/E5ulECwPQRsFkSpuUD9qQG8OaBEER2ihMhJYDaoZlcf5WyB0NNlfouIJDdPUP0dRtM+baBxVX0EbBabRRMJnSpm4sArLjJ7ma82SEgLh81ZFK4xgth9SXiQGMpRSXNGkXkiOcSMU2yCjEz4TuOtUgSTNptFEwmZKmHqbHhhEc4GfMgPH5yq/qwuZJXkNvK/S1nghSOkl6dCghAX6UNGqj4AilTarmkSGuUziRLKBnC06hjYLJHGrqZmaSE4vWTiVhDhzrgJ5G48b0VloOqK1BMwU/iyA30UrJ6Jed5uwcbOwhMiSApMhgYwYci6npYLNTaKNgIt22IWo7jjHLSKOQWKC2ugbS+BhQ8+hUZiZaOdTY7XSfXF/gUGM3M5OcKAJ5KpFpEGjVMwUn0UbBRMbcDIYahYSxfO19xo3prTQfgNBYCI83bMj8RCtH+4do7RkwbExvZMQuKTF6liyEmi1oo+AU2iiYyKFG5WYw9MYIjlTplU3aKIxLy0HDXEdj5I/+Xx7Sweazcri9D9uQ3dgHIlBGofmAjqk5gTYKJrK/vovosEASIoKMHThxjnYfjcfwgDKcBq9+nTmagbRP91Y4K2ONdQx9IAJIKFAxta46Y8f1IbRRMJGimqMsTI8yzqc6RtJ8aK+A/g5jx/UmGothZBDSjC0zHhkaQF6Cle3V+m9/Nj6uOUpIgJ9x6ahjpI8Wcz4yXg8wzZnQRsEkWnpsVLf1sXh6tPGDZ54PSKjZYvzY3kLtdrVNXWL40OfMiGZXzVGGRuyGj+0t7KjuYGFGFAHOlrc4lYTZagFnzUfGjutDaKNgEutLWgBYPplm5eORvFCVbqjaaPzY3kL1RlUvx5pg+NDLZsTQPzjCzsN6tnA6WrptHGrqZmlmjPGDW/zUbKFqvY4rTBJtFEzijb2NpEeHMjvZYJ8qqGYxM1bAoTdgRNfh+RS2LqhcD3mrXTL8yrx4woP8eXmX9mufjjX7GpESrpiT6JoTzLwGjh6G+o9dM76Xo42Cm+noG+S1PfV8WNHGTYvTjI8njLHgNuhtgn0vuWb8qcyuv4B9CAquc8nwIYF+3LAoldf2NPBSUS1NXTa9bmGUbtsQT2+pZm5qJNnxBscTxph1jVqvsO4nMDzomnN4MbqZrJvo6h/irr/upKjmKKCyLu5YPt11J8y9HFIK4fVvqC/BxDlw4fcg1AUxjKlAfwdsfUylK1Z8ADmrIGWRy07375fksqO6g/tfVo0J56dF8exdS7AGO9lIaYry148O8/SH1TR12RiRkl/d4MKeB8GRcMl/wZrvwENpyk14yY8gd5XrzulFiKn2BFNYWCiLiqZee+j7XtzLa3vq+X+X5JAVF87KvHhCAg2qd3Qm+tph86+hYbcKrE4/D774mlrk40uMDMFTl0BTsSoYmDQPLvuZyw3k8IidvXWdbKvq4H/fK+XL52Xy/SuNXRcxFVhf2sIdf97JkunRzE+P4rLZiSzKmOb6E5d/oGIL5e+pFNWvb4NpzpdIn6oIIXZJKQvHO07PFNxAc7eNf+2p5/Zzp/ONi3Lcd+KwGLj8Z+rf25+At++H6k0q3uBL7H8VGvfADX+Cguvddlp/PwuLMqJZlBFNaVMPL+ys5b5VecYVP5wi/O6DctKjQ/nbV84xPtvobORcol5LvwYPz4dtf4ArfuG+809RdEzBDbyxt4ERu+S2pSY+pSz8oppW733BPA1mUfQniMmBWZ81TcLnClPpsQ2zobTFNA1mcKS9n721ndy2NN29BuFkIlMh/0rY/wrYR8zRMIXQRsENbCxrJTs+nMxYJ/vQOkNAsMq2KV2j3Cm+Ql+7cp0VXG9II53JsmxGDOFB/mwubzNNgxm8e6AJgCsKkswVMusaVSa9doe5OqYA2ii4GNvQCDuqOzg/xwXrESZKziqwdarVvL5C5TpAQs6lpsrw97NwTmY0WyvbTdXhbj6qbCM7Ppy06FBzhWRdBAg4vNlcHVMAbRRczL76LgaG7Zyb5QFGIeNctT3iQ6s9D29SbrPkBWYrYVlWDFVtfTR328yW4haklBTXdTE/LcpsKRAyTZVI1yudx0UbBRdTXKdaA85LjTRZCWBNhGmZUONDdWEadqsV3hbzg7sL0lXGzd5a3yiW19Blo71v0DOufYCMZcp9pBd0nhVtFFzM/vouEiKCiI8wqLuUs6QvhfpdZqtwD0PHVG39lIVmKwFU3ww/i2BfvW/0EC4eNX5zUj1gpgCqztVQH7SVma3Eo9FGwcUU13UyJ8VDnpRA5ej3NkFPk9lKXE/zAbAPe4TrCNRK55z48OOzR2+nuL6LAD/BzCQXrVyeKGNl0hv3mqvDw9FGwYX0DgxT1dbHnBQPeVKCk24MHwg2N+xWWw8xCgBzUiLZX9/lE2Uvius6yUu0EuRvvusOgNgcCAjVRmEctFFwIQcbupES5qS6oOjdZEmco7a+cGPUfwxhcRCRYraS48xJjaS9b5D6zmNmS3EpY0HmuZ7iOgIVV0qc4xvXvhO41CgIIS4XQpQKISqEEA+e5uffFkIcFEIUCyHWCiG8ag168Wj3rQJPch8FWSEmW63w9XbGgsweVNZjzJW438vjCofb++mxDXtOkHmMpHmq3Ild97o4Ey4zCkIIP+Ax4ApgFnCLEOLUwi+7gUIp5VzgZeCXrtJjBvvru0iMCCbe6iFB5jGS5nm/+2igF9pKPcp1BKoQoi8Em8ceiDzKdQrq2h/shY4qs5V4LK6cKSwBKqSUVVLKQeAF4NqTD5BSrpdS9o++3QakulCP2ymu7/KsWcIYSfOg64h3t+ts2gfSDsnzzVbyCYID/MhNsLKvvttsKS6luK6LIH8LuQnhZkv5JMdjaj4wU54krjQKKUDtSe/rRvedibuAt0/3AyHE3UKIIiFEUWtrq4ESXUe3bYiq1j7Pmz6Db9wYx4PMnpGOejJzUiLYV9fp1cHmfXVdzE6OwN+sekdnIi4f/INPXB+aT+ER/2NCiNuAQuBXp/u5lPIJKWWhlLIwLi7OveImyf7RtMO5nrCa81SSRp+evXm9QsPHKsDsgnabzjInJZKj/UNeG2weHrGzv8HDgsxj+AWoh6K6qVd+31240ijUA2knvU8d3fcJhBCXAN8HrpFSDrhQj1spHvUZe9QahTFCoiA2F+q82Sjs9rh4whhji7m8Ndhc1txL/+AIC9I90CiAaq7UuNe3CkNOAFcahZ1AjhAiUwgRCNwMvH7yAUKIBcDjKIPgVTWFi+s6SYsOITos0GwppyelEOqLvLO5ua0L2is8Lp4wRn6iFX+L8NpFbLtrVXfBBWluaKQzGVIWwfAxaDlothKPxGVGQUo5DHwDeBc4BLwopTwghPixEOKa0cN+BYQDLwkh9gghXj/DcFMKKSV7jnR65vR5jNRFqpRwZ43ZSozHAxetncxYsHmPl9ZA2nOkk+iwQNKiQ8yWcnrG2rB6s/vUCVzaeU1KuQZYc8q+H57070tceX6zONzeT0OXja9lenA/5NTFalu7E6ZNN1WK4VRvBuGnat14KMuyYnh2Ww22oRGv68RWVHOUBWlRCA9aH/IJpk2H0Bh17RfeabYaj8MjAs3exoflKkPqvBwPDorHz1blhCvXma3EeKo3qSJ4wR60kvwUzsuJZXDYzs7D3pUWXNvRT3VbH8uzPaBU/JkQAjJXQMUHehHbadBGwQWsLWkhdVoI02NMbixyNvz8IftSKH/Xu1oU9ncot0CmZ/ehPiczmiB/y/HOZN7C2kPNAFyQ68EPRAC5l0NfCzTq1NRT0UbBYOo7j7GprJXPLkjx3OnzGHlXQH87VK03W4lx7H8F5AjMunb8Y00kNNCfKwoSeW1PA90278iCkVLy0q46ZidHkB3vYYvWTiXnUvALhD1/N1uJx6GNgkG09gywobSF+1/ai7+fhZsWp43/IbPJvwqsSbD+5zA8aLYa5xnsh21/UEXPxgr/eTBfPn8GfQPDfPsfe9hS0UbDFF23IKWkudvGk5urONDQzW1Lp0AJs9BomHMj7H5OB5xPwaWBZl9ASslD75Tw5KYq7BIsAn50bQGp0zzYdTSGfyCs+h945S74TT5EZcDsz8K53/SoInJnxT4CO56AsnegtQx6GuEL/5wS+gtSIvnuFTP5+duH+OCQysi+58Is7r8s32RljrPzcAffeWkvNe2qWs2SzGg+t2iKVKu5+D+hci08eREEhEHaEvjsH1WHQh9GTLWl9oWFhbKoyHNWI76xt4FvPr+b6xemcmNhKrkJVqZ56tqEM1GyBkreUgXk6nbCZT+DZfeYrcox3vsBfPQIJBRAdCbM/7xyi00hOvoGKWnq5sWdtfxrTwMv3L2UpTNizJY1Lm29A1z4vxuIDgvk9nOnk5dopTAjmkD/KeSA6GmG/S/D0cPKlZRQAHe8DZYp9Ds4iBBil5SycLzj9EzBSR5ZV05+opVf3jAXP4vnP52elvzV6iUl/P0mWP8zmH+ryk7yZNor4aNHYeGX4OqHp8Ts4HREhwVyblYsC9OnsaWyncfWV0wJo/Dkpir6B0f459cXe34M4UxYE048ACXNg9fugdI1MPMqc3WZiPeZQzdypL2fsuZeblqcNnUNwskIARd9X5UW3vuC2WrGZ+dTYPGHC78/ZQ3CyQQH+HHrknQ+rGijudtmtpxxefdAE+fnxE5dg3Aq826BiFQoetpsJaaijYITrC9VfuAL8+JNVmIgSfPUa/8rZis5O/YRKP6HeqLzwKJ3k+XqeclICW8WN5ot5axUtfZyuL2fi/K96Nq3+MG8m6Fqo3eXlR8HbRScYH1pCzNiw5geG2a2FGPJv1pVkexpNlvJmanbqdJpZ14z/rFTiOz4cHITwllf4tmlwNaXqgWaXvVABJC3WqU0l79vthLT0EZhkhwbHGFrZTsrve2mAMi7HJBqYZunUroGLAGQfepv1D4AAA/3SURBVLHZSgzn/Jw4dhzuwDbkuYsK15e0kBMfTlr0FMiymwjJCyA8AUrfMluJaWijMEm2VrUxMGznwnwPX7k5GRIKlG+1/D2zlZyZ8vch41wI9sDS5E5y/mgJjB3VnunC6BsYZnt1Oxd6k+toDIsFclZB5XqfLa2tjcIkWV/SSmigH0s8uejdZBECci5RvlVPvDF6mlTZ46yLzFbiEs7JjCHQz8KmMs/sMvhhRRtDI9L7XEdjZF8CA90+24hHG4VJIKVkXUkLy7NjCfL3rgqXxxm7MWp3mK3k01RtUNusC02V4SpCAv0onD6NDyvazJZyWtaXtGAN8qdwuoenLE+WGStVld2KD8xWYgraKEyCipZe6juPee+TEqiCchZ/z7wxqjao0scJnl/KYrJckBtHSVMPLR6WmiqlZH1pC+fnxhLgaf2XjSIkSpWW98Rr3w146f+qa1k3mhnilfGEMYIjIG0pVHhYFoaUyt+bucIrV52OcX6OKj29qdyzZgsHG7tp7h7wzgSLk8m+BBr3QK9nuvBciffeVS5kfWkL+YlWkiI9tLOUUWRfDE37lA/fU2gthd4mNcX3YmYmRhAbHsjmcs/6Utowmoq6Ms+LH4jgRFabN1UQdhBtFCbI0b5Bdh4+ysUzvfxJCdTTEnhWI56xKf2MlWaqcDkWi+D8nDg+LG/Dbvec+mTvH2xmbmok8dZgs6W4lqT5ykXpgy4kbRQmyNqSFkbskstm+0AlxcQ5KmfbkxbyHHpDpcxOmwLlmZ3kgtxY2vsG2e0hvZybumzsqe30jWvfYoGsi6Firc91Z9NGYYK8tqeelKgQ5qR4X378pxBCzRYq18HwgNlqoLsBard73SrmM3HJzASCAyy8vKvObCkAvFncAOAbRgHUtd/fplbPewK7/uqWGIc2ChOgvLmHzeVt3Lw4zfO7qhlFwfVg6/SMWki7/qK2c24wVYa7sAYHsHpOEq/vqaet11yjPDRi59ltNSzKmOY9BfDGI381BEXAjsfNVgJ1u+CNb7nlPtRGwQF6B4apau3lgVeKiQj259Zz0s2W5D6yLoL42bDhIbB1maejqx62/Z/qlRCTZZ4ON3PPhdkMDNv5zkt7OdTYTVOXze0xhhG75OEPyqlp7+ffVvjO354gKyy6HQ78E6o3m6fDbocP/kuVsp9/q8tPp/spnIWiwx08+Oo+Klp6AdVV7ZFbFhITHmSyMjciBFz5a/jLlfDYOZC+FKJnwJzPQfxM1567/mPY+/yo22iHqoy66n9ce04PIysunP++Zjb/+dr+45k/M2LDePTWhcxKjnDZeTv7B/mv1w+wobSVrmNqVftnF6RwiS8kWJzMBferBlTPXAvJ89W1n3WxqqbqSm/BsU5Vvr6tTGUA1u2Aq36nUsVdjO68dgbqjvaz6rebiLMGcdPiNOLCg1iQHkV2vNXl5/ZIDm+BrY9B6yHVpcovEL74mjISruDga/DSHeAfBFHpMG06rHgAUha65nwezuG2PvbVd9HeO8AfN1YxbJe8fe/5xFmNf0AZGrFz4+Nb2V/fxXULUkmICCIvMYLLCxK9o2/IROlrh49+r9YttJaqlq8X/IfqPeIKuurgz6uhswaCo1Qf9UW3wzlfdcoQOdp5TRuFM/AfL+/lX3saWHffiqnRb9md9LbAny6D4UH4xg4INLh0eOcR+MO5EJ8Pt73ilUXvnKGkqZurH/mQ6xem8tD1cw0f//GNlfz87RIeuWUBV89LNnz8KY3dDq9/E/Y8B3d9AGmLjT/Hs9epmfHnX4KMZYYN66hR0DGF01DV2ssrH9dz2zkZ2iCcjvB4uPYP0F0HH/7O+PE3/AJGBuH6p7VBOA35iRF8/pwMXiyqpbK119Cx+waG+ePGSlbmxWmDcDosFrjiIfX0/vZ/qBX2RlK9CSrXwsoHDDUIE0EbhdPw2w/KCfSz8PULfSioNlEylsHs62Dro8aueG6rgL1/h8Vf9om1CJPlGxdlE+TvxyNryw0d99ltNRztH+Lei3MMHderCLLCRT+Aho+Vm9MopIQP/hsiUtT1bxLaKJzCocZu3tjbwB3LpxPrSwHlyXDRD9QT/cZfGDfmxofAPxjO+3fjxvRCYsOD+OKyDF7f23A8EcJZ+geHeXJTFRfkxrEg3UsroBrFvFsgLh/W/QRGho0Zs+RNqN8FKx+EAPNK6GijcAq/eb8Ma7A/X71AzxLGJSYLFt2hFtW0VTg/XksJ7HsZltwN4V5eW8cAvnLBDIL8/Xh0nTGzhee21dDeN8i9F2cbMp5XY/GDi/8L2itg97POj2cfgbU/gdhcmOf6tNOzoY3CSeyo7uD9g83cff4MIkMDzJYzNVjxH+rJft2PnR9rw89V0Prcbzk/lg9w8mzB2dhCj22IJzZVcX5OLIsyvLBxlCvIu0JVEt7wEAz2OzfWx3+FtlK46D/Bz9yVAi41Cv+/vfuPraq84zj+/lBooQXLkB9BqKUgIg3KVIYIE3+AG0OdixHFTGUGNZsDdc4taJbhljlctixuizMaxjQbgzE0ij822JRsZmH80i1YKhZEoNoKIpYOVkrpd388p5cr046UnnvoPd9XQnrO6ek93yfc9nvPc57n+0iaJmmLpK2S5n3M94sk/T76/lpJw+KMpz0Hm1v4zjObGFzak1svGp5UGF1P74EwcW7oWz2RlareXAWbn4EL50DJqZ0XX55ru1uY/2wVRzo4qc3MePCFavYeaObez43q5AjzmARTHwhVe9c+2vHXaayHl38A5ZNg9FWdFV2HxZaSJBUAjwCXA7XAekkrzGxz1mmzgX1mdoakmcCPgOvjiinbxh37WLPtfZpbWkFiVVU9W3f/myduGU+vwjxdTS0uE+fAhkWw7GaYtiDcAhf1gV79oPATRm+1NIeJOY114RZ89QIYWAkX3ZPb2Lu4/r2LmH9VJfOe3sT1j61hfEU/DOjVo4BLRw3k7KH/O3pr9/4mVlbVU9fQxMHmI+zYe4DVW/bw1YtHMLasb+4b0ZWVXwijroDVPwR1g/LPhhFzvfqGD0wfp7UVDuyBI4fC8Os/zoPDTWGS6ElQPie2eQqSLgQeMLPPR/v3AZjZgqxzVkbnrJHUHagHBlg7QZ3IPAUzY+OOfTz8l5rMUodSeOg/vH8J908fzdTKQR167dSrfx2W3hDe5G269YCKyXDOdXDWldDaEurTVz8PNavCcp9tBp0NNyyBvmW5jz0PLFu/i1+srqG+oQkhDre2YgaXnTWQey4/kzFDSqndd5CFr2znd+t20tzSSkE3UVxYQO+i7swYV8bdU0bSLY2T005U0354anZ4T2cbWBlqh509IwxhrV0XqvxWPwf73zl6XlEpzFh0tFR9TBKfvCbpWmCamd0a7d8EXGBmc7LOeT06pzba3xad84nLTXU0KSxdt5Ofv1TDuw1N9Csp5GsXj2Dm+DJ6F3Wn1UjnTM3OduRwqCjZWAfNB+D9mtAllJ0oINSpHzU9rIlQWganDIZThub1Smq51nDwMIvX7eCxv75Fw38OM7i0J+/tb6KbxLXnD+W2ycMZ3r8kPYUdc2HPm7BvOxxqDKVZ3ngBdv3jo+d07xnKZFRMDs/Pik+FYZNyMh/neJNCl6h9JOl24HaA00/vWDG6AX2KGF/Rj3HD+nHNeUMoLjza9AL/vegcBT2gfOJHj039Xih3/fYrYcRG2QQouyDxh2n5rrS4B3dccgY3Tihn6bqdVNc1Un5qMdeNK+O0vnm+YmBSBpwZ/rWZdCfs2xHuDA7tD+uTDL8Uik7uKrOp6j5yzrm0OhnKXKwHRkqqkFQIzARWHHPOCmBWtH0t8HJ7CcE551y8YruHN7MWSXOAlUABsMjMqiR9H9hgZiuAXwG/kbQV+ICQOJxzziUk1o5dM3sRePGYY9/N2m4CZsQZg3POuePnwz2cc85leFJwzjmX4UnBOedchicF55xzGZ4UnHPOZXS5NZol7QF2dPDH+wOfWEIjT3mb08HbnA4n0uZyM/u/C5V0uaRwIiRtOJ4ZffnE25wO3uZ0yEWbvfvIOedchicF55xzGWlLCo8nHUACvM3p4G1Oh9jbnKpnCs4559qXtjsF55xz7UhNUpA0TdIWSVslzUs6nrhJKpO0WtJmSVWS7ko6plyQVCDpNUnPJx1LLkjqK2m5pDckVUfrmOQ1Sd+I3tOvS1oiqWfSMXU2SYsk7Y5Wp2w71k/SnyXVRF8/Fce1U5EUJBUAjwBfACqBGyRVJhtV7FqAb5pZJTAB+HoK2gxwF1CddBA59DPgT2Z2FjCWPG+7pCHAncA4MxtDKMufjyX3nwCmHXNsHvCSmY0EXor2O10qkgIwHthqZm+ZWTOwFLg64ZhiZWZ1ZvZqtN1I+GMxJNmo4iVpKHAFsDDpWHJBUikwmbAuCWbWbGYfJhtVTnQHekWrNRYD7yYcT6czs78R1pjJdjXwZLT9JPClOK6dlqQwBNiVtV9Lnv+BzCZpGHAusDbZSGL3MPBtoDXpQHKkAtgD/DrqMlsoqSTpoOJkZu8APwF2AnVAg5mtSjaqnBlkZnXRdj0wKI6LpCUppJak3sBTwN1mtj/peOIi6Upgt5ltTDqWHOoOnAc8ambnAgeIqUvhZBH1o19NSIinASWSbkw2qtyLli2OZehoWpLCO0BZ1v7Q6Fhek9SDkBAWm9nTSccTs0nAFyW9TegevEzSb5MNKXa1QK2Ztd0BLickiXw2FdhuZnvM7DDwNDAx4Zhy5T1JgwGir7vjuEhaksJ6YKSkCkmFhAdTKxKOKVaSROhrrjaznyYdT9zM7D4zG2pmwwj/vy+bWV5/gjSzemCXpFHRoSnA5gRDyoWdwARJxdF7fAp5/nA9ywpgVrQ9C3g2jovEukbzycLMWiTNAVYSRissMrOqhMOK2yTgJmCTpH9Gx+6P1s12+WMusDj6sPMWcEvC8cTKzNZKWg68Shhh9xp5OLNZ0hLgEqC/pFpgPvAQsEzSbEKl6OtiubbPaHbOOdcmLd1HzjnnjoMnBeeccxmeFJxzzmV4UnDOOZfhScE551yGJwXn2hFVIb0j2j4tGg7pXN7yIanOtSOqG/V8VJHTubyXislrzp2Ah4AR0QTAGmC0mY2R9BVClcoSYCShSFshYcLgIWC6mX0gaQShbPsA4CBwm5m9kftmOHd8vPvIufbNA7aZ2aeBbx3zvTHANcBngAeBg1FhujXAzdE5jwNzzex84F7glzmJ2rkO8jsF5zpudbRWRaOkBuC56Pgm4JyoQu1E4A+hTA8ARbkP07nj50nBuY47lLXdmrXfSvjd6gZ8GN1lONclePeRc+1rBPp05Aej9Su2S5oBoXKtpLGdGZxznc2TgnPtMLO9wN+jBdR/3IGX+DIwW9K/gCryfBlY1/X5kFTnnHMZfqfgnHMuw5OCc865DE8KzjnnMjwpOOecy/Ck4JxzLsOTgnPOuQxPCs455zI8KTjnnMv4L6baU5HA335sAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plot_population(res)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmUXHWd9/H3t6qrF9KddCARsyBBjBEUhBjjMsiDDyibsrjG5XlcD48LI6iDQ0YfhmHGBTnj43I4o6gc0YPiBhg1mnEbt1GgIYSwRUJAkk6EsCRdTXqp5fv8UXU71ZWq7ttddWu5+bzOyemqW7+69c3t2/dbv+X+fubuiIiIhJFodgAiItI+lDRERCQ0JQ0REQlNSUNEREJT0hARkdCUNEREJDQlDRERCU1JQ0REQlPSEBGR0DqaHcBMLViwwJctW9bsMERE2srtt9/+uLsvrHU/bZc0li1bxsDAQLPDEBFpK2b213rsR81TIiISmpKGiIiEFmnSMLMzzGyLmW01s0srvP5OM9ttZncW/703ynhERKQ2kfVpmFkSuBp4FbADuM3M1rn7vWVFv+vuF0YVh4iI1E+UNY3VwFZ33+bu48ANwLkRfp6IiEQsytFTS4DtJc93AC+pUO71ZnYy8Bfgw+6+vUIZkZZ288ZBrtqwhZ17Rljc38Mlp6/gvBOXNDsskbprdkf4j4Fl7n488AvgukqFzOwCMxsws4Hdu3c3NECR6dy8cZC1N25mcM8IDgzuGWHtjZu5eeNgs0MTqbsok8YgcETJ86XFbRPc/Ql3Hys+/Rrwoko7cvdr3H2Vu69auLDme1NE6uqqDVsYyeQmbRvJ5Lhqw5YmRSQSnSiTxm3AcjM7ysw6gTXAutICZrao5Ok5wH0RxiMSiZ17Rma0XaSdRdan4e5ZM7sQ2AAkgWvd/R4zuwIYcPd1wIfM7BwgCzwJvDOqeESisri/h8EKCWJxf08TohGJVqTTiLj7emB92bbLSh6vBdZGGYNI1C45fQVrb9w8qYmqJ5XkktNXNDEqkWi03dxTIq0mGCW19sa7GMnkWdjbxcfPPkajpySWmj16SiQWzjtxCauWHQrAv7/phUoYEltKGiJ1kh7NAjA8lm1yJCLRUdIQqZP0aGbST5E4UtIQqZOgphH8FIkjJQ2ROgmapZQ0JM6UNETqIJvLs2+8MORWSUPiTElDpA5KO7+Hx9SnIfGlpCFSB6W1C9U0JM6UNETqYKhkxJSShsSZkoZIHQSJojOZIK37NCTGlDRE6iBIGov6u3WfhsSakoZIHQSJYtG8bjVPSawpaYjUQTB6anF/D8NKGhJjShoidRDULhbP62EkkyOTyzc5IpFoKGmI1MHQaIbOjgSH9XYCqLYhsaWkIVIH6dEsc7s76O0qLFGjmW4lrpQ0ROogPZqlrztFX3cKmHzfhkicKGmI1EF6NENvVwd93cWahpqnJKaUNETqoFDT2J80NOxW4kpJQ6QO0qOZYtIoNE+lNWmhxJSShkgdBH0aEx3hqmlITClpiNTBcFnz1JCShsSUkoZIjfJ5Z3i8UNPoTiULkxYqaUhMKWmI1Gh4PIs7zC3WMnq7O7QQk8SWkoZIjYJaRdA01dfdoZqGxJaShkiNghlue7tSxZ8d6giX2FLSEKmRahpyMFHSEKlRUNPYnzRSmkZEYktJQ6RG+2saheapvq4OTVgosaWkIVKjIGnMVfOUHAQiTRpmdoaZbTGzrWZ26RTlXm9mbmarooxHJAoH1DS6UwyPZXH3ZoYlEonIkoaZJYGrgTOBY4G3mNmxFcr1ARcBt0QVi0iU0qMZOhJGd6rw59Tb3UEu74xkck2OTKT+oqxprAa2uvs2dx8HbgDOrVDuX4ErgdEIYxGJTDDDrZkBaKZbibUok8YSYHvJ8x3FbRPMbCVwhLv/NMI4RCJVmOE2NfE8mLRQSUPiqGkd4WaWAD4HfDRE2QvMbMDMBnbv3h19cCIzkB7NTiQKgLnB9OgadisxFGXSGASOKHm+tLgt0Ae8APgvM3sYeCmwrlJnuLtf4+6r3H3VwoULIwxZZOaC5qmAmqckzqJMGrcBy83sKDPrBNYA64IX3X2vuy9w92Xuvgz4M3COuw9EGJNI3Q2VN08FS77qXg2JociShrtngQuBDcB9wPfc/R4zu8LMzonqc0UabXgsO3GPBuwfeqvmKYmjjumLzJ67rwfWl227rErZU6KMRSQqap6Sg4nuCBepgbszPJad1Dw1p1NJQ+JLSUOkBvvGc+TyPqmmkUwYvV2aSkTiSUlDpAblU4gEeru0ep/Ek5KGSA0mFmDqntw9qEkLJa6UNERqMFS2AFNASUPiSklDpAZBTWNuWdLo7U6R1n0aEkNKGiI1qNanUahpqE9D4kdJQ6QGwV3f5c1Tc9U8JTGlpCFSg/3rg1cYPaWkITGkpCFSg/RoloTBnM7kpO193SlGMjkyuXyTIhOJhpKGSA2CadGDBZgCwVTpqm1I3ChpiNSgfIbbQJ9mupWYUtIQqUH5ZIWBIJEMaQSVxIyShkgNCku9Vkoaap6SeFLSEKlBoaZRvXlKw24lbpQ0RGpQmBa9evNUWpMWSswoaYjUoFqfhkZPSVwpaYjMkrsX+zSqN08NKWlIzChpiMzSWDZPJucVaxpdHQlSSVOfhsSOkobILA1VmUIEwMzo605pISaJHSUNkVmamOG268CaBmhNDYknJQ2RWUpXWYApoEkLJY6UNERmqdoMtwHVNCSOlDREZml4mppGX3dK04hI7ChpiMzSdM1TfV0dmrBQYkdJQ2SWpho9Vdiu5imJHyUNkVkKEkJvldFTvd2Fmoa7NzIskUgpaYjMUrAAUzJhFV/v606RyzsjmVyDIxOJjpKGyCxVmxY9oJluJY6UNERmKahpVBO8pqQhcaKkITJL6bGpaxpzg+nRNexWYiTSpGFmZ5jZFjPbamaXVnj9fWa22czuNLM/mNmxUcYjUk/VFmAKqHlK4iiypGFmSeBq4EzgWOAtFZLCt939OHc/Afgs8Lmo4hGpt+Eqa2kEeoMlX3WvhsRIlDWN1cBWd9/m7uPADcC5pQXcfajk6RxAYxOlbQxNW9NQ85TET6ikYWbPNbNfmdndxefHm9knpnnbEmB7yfMdxW3l+/6gmT1IoabxoXBhizRfejTD3KlqGuoIlxgKW9P4KrAWyAC4+13AmnoE4O5Xu/vRwD8CFRORmV1gZgNmNrB79+56fKxITcazecay+ambp5Q0JIbCJo1D3P3Wsm3T/SUMAkeUPF9a3FbNDcB5lV5w92vcfZW7r1q4cOG0wYpEbboZbgGSCaO3S1OJSLyETRqPm9nRFPsczOwNwK5p3nMbsNzMjjKzTgo1k3WlBcxsecnTs4EHQsYj0lTTTSES6O3q0Op9EitTn/H7fRC4BniemQ0CDwFvn+oN7p41swuBDUASuNbd7zGzK4ABd18HXGhmp1Fo9noKeMcs/x8iDTXdDLcBTVoocRMqabj7NuA0M5sDJNw9HfJ964H1ZdsuK3l80QxiFWkZYZqnCq8raUi8TJk0zOwjVbYD4O66r0IOSumxcDWN3u4Ue0fUPCXxMV1No6/4cwXwYvb3SbwWKO8YFzloBLWHuSFqGjue2teIkEQaYsqk4e7/AmBmvwNWBs1SZnY58NPIoxNpUfubp6bp09DoKYmZsKOnDgfGS56PF7eJHJQmRk+F6AgfVtKQGAk7euqbwK1mdlPx+XnAddGEJNL60qMZelJJUsmpv3f1dacYyeTI5PLTlhVpB2FHT33SzH4GvKK46V3uvjG6sERaW3qayQoDwX0cw6NZ5s/pjDoskciFShpm9izgceCm0m3u/khUgYm0svRodtqmKdjf5zE8pqQh8RC2eeqn7J+Btgc4CtgCPD+KoERa3dBoZtp7NGD/fRxDmulWYiJs89Rxpc/NbCXwgUgiEmkDw2PZKWe4DUzUNNQZLjExq545d78DeEmdYxFpG2H7NLR6n8RN2D6N0jvDE8BKYGckEYm0gfRohr6u6ZunJqZH16SFEhNh+zT6Sh5nKfRx/LD+4Yi0h/A1jUJiUfOUxEXYpHGvu3+/dIOZvRH4fpXyIrGVzeXZN54L2RFe+BMbUtKQmAjbp7E25DaR2BsOOVkhQFdHglTS1KchsTHdLLdnAmcBS8zsiyUvzWX6lftEYinsFCJQmBG6rzulhZgkNqY763cCA8A5wO0l29PAh6MKSqSVBfdchBlyC1pTQ+JlulluNwGbzOx6d9dZL0Lpqn3T92lAcclXJQ2Jiemap77n7m8CNpqZl7/u7sdHFplIixoOudRrQDUNiZPpzvpgOdbXRB2ISLsI7rkIX9NIaSEmiY3pmqd2FX/+tTHhiLS+9AxrGnO7OyZGXIm0u+map9Lsn6gQwIrPDXB3nxthbCItaaZJQ81TEifT1TT6pnpd5GA0NJqhsyNBV0cyVPneYk3D3TGziKMTiVbYO8KDmW1PolDT+IMWYZKDVXo03Ay3gb7uFLm8M5LJcUhn+PeJtKJQd4Sb2WUUlnc9DFgAfMPMPhFlYCKtKj2anZiIMAzNdCtxEvbMfxvwQncfBTCzzwB3Av8WVWAirSodcgGmwMRMt6NZDlcvoLS5sHNP7QS6S553AYP1D0ek9Q2HnOE2MLeYYNJavU9iIOyZvxe4x8x+QaFP41XArcF8VO7+oYjiE2k56dEsyxYcErp8r5qnJEbCJo2biv8C/1X/UETaw0ybpyaWfNW9GhIDYdcIvy7qQETaRdgFmAJ9ap6SGAk7euo1ZrbRzJ40syEzS5vZUNTBibSafN4ZHs/OuiNcpN2F/br0eeB1wGZ3P2DiQpGDxfB4Fvfw06KDkobES9jRU9uBu2eaMMzsDDPbYmZbzezSCq9/xMzuNbO7zOxXZnbkTPYv0mgTCzDN4D6NZMLo7dJUIhIPYc/8jwHrzey3wFiw0d0/V+0NZpYErqYw0moHcJuZrXP3e0uKbQRWufs+M3s/8FngzTP8P4g0TNAvMZPmKSiuqaHV+yQGwtY0Pgnso3CvRl/Jv6msBra6+zZ3HwduAM4tLeDuv3H3YM7oPwNLwwYu0gwzXUsjoEkLJS7CnvmL3f0FM9z3EgrNWoEdwEumKP8e4Gcz/AyRhprpDLeBXiUNiYmwNY31ZvbqqIIws7cDq4Crqrx+gZkNmNnA7t27owpDZFpDs2ye6utOkdZ9GhIDYZPG+4Gfm9nIDIbcDgJHlDxfSoWpR8zsNODjwDnuPlb+OoC7X+Puq9x91cKFC0OGLFJ/QW1hJqOnIGieUp+GtL+wN/f1mdmhwHImz0E1lduA5WZ2FIVksQZ4a2kBMzsR+Apwhrs/FjpqkSbZ3zw1w5qGRk9JTIRKGmb2XgrrhS+lMLvtS4H/Bk6t9h53z5rZhcAGIAlc6+73mNkVwIC7r6PQHNULfL+4OM0j7n5ODf8fkUilRzN0JIzuVNhKekFfd8dEJ7pIOwtbx74IeDHwZ3d/pZk9D/jUdG9y9/XA+rJtl5U8Pm0GsYo0XTCFyExX4OvrTjGSyZHJ5UklZ5ZwRFpJ2LN3tGQtjS53vx9YEV1YIq0pPZqZmLV2JoKbAZ9WZ7i0ubBn/w4z6wduBn5hZk8Bf40uLJHWlB7N0tc1s/4MmLx6X/8hnfUOS6RhwnaEn198eLmZ/QaYB/w8sqhEWlR6bGYz3AaC9wxpBJW0uRmf/e7+2ygCEWkH6dEsS/p7Zvy+YLSVOsOl3alHTmQG0qOZGd+jAZObp0TamZKGyAzMdAGmwMT06Jq0UNqckoZISO7O8NjMFmAKqHlK4kJJQySkfeM5cnmvsSNcSUPam5KGSEgTCzDNIml0dSRIJU19GtL2lDREQprtAkwAZkZfd0oLMUnbU9IQCSmY2nw2zVOAlnyVWFDSEAlpttOiBzRpocSBkoZISLU0TxXep5qGtD8lDZGQZrvUa6C3K6VpRKTtKWmIhFRrTWNudwfDmuVW2pyShkhI6dEsCYM5nclZvV/NUxIHShoiIaVHs/R2zXwBpkBvsabh7nWOTKRxlDREQhoazcy6aQoKzVq5vDOSydUxKpHGUtIQCWl4lpMVBiYmLVQTlbQxJQ2RkGY7w21A06NLHChpiISUHquteWpu8b1pDbuVNqakIRJSrTWNXtU0JAaUNERCqlfzlO7VkHampCESgruTrsPoKVDzlLQ3JQ2REMayeTK52S3AFNDoKYkDJQ2REII5o/q6lDTk4KakIRLC/skKZ988lUwYczqTShrS1pQ0REIYrnGG24BW75N2p6QhEkI9ahqF92vSQmlvShoiIeyfFr22mkavkoa0OSUNkRBqXYAp0NedmlhrXKQdRZo0zOwMM9tiZlvN7NIKr59sZneYWdbM3hBlLCK1GKpxAaZAoXlKfRrSviJLGmaWBK4GzgSOBd5iZseWFXsEeCfw7ajiEKmHoKbRW8OQWygM2VXzlLSz2v4CprYa2Oru2wDM7AbgXODeoIC7P1x8LR9hHCI1CxZgSiZmtwBToK+7Y2Iklkg7irJ5agmwveT5juI2kbaTHs3UXMuAQvPWSCZHJqfvSdKe2qIj3MwuMLMBMxvYvXt3s8ORg9DwWG2TFQaCxPO0OsOlTUWZNAaBI0qeLy1umzF3v8bdV7n7qoULF9YlOJGZqHWG24AWYpJ2F2XSuA1YbmZHmVknsAZYF+HniUSm1hluA0HSGNIIKmlTkSUNd88CFwIbgPuA77n7PWZ2hZmdA2BmLzazHcAbga+Y2T1RxSNSi/rVNAqJR53h0q6iHD2Fu68H1pdtu6zk8W0Umq1EWtrQaLauNQ01T0m7aouOcJFmS49mmFvHjvC0Ji2UNqWkITKN8WyesWxezVMiKGmITCuY9qM+92kEHeFKGtKelDREpjE8Vp9p0QG6OhKkkqY+DWlbShoi06jXDLcAZqaFmKStKWmITKNeM9wGejVpobQxJQ2RadSzphHsRx3h0q6UNESmESSNuXWqaWjJV2lnShoi06jXUq+B3q6UphGRtqWkITKNiQWY6pQ05nZ3TIzIEmk3Shoi00iPZuhJJUkl6/PnouYpaWdKGiLTSI9m61bLgEKNZXgsi7vXbZ8ijaKkITKNdJ0WYAr0dafI5Z2RTK5u+xRpFCUNkWmk6zTDbWBi0kI1UUkbUtIQmUa9ZrgNaHp0aWdKGiJTuHnjIJt37OX3DzzO333m19y8cVYrFk8S3O+R1rBbaUNKGiJV3LxxkLU3biabL3RYD+4ZYe2Nm2tOHL2qaUgbU9IQqeKqDVsO6KweyeS4asOWmvYbNE/pXg1pR0oaIlXs3DMyo+1h/enBJwD4wPV31K3JS6RRlDREKvjvBx/HrPJri/t7Zr3fmzcOcuXP7594Xq8mL5FGUdIQKeHufOOPD/G/vn4rC3q76OqY/CfSk0pyyekrZr3/qzZsYTSTn7StHk1eIo2ipCFSNJbN8bEf3MXlP76XV654Br/66P/gytcfz5L+HgxY0t/Dp193HOeduGTWnxFVk5dIo9Rv8LlIG3t0aJT/863buXP7Hj506nIuPnU5iYRx3olLakoS5Rb39zBYIUEcOqezbp8hEiXVNOSgt/GRp3jtl/7AXx5N8+W3r+Qjr3ouiUSVDo0aXXL6CnpSyUnbDHji6XE++r1N7Nk3HsnnitSLahpy0Ll54yBXbdjCzj0jzDskRXokw+L5PXzzPS/nec+cG+lnB7WW4PMX9/dw8WnLefiJp/nyb7fx2788xuXnPJ+zj1uEVeuJF2kia7eZNletWuUDAwPNDkPaVHDDXun9FwmDfz33BbztpUc2MTK4d+cQ//jDu9g8uJfTjjmck5Yfxld/99BEcrnk9BV1bSqTg4uZ3e7uq2rej5KGxEVpDaL8IpvN5dk8uJd3XHsrQxXuxF7S38MfL/2fjQ75ANlcnmv/+BCf/fn9ZCcPsqInlay5I14OXvVKGmqeklgor0EM7hnhH394F7++/1HSo1lue/ipKe/AbpXRSx3JBBecfDRf/8NDPDo0Num1kUyOK39+f8WkMVXCFKknJQ1pirAXuenKuTt7RzJ8cv19B0z5MZbNs27TLp69cA7nnrCYlx19GP/20/v4297RAz6nlhv2ovBYWcII7No7yhmf/x0nPms+LzpyPiuf1c+m7Xv4p5vunpQw1964GUCJQ+pOSeMgNZNvpvW6wJeWK68VVLrI3XTHDtbetHniZrjBPSP8w/c38YPbt9PVkWTHUyPseGofT49XX8zIgF9/9JSJ59mcH9CnUesNe1GoNjS3r7uDZ8zt5iebdvKdWx8BCn0y+bJW5sINgwfWSqL4vcvBpe36NLoWLfdVF32loRe5uO2zUmdwtfbysGUrletOJfjn1x7LK1cczvBYln3jWZ4ey3Hht+/giacPHFranUpw/NJ+9u7L8NS+cR5LV/62bcAxi+ayZH4PS+f3sHT+IVz9m608WWGflfoq2uFiON1xz+edrbuHuf2vT00k3EqWHXYIi+b1sKi/m+HRDL/ZsptMbv/ffHcqwafPP47zVy6d0eeXl23mOS/htEVHuJmdAXwBSAJfc/fPlL3eBXwTeBHwBPBmd394qn12LVrui97x+Ugucu2yz+5Ugn957fM56/hF5POQcyfvzvq7dvKp9fczWtKD2tWR4MJXPoeXP2cB2VyebN750Hc2Vrxoz+tJcfFpy8nmnPFcnmzO+drvt5Gu0BfQ3ZHgJc8+jNFMjrFsnrsH905MIV6L1UcdSn9PivmHdPLdge0Vyxjw0GfOnrRtJse+XYS9cP7dZ35dsVbS29XBKSsWsmvvKLv2jLCzQrNcYGFfF/09Keb1pOg/JMUftz5RcTnaw+Z08qW3nsiczg7mdCX5/QOPc+XP7580NUoj/45mcpyandyavc83nX7S+PhjD3VV3NEMRJY0zCwJ/AV4FbADuA14i7vfW1LmA8Dx7v4+M1sDnO/ub55qv0HSgMLF85TnPgPHcQcHfv+X3ZMumhPv60jwsqMPmyh3y7YnGKtSbuWz5uMUjos7bHxkD+O5A8umksZxS+bhxXL37Nw76VtcoCNhLFswh3zeJy7wO58aJVfh2BvQ05kk704+T8XPbRUvXDqPrlSSro4Ev3/g8arlPnX+cczpSjKns4NDupJc9J072T18YC2ivFZQ7WJYbaTTwfrtNOxF9qhLf0q1v/a3rD6CPfsy7B3JsGdfhnt3DdUUU8JgyfweOpMJujqSdHYkuHfnUMXzeU5nkjWrn0VnR4KujgTX/uGhiiPc5h+S4tOvO45kIkEyAclEglu2PcHX/vAQ42VflC4+bTmvfv4z6UgYyYTxy/se5TNlX6i6Uwk+dd4LOH/l0kn3xLTal8R67XPb1/6esV0P1HzzT5RJ42XA5e5+evH5WgB3/3RJmQ3FMn8ysw7gb8BCnyKo0qQBsOLwvuK+Cs/v/1u6akzHL52HFQtv2r6narkXL5uPYYUrOHDrQ09WLfuK5QsmHk914SzcrAXJhJE048YpZjV970lHkUwYZsaXf/tg1XKfOPsYzIxkcb//90f3VC37jXe9mFQyQSqZ4IPX31Hxov3Mud387KJX0JG0ibInf/Y3oS7cM7nAR/WN82AWJmHO5HdUrezCvi6+uObEQlPjeI4PfWdj1ZjOP3EJ49k8Y9lCbXSqv485nUnGc/mKX7oaIVlMLh0JY2Q8VzG5JhPGkv4eEgYJMx55cl/F2nUqaRy7eN5EuYTBpu17KybMro4EJz1nAVYs97sHdh8woSUUzvuzj180sU8z40d3DrKvQn/enK4kb139LMwMA66/5RGGx7Lsuu7iuiSNKDvClwCl7Qs7gJdUK+PuWTPbCxwGVD+7St/c38OGD588adtUfxjrLjwpVLnvv+/loff5rfe8JFS5q9+2ctK2Wx56smrZT7zm2InnP960s2q5977i2ZO2ffm326qWPWXFMyaef/zsYypejC8983nML5sD6ZLTV4TqOA5bDirfFV3pIhe2nBBqjqyZ/I6qlf34WcfwsqMPm9h25c/ur3rO/b83nzBpW5iklc87J13564pNac/o6+K6d68ml3eyeSeXz/OG//hT1drTF9acMFH2Yz+4q0opuOjU5ZP2+dXfP1SxXC7vvOjI+YVWAIdtjz9dsVwm5/T3pMh7oQUk7161xWAsm+fR9Cj5fKFcpYQBhYENf3rwieJnO7k8FRMGwNNjOa6/5ZGJz6/UolKLthg9ZWYXABcAdD7zOcDMT/ZaLnJx2+dMLsZRXeDDTgRY7wkDD2ZR/N7rfX4mEsbHznhexXL/dNYxHLNo8jQv1UaZLenv4dwT9sf6hV8+ULXch1/13Enb1m/+W6hEeMdfn6pa7rp3r560baqE+ZO/f0WocmFrg2FbAWYryqQxCBxR8nxpcVulMjuKzVPzKHSIT+Lu1wDXQKF5akmDL3Jx22dQNuzFWBf4+Kj3773Z53zcvtBFuc96ibJPo4NCR/ipFJLDbcBb3f2ekjIfBI4r6Qh/nbu/aar9ahoRESnV7FFJ7bLPlh89BWBmZwGfpzDk9lp3/6SZXQEMuPs6M+sGvgWcCDwJrHH3bVPtU0lDRGTm2mLuKXdfD6wv23ZZyeNR4I1RxiAiIvWjRZhERCQ0JQ0REQlNSUNEREJT0hARkdCUNEREJLS2mxrdzNLAlmbHEcICQk6H0mSKs37aIUZQnPXWLnGucPe+WnfSFtOIlNlSj7HGUTOzAcVZP+0QZzvECIqz3topznrsR81TIiISmpKGiIiE1o5J45pmBxCS4qyvdoizHWIExVlvB1WcbdcRLiIizdOONQ0REWmSlk0aZnaGmW0xs61mdmmF17vM7LvF128xs2VNiPEIM/uNmd1rZveY2UUVypxiZnvN7M7iv8sq7asBsT5sZpuLMRwwisIKvlg8nneZ2cpK+4kwvhUlx+hOMxsys4vLyjTtWJrZtWb2mJndXbLtUDP7hZk9UPw5v8p731Es84CZvaPBMV5lZvcXf6c3mVl/lfdOeX40IM7LzWyw5Hd7VpX3TnldaECc3y2J8WEzu7PKext5PCtehyI7P9295f5RmEr9QeDZQCewCTi2rMwHgC8XH68BvtuEOBcBK4uP+yisH1Ie5ynAT1rgmD4MLJji9bMtyO9BAAAEpUlEQVSAn1FYGf2lwC1N/v3/DTiyVY4lcDKwEri7ZNtngUuLjy8FrqzwvkOBbcWf84uP5zcwxlcDHcXHV1aKMcz50YA4Lwf+IcR5MeV1Ieo4y17/d+CyFjieFa9DUZ2frVrTWA1sdfdt7j4O3ACcW1bmXOC64uMfAKeaWc2Lps+Eu+9y9zuKj9PAfRTWPW9H5wLf9II/A/1mtqhJsZwKPOjuf23S5x/A3X9HYc2XUqXn4HXAeRXeejrwC3d/0t2fAn4BnNGoGN39P909W3z6ZworaDZVlWMZRpjrQt1MFWfxWvMm4DtRfX5YU1yHIjk/WzVpLAG2lzzfwYEX44kyxT+KvcBhNEmxeexE4JYKL7/MzDaZ2c/M7PkNDWw/B/7TzG63wprr5cIc80ZZQ/U/xlY4loHD3X1X8fHfgMMrlGml4/puCrXJSqY7PxrhwmIz2rVVmlJa6Vi+AnjU3R+o8npTjmfZdSiS87NVk0ZbMbNe4IfAxe4+VPbyHRSaWV4IfAm4udHxFZ3k7iuBM4EPmtnJTYpjSmbWCZwDfL/Cy61yLA/ghbp+yw5FNLOPA1ng+ipFmn1+/AdwNHACsItC008rewtT1zIafjynug7V8/xs1aQxCBxR8nxpcVvFMlZYj3we8ERDoithZikKv6jr3f3G8tfdfcjdh4uP1wMpM1vQ4DBx98Hiz8eAmyhU9UuFOeaNcCZwh7s/Wv5CqxzLEo8GTXjFn49VKNP042pm7wReA7ytePE4QIjzI1Lu/qi759w9D3y1yuc3/VjCxPXmdcB3q5Vp9PGsch2K5Pxs1aRxG7DczI4qfvNcA6wrK7MOCHr63wD8utofRFSK7ZpfB+5z989VKfPMoK/FzFZTOOYNTW5mNsfM+oLHFDpH7y4rtg7431bwUmBvSdW2kap+g2uFY1mm9Bx8B/CjCmU2AK82s/nFJpdXF7c1hJmdAXwMOMfd91UpE+b8iFRZ/9n5VT4/zHWhEU4D7nf3HZVebPTxnOI6FM352Yje/VmOCDiLwiiAB4GPF7ddQeHkB+im0ISxFbgVeHYTYjyJQpXvLuDO4r+zgPcB7yuWuRC4h8JIjz8DL29CnM8ufv6mYizB8SyN04Cri8d7M7CqCXHOoZAE5pVsa4ljSSGR7QIyFNp930OhD+1XwAPAL4FDi2VXAV8ree+7i+fpVuBdDY5xK4U26+D8DEYcLgbWT3V+NDjObxXPu7soXOwWlcdZfH7AdaGRcRa3fyM4J0vKNvN4VrsORXJ+6o5wEREJrVWbp0REpAUpaYiISGhKGiIiEpqShoiIhKakISIioSlpiIhIaEoaIiISmpKGSI3M7Bwz+2HZtveb2ZeaFZNIVJQ0RGr3SeCfy7Y9CBzThFhEIqWkIVIDM3shkHD3u83sSDN7f/GlFC08663IbClpiNTmBOD24uNXAcuLj4+lMPeQSKwoaYjUJgH0mlmSwnTZfWbWA7wT+HYzAxOJgpKGSG3WU5jV9E7gy8DzgQHgGi8uwSkSJ5rlVkREQlNNQ0REQlPSEBGR0JQ0REQkNCUNEREJTUlDRERCU9IQEZHQlDRERCQ0JQ0REQnt/wMiZknde7pBLwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def plot_spectrum(pulse, tlist, xlim=None):\n",
"\n",
" if callable(pulse):\n",
" pulse = np.array([pulse(t, None) for t in tlist])\n",
"\n",
" dt = tlist[1] - tlist[0]\n",
" n = len(tlist)\n",
"\n",
" w = np.fft.fftfreq(n, d=dt/(2.0*np.pi))\n",
" # the factor 2π in the normalization means that \n",
" # the spectrum is in units of angular frequency,\n",
" # which is normally what we want\n",
" spectrum = np.fft.fft(pulse) / n\n",
" \n",
" # we assume a real-valued pulse, so we throw away\n",
" # the half of the spectrum with negative frequencies\n",
" w = w[range(int(n / 2))]\n",
" spectrum = np.abs(spectrum[range(int(n / 2))])\n",
"\n",
" fig, ax = plt.subplots()\n",
" ax.plot(w, spectrum, '-o')\n",
" ax.set_xlabel(r'$\\omega$')\n",
" ax.set_ylabel('amplitude')\n",
" if xlim is not None:\n",
" ax.set_xlim(*xlim)\n",
" plt.show(fig)\n",
"\n",
"\n",
"plot_spectrum(H[1][1], tlist, xlim=(0, 20))"
]
}
],
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment