WeakZeeman.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extremely weak-field Zeeman effect\n",
    "**PHYS 475, Lecture 16, 2019.02.22**\n",
    "\n",
    "Instructions:\n",
    "(To answer questions, add a \"markdown\" cell and type your answers there.)\n",
    "\n",
    "1. Check the constants to verify that they are correct for the hydrogen ground states\n",
    "2. Derive the value of the atomic unit for magnetic field.  (The answer is given below.)\n",
    "2. Below in the cell beginning \"next steps\" is a function called np.kron().  What does this do?  Where have we seen this before?\n",
    "3. Find the eigenstates in the small-B limit and the large-B limit.  (Numerically, what are \"small B\" and \"large B\"?) Explain how to get these results from perturbation theory.\n",
    "   * Bonus: If you have time, _do_ the perturbation calculations in the code and compare with the \"exact\" values done here.\n",
    "4. Modify the code such that B is in the x direction.  How does this change the eigenvalues?  How does this change the eigenstates?\n",
    "5. The most common isotope of the alkali metal potassium (K-39) has nuclear spin I=3/2 and ground state hyperfine splitting of 462 MHz (in frequency units).  For K-39 the nuclear g-factor is gN = -1.4e-4.  (*Error:* this number is a factor of mp too small.) Adapt the code to calculate the hyperfine states of K-39.  What is the value of A for this state?\n",
    "6. Create a plot showing the B dependence of the energy eigenvalues of K-39 for fields in the z direction.\n",
    "7. The homework problem looks at deuterium.  What changes in the code to do this atom?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your name here*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:13.958471Z",
     "start_time": "2019-02-24T15:07:13.441867Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "np.set_printoptions(precision=3, floatmode='maxprec')\n",
    "#assuming atomic units throughout unless specified"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:13.973945Z",
     "start_time": "2019-02-24T15:07:13.962968Z"
    }
   },
   "outputs": [],
   "source": [
    "#constants\n",
    "alpha=1./137.; # fine structure constant\n",
    "mp = 1836.; # proton mass\n",
    "muB = 0.5; #Bohr magneton\n",
    "muN = muB/mp; #nuclear magneton\n",
    "mu0 = 4*np.pi*alpha**2; #permeability of vacuum\n",
    "gs = 2.0; # electron spin gyromagnetic ratio\n",
    "psi0 = 1/np.sqrt(np.pi); #value of wavefunction at origin\n",
    "h = 4.136e-15; #Planck's constant in eV s\n",
    "EH = 27.2; #Hartree energy in eV\n",
    "B0 = 2.35e5; #atomic unit of magnetic field in T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.001394Z",
     "start_time": "2019-02-24T15:07:13.978938Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = 2.16e-07 (a.u.) = 1422.41 MHz\n"
     ]
    }
   ],
   "source": [
    "# Properties of this atom (hydrogen)\n",
    "S=1/2; I=1/2;\n",
    "gN = 5.59;# nuclear gyromagnetic ratio (depends on atom and isotope)\n",
    "#hyperfine amplitude coefficient (depends on atom and state)\n",
    "A = 2*mu0/3*gs*muB*gN*muN*psi0**2 # McIntyre Eq. 11.10\n",
    "print(f\"A = {A:.2e} (a.u.) = {A*EH/h/1e6:.2f} MHz\") # display A in MHz to compare with McIntyre Fig. 11.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.023850Z",
     "start_time": "2019-02-24T15:07:14.011873Z"
    }
   },
   "outputs": [],
   "source": [
    "#degeneracy of S and I.\n",
    "#Must round so that these are stored as integers (the extra 0.2 is to ensure the rounding goes to the right number)\n",
    "degS = np.int(2*S+1.2); degI = np.int(2*I+1.2); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.057288Z",
     "start_time": "2019-02-24T15:07:14.034837Z"
    }
   },
   "outputs": [],
   "source": [
    "#create lists of possible Sz, Iz values in decreasing order\n",
    "Szlist = np.arange(S,-S-1,-1);\n",
    "Izlist = np.arange(I,-I-1,-1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.080253Z",
     "start_time": "2019-02-24T15:07:14.061281Z"
    }
   },
   "outputs": [],
   "source": [
    "#Generating the individual spin matrices in the uncoupled basis\n",
    "S2 = S*(S+1)*np.matrix(np.eye(degS)); #S-squared matrix\n",
    "I2 = I*(I+1)*np.matrix(np.eye(degI)); #I-squared matrix\n",
    "Sz = np.matrix(np.diag(Szlist)); #Sz matrix\n",
    "Iz = np.matrix(np.diag(Izlist)); #Iz matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.120769Z",
     "start_time": "2019-02-24T15:07:14.084251Z"
    }
   },
   "outputs": [],
   "source": [
    "# list of the basis state  labels |Sz,Iz>\n",
    "basislist = np.column_stack((np.diagonal(np.kron(Sz,np.eye(degI))),np.diagonal(np.kron(np.eye(degS),Iz)))); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.157881Z",
     "start_time": "2019-02-24T15:07:14.124264Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.5  0.5]\n",
      " [ 0.5 -0.5]\n",
      " [-0.5  0.5]\n",
      " [-0.5 -0.5]]\n"
     ]
    }
   ],
   "source": [
    "print(basislist) # display the row labels (Sz, Iz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.199799Z",
     "start_time": "2019-02-24T15:07:14.161872Z"
    }
   },
   "outputs": [],
   "source": [
    "# Generate other spin matrices\n",
    "Sp=np.matrix(np.diag(np.sqrt(S*(S+1)-(Szlist*(Szlist+1))[1:]),1)); #S+ ladder matrix\n",
    "Ip=np.matrix(np.diag(np.sqrt(I*(I+1)-(Izlist*(Izlist+1))[1:]),1)); #I+ ladder matrix\n",
    "Sm=Sp.H; #S- ladder matrix\n",
    "Im=Ip.H; #I- ladder matrix\n",
    "Sx=(Sp+Sm)/2;\n",
    "Ix=(Ip+Im)/2;\n",
    "Sy=(Sp-Sm)/(0.+2j); #Sy matrix (\"j\" is sqrt(-1))\n",
    "Iy=(Ip-Im)/(0.+2j);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.235133Z",
     "start_time": "2019-02-24T15:07:14.206287Z"
    }
   },
   "outputs": [],
   "source": [
    "# next steps: build hyperfine operator matrix and zeeman operator matrix\n",
    "# We are assuming l=0 states here (otherwise HF is more complicated)\n",
    "\n",
    "# hyperfine operator = A*(S dot I)\n",
    "HF = A*(np.kron(Sz,Iz)+0.5*(np.kron(Sm,Ip)+np.kron(Sp,Im)));\n",
    "# Zeeman operator, except for the magnetic field magnitude\n",
    "ZnoB = gs*muB*np.kron(Sz,np.eye(degI))+gN*muN*np.kron(np.eye(degS),Iz); #multiply this by B to get the Zeeman operator\n",
    "# magnetic field magnitude in atomic units:\n",
    "B=0; # Change this value to test different magnetic field strengths\n",
    "# Total hamiltonian\n",
    "H = HF+B*ZnoB;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.268354Z",
     "start_time": "2019-02-24T15:07:14.239126Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 5.407e-08  0.000e+00  0.000e+00  0.000e+00]\n",
      " [ 0.000e+00 -5.407e-08  1.081e-07  0.000e+00]\n",
      " [ 0.000e+00  1.081e-07 -5.407e-08  0.000e+00]\n",
      " [ 0.000e+00  0.000e+00  0.000e+00  5.407e-08]]\n"
     ]
    }
   ],
   "source": [
    "# print the Hamiltonian matrix\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.299363Z",
     "start_time": "2019-02-24T15:07:14.272346Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The eigenvalues are [-1.622e-07  5.407e-08  5.407e-08  5.407e-08].\n",
      "The eigenvectors are the columns of the matrix:\n",
      "[[ 0.     1.     0.     0.   ]\n",
      " [ 0.707  0.     0.     0.707]\n",
      " [-0.707  0.     0.     0.707]\n",
      " [ 0.     0.     1.     0.   ]].\n"
     ]
    }
   ],
   "source": [
    "#calculate and display the eigenvalues and eigenstates of H for this value of B\n",
    "evs, ess = np.linalg.eigh(H) \n",
    "print(f\"The eigenvalues are {evs}.\")\n",
    "print(f\"The eigenvectors are the columns of the matrix:\\n{ess}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.324308Z",
     "start_time": "2019-02-24T15:07:14.305343Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1422.411016622688"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# maximum frequency splitting in MHz\n",
    "(max(evs)-min(evs))*EH/h/1e6 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.358742Z",
     "start_time": "2019-02-24T15:07:14.329796Z"
    }
   },
   "outputs": [],
   "source": [
    "# Generate data for the plots by stepping through B values.\n",
    "Bmin = 0; # minimum B in atomic units\n",
    "Bmax = 5e-7; # maximum B in atomic units\n",
    "Bsteps = 101 # number of steps\n",
    "Blist = np.linspace(Bmin, Bmax, Bsteps) # generate list of B values to evaluate\n",
    "Elist = []; #empty list to store energy eigenvalues\n",
    "for B in Blist:\n",
    "    H = HF+B*ZnoB;\n",
    "    evs, ess = np.linalg.eigh(H)\n",
    "    Elist.append(np.sort(evs))\n",
    "Earray = np.array(Elist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:14.996028Z",
     "start_time": "2019-02-24T15:07:14.363233Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl0XFl16P/vLqmkGjRPlmVNnmVbtjzIktwTzdhAoLshzdA0Dd0MDYQQeCG/8CN5KyHJS0LWey+PvF8ISSfk8QgJhIShCSHMdDeElmfZludRtjxpHktSTef3x7kqldSyLckqVUnan7Xuqqpbt27tkqpq1znn7nPFGINSSil1J65kB6CUUmpx0IShlFJqRjRhKKWUmhFNGEoppWZEE4ZSSqkZ0YShlFJqRjRhqDkTkY0iclhEBkXkN0TEKyL/JiL9IvIvCX7u4yLy4Dzt60ERaZ+nfT0lIr9ItX0tZyJySURek+w4lgJNGPNARJ4QkaFpFiMiv5fs+BLot4HnjTHZxpj/DTwGrAAKjTFvm+tOReQzIvKV221jjNlijHl+jvs3IrJuTsGpZWEuPyKWw/tKE8Y8MMb8ozEmK34BPgHcBP42yeElUhVwfMrtM8aYcJLiWTZEJD3ZMcyViKQlOwY1R8YYXeZ5AXYAg8CDcetygS8C14GrwH8D0uLufx9wEugFfgBUxd1ngF8Dzjr7/SNgLfASMAB8Hchwts0Hvgt0Ovv6LlAet6/nncf/p7OvHwJFt3ktbwJagD7gl8A2Z/1PgQgwCgwBXwWCQMi5/f4ZvK4twI+AHmxy/R3g9VP2c+QWcV0CXuNc/4zzN/iy85qOA/W3eNyLzt9z2Nn/O4AHgXbgk0CH8z96Ou4xmcD/AC47cf414L3F/p8CfuFs3wtcBN7g3Pc24OCU7T8JfNu5Xgh8x/mf7nP+T7+Y8j74qPM+uOisuwfYD/Q7l/fEbb/aeb2DwI+BzwNfibu/yfmf9gFHmPx+ne375Ledv9s14ANOrOuc+74EfAH4nvN3fw328/Bl7Pu0DfivgCvu/xkfZ7Wzv/SZxAY86eyzG/jd+PfKNHG/ETjh7Ocq8FuAHxgBos57ZAgoAxqwn7k+57X+JROfu5e9r273+XHu+5TznIPAaeDVyf7uuuN3W7IDWGoLkAecBz41Zf23gb9x3owlzhfCh5z7HgXOAZuAdOfD88u4xxrsF0kO9kt2DPgJsMb54J0A3utsWwj8KuADsoF/wflCcu5/3olvA+B1bn/2Fq9lJ/YLtBFIA97rfPgy4/b1gbjtp37Qb/m6nNiuY78wPc7txun2c4vYYl8Czvajzoc/DfhToPk2j419mTm3HwTCwB8Cbmc/ASDfuf9zzt+/wInz34A/vcW+n8Imuw86sXwE+yUq2MTTA2yK2/4w8KvO9a9hE58fqMV+mUxNGD9y4vA6l73YL8h04HHndqGz/UvYxJUB3IdNRF9x7luF/UJ9I7an4bXO7eI5vE9eD9zAvjd9wD/w8oTRD9zrPJcHmyyec/6e1cAZJn5kTPr/M33CmDY2YDP2C/sB5+/9587/9lYJ4zpwv3M9H9gZ955on7LtLmySTXdiOgl84jbvq1t+foCNwBWgLO41rk3299cdv9+SHcC8vyD4e+ef1DoP+3ol9tfB+DIKPHqb7cX5EDwHSNz6FdgveW/cuseBnznX/2P8w+LcdmG/sKri3oj3xt1/kLiEBPxP4HO3iGk70Bt3+3ngv8bd/jXg+7d47BeAP5qy7jTwirh93S5h3PJ1Oa//8C2ed9J+brHNJSYnjB/H3bcZGLnNY6dLGCM4X0jOug7ny0GwvxrXxt23B+cX/jT7fgo4F3fb5zxfadzf9I+d61uwX/CZ2C+UEFAT99g/4eUJ41Vxt58E9k15/pecGCqxX5S+uPu+wkTC+BTwD1Me+wMmfnjM5n3y98QlUGAdL08YX467Pw37edgct+5D2PGw6d5H1bw8YUwbG/B7wNfi7vNjW6y3ShiXnefOmbL+QaYkjGke+wngW7d5X93y8+P8jTqwrS337Z4nlZalOIbxJewvnrtmjPmZMWa7MWY78Crsl90Pb/OQT2F/Gb7XOO8ORxX2l+t1EekTkT5sa6Mk7v6/iLuvB/tFtSpuHzfjro9MczsLQER8IvI3ItImIgPYpnLelH7jG3HXA+OPnUYV8MnxuJzYKrDN85m43euqwP5KnC9TX5Nnlv383Wby2Mv436UY+6V/MO51fN9Zf8dYjDEB5+r43/j/Au8SEcF+4X/dGDPm7C8d+6tzXNs0+46/v2yabdqwf98yoCfu+ac+tgp425T/7X3AyuleB7d/n5RN2feVabaJX1eEbfXExz4e90zdKrZJsRhjhrEtp1v5VWwrq01EXhCRPbfaUEQ2iMh3ReSG89n6E+e13MotPz/GmHPYhPMZoENEviYiM/1cJc2SSxjGmBexX0wxIrJWRL4vIgdF5OciUjOHXT8G/MeUD2D8czyI7S99zBjTN+XuK9hfVEXGmDxnyTHGbIm7/0Nx9+UZY7zGmF/OIc5PYpu7jcaYHGzTHOwX9Wxdwf4ajo/LZ4z56iwef6vXdQU7DjMdc4v1ydCFTchb4l5DrrEHNsyaMaYZ+4v3fuBd2O4bsH35YewXyrjK6XYRd/0a9kspXiW2K+s6UCAivrj74vd9BdvCiP/f+I0xn53ta3Keq/wWzzNd3F3Y1lR87ONxg23RxcddOstYYs/vvP7CW21sjNlvjHkE++Pt29guwanxjvsCcApY73y2fofbf65u+/kxxvyTMeY+7N/BAH820xeZLEsuYdzCs8DHjDG7sINafzWHfbwTO7D7MiKyEtv//AljzOGp9xtjrmNbJv9TRHJExOUksVc4m/w18GkR2eLsL1dE5npYajb2C65PRAqA35/jfsAe4fVhEWkUyy8ivyIi2TN8/O1e13eBUhH5hIhkiki2iDQ6990EqkUkUe/Pm9jxnzsyxkSxf4f/JSIlACKySkQeuovn/zJ2wDRsjPmF8zwR4JvAZ5xW4mZsn/ftfA/YICLvEpF0EXkHtjvuu8aYNuCAs78M55fzm+Me+xXgzSLykIikiYjHOZS0/OVPc0dfB54WkU3OF/RtDyV3XuvXgT92/u9VwG86MYHt/n1ARCpFJBf49Cxi+VfgTSJyn4hkYMelpn0fOX+XJ0Qk1xgTwo7xRJy7bwKFzvOPy3a2GXJ+dH5kyi6nvq9u+fkRW8P0KhHJxHZ1j8Q9d8pa8glDRLKwR5L8i4i0YLuCVjr3vVVEWqdZfjBlHyuBrdg+3ul8EDtO8Rfy8lqMv3a2eQ+2GX4C22/9r+NxGGO+hf118TWnqdsKvGGOL/lz2IHALqAZ230yJ8aYA9jX9pdOzOew/eMzffwtX5cxZhA70PpmbPfCWeyYEdiBeoBuETk01/hv4zPA/3W6Cd4+g+0/hX3tzc7r+DG2FTdX/4DtuvyHKet/Hdu1cgPbtfp/brcTY0w39iicT2K7XX4beJMxpsvZ5AnseEs39qi8f8a2dDHGXAEewf5K7sT+Gv5/mMN3gjHmP4D/DfwM+3d6yblr7DYP+xi2JXEBe1TZP2HHQjDG/MiJ9Sh2vO67s4jlOPZIsn/CtjZ6sUfA3cqTwCXn//ph4N3Ofk5hfyBecN4nZdgfm+/CHtX0t06M8T5D3PvqDp+fTOCz2M/pDWwL53dm+jqTRSZ3tS8NIlKN/ZVVKyI5wGljzMrbP+q2+/s4tkvimXkKUS1jIuLFDnjuNMacXcDn/WfglDHmblqdM3meTdgfB5lGa3KWlCXfwjDGDAAXx7tCnKZh3Sx38zi36I5Sag4+AuxPdLIQkd1O16dLRF6PbVF8O0HP9Raniycf26r8N00WS8+SSxgi8lVsk3ijiLSLyPuxTfP3i8gRbFHXI7PYXzV2EO2F+Y9WLTcicgn4OLYbKdFKsYegDmG7jD4y3RjbPPkQtmvrPLYvfmr/vloClmSXlFJKqfm35FoYSimlEmPRTmA2naKiIlNdXZ3sMJRSalE5ePBglzHmdsWowBJLGNXV1Rw4cCDZYSil1KIiItPNKvAy2iWllFJqRjRhKKWUmhFNGEoppWZEE4ZSSqkZ0YShlFJqRjRhKKWUmpGUTRjOdMv7ROSIiBwXkT9IdkxKKbWcpWzCwE6N/CpjTB32NKOvF5GmJMeklFIpIzIwwOBPf8rNz/4ZJhRK+POlbOGec4rTIeem21l04iul1LIVDQQIHDxEYG8zw3v3MXr8OESjiMdD3mO/Sua6dQl9/pRNGABiz0N9EHvC9M8bY/YmOSSllFow0WCQ0SNHGG7ey3BzMyNHj0IoBG433m3bKPrIR/A3NeKpq8OVkZHweFI6YTinctwuInnAt0Sk1hjTGr+NiDwDPANQWTndKZCVUmpxMOEwoydOMNy8l0BzM4FDhzCjo+By4dmyhcL3vgdf0x58O3fg8vnuvMN5ltIJY5wxpk9Engdejz2TV/x9z2LP2U19fb12WSmlFg0TjTJ29iyB5mabJPbvJzpke+IzN2wg7+1vw9/YiG/3btJycpIcbQonDBEpBkJOsvACr8GeyUsppRYlYwzBS5cI7N1rE8TevUR6ewFwV1WS84Y34N/ThK+hgfSioiRH+3IpmzCAldgTqqdhj+b6ujFmxieDV0qpVBC6ds1JDrYVEb55E4D0FSvIeuB+fE178Dc24C4rS3Kkd5ayCcMYcxTYkew4lFJqNsJdXQzv3UugeS/De/cSunwZgLT8fHxNjfgbm/A3NeKuqkJEkhzt7KRswlBKqcUg0t9PYP/+WCti7Ow5AFxZWfgaGih49xP4GhvJXL8ecaVy6dudacJQSqlZmFQL8VIzoydOgDGIx4Nv505yHn4Yf1MTnk2bkPSl9RW7tF6NUkrNs2gwyEhLS6yLaeTIEQiHwe3GV1dH0Uc/ir+xYcFqIZJJE4ZSSsUx4TCjx4/HupgCBw9hxsZsLURtLYVPP4WvsSlptRDJpAlDKbWsmWiUsdOnJ4rlDhwgOjwMQObGjeS/8x02QeyuJy07O8nRJpcmDKXUsmKMIXjx4uRaiL4+ADKqqsh585vwNzm1EAUFSY42tWjCUEoteaGrVxker6Zubibc2QlA+sqVZL3ylfgaG/A3NeEuLU1ypKlNE4ZSaskJd3YyvHdfrFgudOUKAGmFhXaqjaZGmyAqKhZdLUQyacJQSi16kb4+hvfvd45kaiZ47jwArpwcfA27KXjPe/A1NthaCE0Qc6YJQym16ESGhhk5dDDWxTR68qSthfB68e3aRd5b3oKvsQnPphokLS3Z4S4ZmjCUUikvOjbGyOEWhvc2E2jey8ixYxAOI2433vFaiD1NeLduRZZ4LUQyacJQSqUcEwox0toaO5Jp5NAhTDBoayG21lL4vvfhb2rEu2MHLq832eEuG5owlFJJZ6JRxk6dsmeW29vMyP4DRAMBwKmFePxxfE2N+Oq1FiKZNGEopRbceC3E8EsvEWjeS2DfPiL9/QBkrF5NziMP429swtfYQHp+fpKjVeM0YSilFkSw/WrsMNdJtRBlK8l69avxNzXia2zEvWJFkiNVt6IJQymVEKGODgJ798UGqkPt7QCkFRXhb2y0Z5ZrbMRdXq6Hui4SmjCUUvMi3NtLYN9+O1C9dy/B804tRG4uvt31FDz1FP6mRjLWrtUEsUhpwlBKzUlkaIjAgQOxab/HTp2ytRA+H776XeS99a34mhrx1GgtxFKhCUMpNSPR0VFGDh+OjUGMtLZCJIJkZODdvp3i3/gYvsZGWwvhdic7XJUAmjCUUtMyoRAjx1oZbrZHMo0cPowJhSAtDW9tLYUf/AD+pia827fj8niSHe6yFIqGONV9ipbOFt5V8y7SXIltyWnCUEoBYCIRRk+dis3HFDhwEBMIgAiZm2rIf/e7bbHcrnrSsvzJDndZGg4Nc7TzKIc7DnPo5iGOdh1lJDwCQNPKJtbnr0/o82vCUGqZMsYQPH8+dma54X37icbVQuQ9+og9cVDDbq2FSJKukS4O3TxkE0THIU73nCZiIrjExYb8Dbxl3VvYuWInO0p2UOIrSXg8mjCUWiaMMYTa2xlubo4NVEe6ugBwl5WR/epX20NdGxpxr0j8l4+azBjDxYGLHL5pk8PhjsNcGbTTsnvSPGwt3sr7t76fXSW72Fa8jayMrAWPMWUThohUAF8GSoEo8Kwx5i+SG5VSi0vo5s2JM8s1NxO6dg2AtOIi/E1NtliuqYmM8vIkR7r8hCIhTvSciCWIlo4Wesd6AcjPzGdHyQ7esfEdbC/ZzuaCzbjTkn8gQcomDCAMfNIYc0hEsoGDIvIjY8yJZAemVKoK9/ZOKpYLXrwI2FoIf8NuCt7/PvxNTWSsWaO1EAtsIDjAkY4jse6l1q5WxiJjAFRmV3J/+f3sWrGLHSU7qM6pTsn/T8omDGPMdeC6c31QRE4CqwBNGEo5IkNDBGInDnJqIQCXz4d3dz15jz2Gf08TmTU1iMuV5GiXD2MM14evx1oOhzoOca73HAZDuqRTU1DD2ze+nR0lO9hRsoMib1GyQ56RlE0Y8USkGtgB7J3mvmeAZwAqKysXNC6lFlp0ZCRWCzG8t5nR1uMTtRA7d1L8iY/ja2jEu7VWayEWUDga5kzvGQ53HI4tHYEOAPxuP9uKtvG67a9jZ8lOaotq8bl9SY54blI+YYhIFvAN4BPGmIGp9xtjngWeBaivrzcLHJ5SCWWCQUZaW2MD1bFaiPR0WwvxzAfxNzbh3bEdV2ZmssNdNoZDwxzpPEJLRwuHOw5zpPNI7PDWUn8pu0p2sb1kOztX7GR93vqE10cslJROGCLixiaLfzTGfDPZ8SiVaCYSYfTkqYlZXQ9O1EJ4Nm8m/8knbS3Ezl1aC7FAjDHcGL4Razm0dLZwpvcMURPFJS7W563nkbWPsHPFTrYXb2dl1spkh5wwKZswxI74fBE4aYz582THo1QiGGMYO3s2NgYR2L+f6IBtSGeuX2fnY2pswN/QQFpubpKjXR7C0TCne0/HWg/x3UvedC/birfxzLZn2FG8I2mHtyZLyiYM4F7gSeCYiLQ4637HGPO9JMak1F0xxhC6fHmiWG7vPiLd3QC4y8vJeeh1+Bqb8Dc2kF5cnORol4eB4ECserqlo4VjXcde1r1UV1LHjpIdbMjfQLorlb82EytlX7kx5hdA6h1XptQshW7cmFQsF75+HYD0khL899yDv8meFyKjfFWSI136jDFcGbxCS2dLLEGc7zuPwZAmabHq6R0lO9hesp1Sf2myQ04pKZswlFqswj09k4rlgm1tAKTl5eFrbMT3wQ/gb9pDxurUPNZ+KRmLjHGi+0QsORzpPELPaA8A2e5stpVs46Hqh9hRsoOtRVsX7dFLC0UThlJ3KTI4SGD/AdvF9FIzY2fOAODy+/Ht3k3e4+/E39hI5saNWguRYB2BjlhiaOls4UT3CcLRMAAV2RXct+o+tpdsZ3vxdtbmrcUl+v+YDU0YSs1SdGSEwKFDsS6m0dZWiEaRzEy8O3dQ/F/+C/6mRjxbtiDp+hFLlPHah5aOFlo6WzjScYRrw3bqkwxXBrVFtTy56Um2l2ynrriOQm9hkiNe/PTdrNQdmGCQkaNHY11MgSNHYLwWYts2ij78IXzj54XIyEh2uEtW72gvRzqPxOofjncfjw1Ol3hL2F6ynSc2PUFdSV3KzL201GjCUGoKE4kweuJEbKA6cPAgZnTU1kJs2ULhe99jxyJ27sTl11qIRIhEI5zrOxdLEEc6j9A2YMeCxqfWeOv6t7K92LYeSv2lOh60ADRhqGXPRKOMnT03USy3fz/RwUEAMtattfMxNTXi271bayESpG+0j6NdR2npaOFo51GOdR0jEA4AUOApoK64jresewt1xXVsKdqCN92b5IiXJ00YatkxxhBqa4vNxxTYu49Ijz1yxl1ZSc7rH8LX1IS/sZH0osUxKdxiEo6GOdd3jqOdRznSeYSjnUe5NHAJIHZo68NrH6aupI664jrKs8q19ZAiNGGoZSF0/bqtpJ6mFiLr/vvwNTTib2rEvUprIeZb90g3RzuPcrTraKz1MD72UOApYFvxNh5Z94htPRRu0UNbU5gmDLUkhbu7Y7UQw3ubCbVdBiAtPx9fYyP+Zz5oTxxUrbUQ8ykUCXGq5xRHuyZaD1eHrgJ27GFjwUYeXfcodcV1bCvepq2HRUYThloSIgMDBPbvjx3JNHb2LACurCx8DQ0UPPEEvsZGMtev11qIeWKM4drwNY51HuNI5xGOdR3jZPdJgtEgACW+EuqK63i85nG2Fm1lc+FmPOmeJEet7oYmDLUoRQMBAgcPxQaqR0+csLUQHg++nTvJefjN+Jua8GzapLUQ82Q4NExrVyvHuo7ZLqbOo3SP2nmwMtMy2VK4hcdrHqeupI6tRVt1Wo0lSD9JalGIBoOMtLTExiBGjh61tRBuN966bRR9+MP49zThqavTWoh5EI6GOd93nqNdRznWeYxjXcdicy4BVOdUc0/ZPWwr3sbW4q1syN+A26V1D0udJgyVkkw4zOjx47FZXQMHD2HGxsDlsrUQT70XX0Mjvl07cfl0kPRujJ9O9FjXMVq7WjnaeZSTPSdjA9N5mXlsLdrK66pex7bibdQW1ZKbqYcXL0eaMFRKMNEoY2fOTBTLHThAdGgIgMz168l7x9vtrK719aTl5CQ52sWtf6yf413HYwniWNexWNeS2+VmU8Em3rr+rWwt2sq2om2UZ+vAtLI0YaikMMYQvHTJHsn0UjOBvXuJ9PUBkFFVRc4b34h/TxO+hgbSC3UOoLkaDY9yqucUrV2ttHa30trVGquYBlidu5p7V93L1qKtbC1yupZ0Sg11C5ow1IIJXbs2ceKg5r2Eb94EIL20lKwHH7SHuzY14l65dE9xmUihaIjzfedtcuhq5Xj3cc71niNs7GytJd4SaotqeXTdo2wp3EJtUS3ZGdlJjlotJpowVMKEu7omFcuFLju1EAUF9rSjTXvwNzbgrqrSLo9ZipoolwYucbzreCw5nOo5xVhkDICcjBxqi2p5uvZpaotqqS2qpcRXkuSo1WKnCUPNm0h//0QtxN5mxs6eA5xaiN27KXj3E/gam8jcsF4TxCyMnyXuRPcJjnfbBHGy5yTDoWHAnmd6U8Em3r7x7dQW2uRQkV2hf2M17zRhqDmLDg8TOHQoNlA9euIEGGNrIXbtIufND9tDXbUWYsaMMbQPtXOi+0QsQZzoPsFg0E6G6Ha5qSmo4U1r3hTrVlqTu4Y0V1qSI1fLgX6K1YxFx8YYaTkSG4MYOXoUwmFwu/HV1VH00Y/ib2zAW1eHaC3EHRljaB9s53jP8ViCONl9koHgAADprnQ25G/goeqH2FK4hS2FW1iXt04HpVXSaMJQtxRfCzHc/BIjhw5P1ELU1lL49NP4mpzzQnh1uunbiUQjtA22cbL7pF167OVgyLYc0l3prM9bz2urXsuWoi1sLtzM+rz1ZKRp4lWpQxOGijHRKGOnT0+cWe7AAaLDtp88s6aG/He+A994LUS2Hl1zK6FIiPP95znZfZJTPac42WMvxwvhMlwZbMjfwBtWv4GawhpNDmrRSOmEISJ/D7wJ6DDG1CY7nqXGGEPw4sWJYrm9e4n09wOQsXq1nY+psdHWQhQUJDna1DQcGuZ0z2lO9ZyKLef6zhGKhgA7IF1TUMOj6x5lc+FmNhVsYk3eGp1GQy1KKZ0wgC8Bfwl8OclxLBnB9qsTZ5Zrbibc2QlAetlKsl79antmuYYG3KU6cVw8YwwdgQ5O956elCAuD16ObZOfmU9NQQ3v3vxuNhVsoqaghqqcKlyis+OqpSGlE4Yx5kURqU52HItZqKODwN599sxyzXsJtbcDkFZYiL+xwZ5Zbs8e3OU6/cO4YCTIhf4LnO45zene05zpOcPp3tP0jfXFtinPKmdT4SYeXvswNQU11BTUUOIr0b+hWtJSOmHMhIg8AzwDUFlZmeRoki/S18fwvn2xYrng+fMAuHJy8DXspuA978G/p4mMdeuW/ZfbeKvhTO8Zzvad5UzvGU73nOZS/6VYdXRmWibr89bz6spXsyF/AzUFNWzI30BWRlaSo1dq4S36hGGMeRZ4FqC+vt4kOZwFFxkaZuTggdiZ5cZOnrK1EF4vvvp68t7yKL6mPXg21SBpy/dY/aHgEOf6znGm9wzn+s5xtvcsZ/vO0j/WH9um1F/KhvwNPFjxIBvzN7IhfwOVOZWkuxb9x0SpeaGfhEUmOjbGyOEWhptfItC8l5FjxyASQdxuvDt2UPSxX8ff2Ih369ZlWQsxGh7lQv8Fzved52zfWc71nuNc3zmuD1+PbeNL97Eufx2vrXotG/I3sD5vPevz1+uU3UrdgSaMFGdCIUZaW2Pnpx45dAgTDEJaGp7aLRR+4AP4mxrx7tiBy7N8Tn85Eh7hUv8lzvef53zfec71neN833naB9tjJ/lJd6WzOnc120u28/b8t7M2dy3r89dTllWmA9FKzUFKJwwR+SrwIFAkIu3A7xtjvpjcqBLLRKOMnToV62Ia2X+AaCAAQOamTeQ//ji+PU4tRNbS70fvH+vnYv9FLvZf5EL/Bbv0XeDq0NWJxCDpVOZUxqbMWJu3lvV566nIqdDDV5WaRymdMIwxjyc7hkQzxhC8cGGiFmLfvsm1EI88jL+xCV9jA+n5+UmONjEi0QjXhq5xceAil/ovTVz2X4yd2AdswVtVbhW1RbU8vO5h1uSuYV3eOiqzK3W6DKUWQEonjKUq2N5OoLk51oqIdHYBU2ohGhtxr1iR5EjnjzGGrpEu2gbauDx4mUsDl2jrb4vdHi90A8jNzKU6p5oHyh9gde5qVueuZk3uGlZlrdJJ9pRKIk0YCyB0s4PAvr2xVkTo6lUA0oqK8Dc04NvThL+padHXQkRNlI5AB1cGr3B54LK9HLwcux0IB2Lbul1uKrIrqMqp4oGKB6jOqaY6p5rVuavJ9yzNlpRSi50mjAQI9/YS2Lc/VlEdvHABAFduLv6GBgqeftrWQqxZs+gSxEBwgGtD17g6eJX2oXauDl2lfbCdK4NXuDp0dVJLId2VTnlWORXZFexasYvK7Eqqcqo7hdLOAAAgAElEQVSozKmkzF+mrQWlFhlNGPMgMjRE4MCBWLHc2ClbC+Hy+fDurifvscfwNTbY80K4UvfonKiJ0jPaw43hG1wfvs61oWuTLq8OXY2dl2FctjubVdmrWJ+/nldWvJLybJsgKrIrKPWXag2DUkuIfprnIDo6ysjhw7H5mEZaW20tREYG3h07KP74b+BrbMRbW4u4U2MwNhwN0z3STUegg45ABzcDN2OXN4ZvcGP4BjcDNye1EMDWLJRllVGWVUZdcR3lWeWUZZWxKnsV5VnlWrug1DKiCWMGTDDISGtrbAxi5PBhTCgEaWl4t22j8IMfwN+0B++O7bgyMxcuLmMYCA7QM9pD90g33aPdscuukS66RrroDHTSOdJJz2gPUROd9Ph0VzorfCtY4VvB1qKtvNb/Wkr9pZT6SynLKmOlfyU5GTmLrttMKZUYmjCmYSIRRk+empjV9eBBTCAAIng2bSL/3e+2xXK76knL8t/98xnDWGSMweAgg8FBBoIDDAQH6B/rp2+sL3bZN9ZH32gfPWM99I720jfaF5vzKJ5LXBR6CinyFlHkLWJz4WaKvEUUe4sp8ZVQ4i9hhW8FBZ4CLWBTSs3YnBOGiPiBUWNMZB7jSYrwUCej7R2M7DvAyN59jOzfT3TAORPa6mq8b34D6bt3IDtqieb4GYmGGYiGCI5dIhQIEYwEGYuMEYwEGY2MMhoejV0GwgFGwiMEQgEC4QCBUIDh0DDDoWGGQkMMh4YZCA4Qjr78i3+cIORm5pKXmUdeZh7lWeVsK9pGXmYeBZ4CCr2FsctCTyF5mXk6oKyUmnczThgi4gLeCTwB7AbGgEwR6QS+BzxrjDmbkCgT7Me/+kqq2mzeu5kHx6uFY9UujlcJfVntQDt0Pwc/nv2+BcGb7sWb7sXv9uN3+/G5fZT4SljtXk2WO4vsjGyyMrLIdmeTk5lDTkYO2RnZ5GXmkZuZS3ZGtrYElFJJN5sWxs+wX5mfBlqNsR3iIlIAvBL4rIh8yxjzlfkPM7EyX7WL090X6M/vIugZIhNoyMyhqWA1UriW9ML1pPuKSHOlke5KJ92VjtvlJsOVQUaaXTLTMslMyyQjLQNvuhdPmgdPuofMtEwdA1BKLQlizMxmBBcRtzEmdLfbJFJ9fb05cODA3HdgDHSfg4sv2uXSzyHgTE1RsAZWP2CX6gcgq3h+glZKqSQTkYPGmPo7bjfThHGHJys1xty46x3dpbtOGFNFo9BxwkkgL0DbL2FswN5XshlWv8ImkKp7wJs3f8+rlFILaKETxr8bY37lrnd0l+Y9YUwVCcP1IzZ5XHwBLu+F8AiIC1Zun2iBVDZBxt0fPaWUUgthQRNGqkh4wpgqPAbtB5wE8nNo3w/RELjcUL57IoGU10P6wtVnKKXUbGjCSIbgMFx+ySaPiy/C9RYwUUj32lbHeAJZuR3StARGKZUaZpowZv2tJSK/N916Y8wfznZfS06GH9a9xi4AI3123GN8DOQnf+Bslw3V904kkJItkMJzTCmlFMytcG847roHeBNwcn7CWWK8eVDzRrsADHXaI6/Gu7DOfN/ZrgBW3+8kkFdA4TrQQ3GVUinmrrukRCQT+I4x5qH5CWnukt4lNVv97U7rw0kiA/Y8GWSvnGh9rH4A8iqTG6dSaklLWJfUNHzAmnnYz/KTWw7b32UXY6DnwkQNyLmfwNF/ttvlV0+0PlY/AFklSQ1bKbU8zWUM4xgw3ixJA4oBHb+4WyJQuNYu9U/bBNJxcmL84/hzcOjLdtviTXFFhPeCV89Qp5RKvFl3SYlIVdzNMHDTmGmmTE2CRdclNRvRSFwNyIvQ9pKtAUFgZd1EC6SyCTKzkh2tUmoR0UrvpS4ctHUfl5xDeK/sc2pA0mFVfVwNyG5we5IdrVIqhWml93ITDMCVZps8LrwQVwPigYrGiRZI2Q6tAVFKTbIkCvdE5PXAX2DHSv7OGPPZ222/rBPGVKP9tgbkwgu2FXKz1a7PyLZzX423QFbUag2IUstcQo+SEpF8YD22DgMAY8yLc9nXbZ4jDfg88FqgHdgvIt8xxpyYz+dZsjy5sPENdgEY7pqYgffii3D2B3a91oAopWZoLkdJfQD4OFAOtABNwEvAq+Y3NBqAc8aYC87zfg14BNCEMRf+Iqh9q10A+q9OJI8LL8CJ5+x6rQFRSt3CXFoYH8eeca/ZGPNKEakB/mB+wwJgFXAl7nY70Dh1IxF5BngGoLJSv9xmLHcV1L3TLsZA78WJGpDzP5u+BqT6fshekdSwlVLJM5eEMWqMGRURRCTTGHNKRDbOe2QwXb/IywZcjDHPAs+CHcNIQBxLn4g9QVTBGtj1lE0gnacmEsiJuBqQoo1xNSD3ga8gqaErpRbOXBJGu4jkAd8GfiQivcC1+Q3LPg9QEXe7PEHPo6YSgZJNdmn8kK0BuXF0IoG0/BPs/1tAoHTrRAukag9kZic7eqVUgtzVUVIi8gogF/i+MSY4b1HZfacDZ4BXA1eB/cC7jDHHb/UYPUpqgYSDcO3QxPhH+z6IBEHSYNWuiRZIRaPWgCi1CMz7YbUiIuYOG89km9kQkTcCn8MeVvv3xpg/vt32mjCSJDQCV/ZOJJBrh8FEIC0TKhom5sBatRPS3MmOVik1RSISxvPAN4DnjDGX49ZnAPcB7wV+Zoz50lwCng+aMFLE6IA9kdSFF+DSi3DjmF2fkQWVe2zyWPMKWLFVa0CUSgGJSBge4H3AE8BqoA9bh5EG/BD4vDGmZc4RzwNNGClquNsewjt+GG/XGbvem28HzsdbIEUbtAZEqSRIaKW3iLiBImDEGNM3h/gSQhPGIjFwfeJEUhdehH6nwZq1YnINSH51UsNUarlYElODzJYmjEXIGOi9NDGN+8Wfw3CHvS+vcnINSM7KpIaq1FK1kCdQUmruRKBgtV12vdepATk9MY37ye/C4a/YbYs2xNWA3K81IEotMG1hqNQWqwFxxj/afgmhYWwNSO3E+EflHvDkJDtapRalRAx6FwM5wGVjTOgu40sITRjLQCQEVw9NdGFd2QeRMacGZOeUGhBvsqNValGYt4ThzNV0PzCCPTJqFdAP/HdjzMV5iHXeaMJYhkIjNmmMj39cPag1IErN0nwmjJ3GmENT1vmBMmPM2bsLc35pwlCMDdrT115ypjG5fhQw4PbbqUvGWyCl28CVluxolUoJepSUUgCBHmj7z4kq9K7Tdr0n1w6cjyeQ4hqtAVHLVsKOkhKRTxlj/mxuYSm1wHwFsOnNdgEYvOEMoD9vL0991673F0+pAVmtCUSpKWbSJfX1+JvAdmPM+oRGNUfawlCzFqsBcY7CGrph1+dWTE4gOWVJDVOpRJrPFsaAMeYDcTv+wpQn+gHwVmPMsHPbY4wZnW3ASiVFfrVddr7H1oB0nZmYxv3096DlH+12hesmFxH6C5MZtVJJMZMWRo0x5lTc7XxjTO94YhCR7wG9wAeNMQER+StjzK8lOO5paQtDzatoFG62ThzC2/ZLCA7Z+1aMnwfkfqi6x46JKLVIzWcL48MiEsWeS7sPWC0iq4F/BP4TuAT8D+BZEfkg058pT6nFx+WCldvscs+v2xqQay0T4x8HvgjNn7c1IGU7JteAZPiSHb1S825GR0mJSD6wHSgAzhhjjsXd9wVjzEdEZB323N4hY8xTCYr3trSFoRZUaNSePGp8/OPqAYiGIS0DyhsmEsiqXZCekexolbqlRFR6b5nubHciUmiM6Xau1wA/Mcasmm3A80EThkqqsUG4vNe2QC684JwHxIDbN3EekNUPwMo6rQFRKSURCeOQMWanc/0Dxpi/i7vPZ4wJONezjTGDc4z7rmjCUCnldjUgVfdNJJCSTXoIr0qqRNRhxL+jfw34u7jbPwd2ASQrWSiVcl5WA3LTngfkwvP28vS/2/X+4slFhAVrNIGolDSbhBHfFJn6btbzbCp1J9krYOtjdgHobbOtj/EzER7/pl2fUz65BiQ3KT28Sr3MbBJGqYg8BRzh5Qlj6cwvotRCya+C/Cdh55O2BqT73EQNyJnvw5F/stsVrJ18HpCs4uTGrZat2SSMzwD1wNNAuYgcB045S9H8h6bUMiICRevtsvv9tgak48REDUjrN+Dg/7HblmyZSCBV94A3L7mxq2VjzpMPikg5sA3YCmw1xrx7PgObCx30VktWJAzXjzg1IC/C5WYIj4K4YOX2iQRS2QQZ/mRHqxaZRT1brYi8Ddui2QQ0GGNmlAU0YahlIzxmzwMyPv7Rvt/WgLjcUL57IoGU10N6ZrKjVSlusSeMTUAU+BvgtzRhKHUHwWG4/NLEIbzXjwAG0r221TE+D9bKOkib9STVaolL2PTmC8EYcxJA9NBCpWYmww/rXmMXgJFeO/fVhRdsK+Qnf2DXZ+ZA1b1xNSCb7RQoSs1ASiaM2XBOIfsMQGVlZZKjUSpFePOh5lfsAjDUMfkQ3jP/Ydf7iuwEiuMtEK0BUbeRtC4pEfkxUDrNXb9rjHnO2eZ5tEtKqfnXd8UpInzBJpDBa3Z9zqopNSDlyY1TLYiU75IyxrwmWc+t1LKXVwHb32UXY6Dngq1Av/ginP0hHPmq3a5gTVwNyANaA7LMLfouKaXUXRKBwrV2eVkNyIvQ+k04+CW7bcnmuBqQe7UGZJlJ1aOk3gL8f0Ax9hwcLcaYh+70OO2SUioBYjUgL8TVgIw4NSB1cTUge7QGZJFa1IfVzpUmDKUWQHgM2g9MtEDa90M05NSA1MfVgOzWGpBFQhOGUmphBIfhyt6JQ3ivHQYTdWpAGuNqQLZrDUiKSvlBb6XUEpHhh7WvsgvAaD9c+s+JQ3h/8ofOdtlQHV8DskVrQBYZTRhKqfnlyYWaN9oFYKgT2n4xcQjvme/b9d6CyTUgheu0BiTFacJQSiVWVjFseYtdAPrbJ86DfvFFOPGcXZ+9cnINSJ4W4qYaHcNQSiXPeA3IePK4+CIEuux9+dUTrY/q++0JqFRC6BiGUir1xdeA1D9tE0jHSecQ3p/D8efg0JfttsU1cQnkXjv9iVpQ2sJQSqWuaMSpAXFOJHW5GUIBQGDltokEUrkHMrOSHe2ipYfVKqWWnnAQrh6caIG074NIEFzpsGpXXA1IA7g9yY520dCEoZRa+oIBuNI8MYh+7TCYCKR7oKLBSSAPQtkOrQG5DU0YSqnlZ7Qf2l6amMbkZqtdn5Ftz38+3gJZUas1IHF00Fsptfx4cmHj6+0CMNzlFBD+3CaRsz+w67359sir8TGQovVaAzIDmjCUUkuXv2hKDcjViQr0iy/Cye/Y9Vmlk2tA8quSF3MK04QBnOsYZHA0jIgg2B8aLhFcIqS57JLuXLrTXKSn2cvMdBfuNBdpLv1lotSikLsK6t5pF2Og9+LE+MeF5+HY1+12eVUTrY/V90P2dOd6W340YQC//53j/Oe57jk/Ps0leNJdZLrT8KS78GSk4UlPw5uRhs9Z/Bnp+DLtpT8znazMdLI86WRnppPjdZPtSSfH4ybH6ybHk056mvavKpVQIvYEUQVrYNd7bQLpPDXRfXXyO3D4H+y2RRudBHK/7cryFSQ39iTRQW+g5UofvcNBDPZvEY1C1BiiBiJRQ8QYolFDKBIl7FwGw1GCkSihsGEsHCEYjjIWjjISijDqLIFghJFQhMBYhOFgmEAwwtBYmGA4eseYsjPTyfW5yfO5yfdlkOfLIN+5XuC3S6E/g8KsTAqzMsj3ZWhLR6n5FI3AjWMTh/C2/RJCw4BA6daJFkjVHsjMTna0d0WPkkphoUiUodEwQ2NhBkZDDI6GGRwNMzASYmA0RP9IiL5AiIGREL2BIL2BEH3OZf9IaNp9ugQK/BkUZWVSnD2xlGR7WJEzcbkix4PHnbbAr1ipJSAScmpAnPGPK3ttDYikTa4BqWgAtzfZ0c6KJowlKhyJ0hsI0TMcpHt4zF4OBekaGqNraIzOwSCdQ2N0DY7ROThGMPLy1kyu101pjofSXA8rcz2szPXayzwPZXleynK9eDM0qSh1W6ERmzTGE8jVQ7YGJC3TqQF5hU0gq3ZCmjvZ0d6WJgyFMYa+QIiOwTFuDozGlhsDo9zot+uu94/SNTT2sscW+DMoy/NQluulPN9Heb7XWXyUF3jJ8aT2B0CpBTc6AJdfcgbQX4Cbx+x6t39yDUjpVnCl1g8yTRhqxoLhKDcHRrnWN8K1/hGu9Y1ytW+Eq70jscuRUGTSY3K9bioKvFTk+6gs8FFR4KOq0F5flefVQXulhrvtIbzjh/F2nbHrPXlQfd9EC6R4Y9JrQDRhqHljjKE3EKK9N0B77wjtvQEu9wS40jPCld4A7T0jk7q+0lzCqjwvVYU+qgv9VBf5qS70UV3kpyLfR0a6JhO1DA1cd5LHC3DhRei/bNdnrYgrInzATuu+wAlEE4ZaMNGo4ebgKJe7A7R1B2jrGaatO8Cl7mHaugIMjoVj27oEyvN9rC7ys7rIz5piP2uKslhd7GdljgeXHumlloveS5PPAzJ0067PrbSJY41zHpCclQkPRROGSgnjrZOLXcNc6hqmrXuYC13DXOoe5mLnMMPBia4urzstlkTWFmextiSLtU5C0UF4taQZY7usxqdxv/hzGO2z9xVtmGh9JKgGZFEnDBH578CbgSBwHnjaGNN3p8dpwlhcjDHcHBjjYtcwF7qGuNA5zIXOIc53DtPeGyDqvDVFYFWel3UlWawrzmL9iizWlWSzriSLXK8OvqslKBqxEyeOD6BffgmCQ9gakNqJ8Y/KPeDJueunW+wJ43XAT40xYRH5MwBjzKfu9DhNGEvHaCjCpe5hzncMc65jiHOdQ5zrGOJC5xBjcYWPK3IyWe8kjw0rstmwIov1K7I1kailJRKyh+2Oz8J7ZR9ExpwakJ02eTR8aM6nsV3UCSOeiLwFeMwY88SdttWEsfRFooarvSOc7RjkbIdNImdv2uuBuO6t0hwPG0qz2VCSxYbSbGpKs1lfkq1dW2ppCI3aGpDxI7CuHoTfPAlZJXPa3VJKGP8G/LMx5iu3uP8Z4BmAysrKXW1tbQsZnkoR0ajhap9NJGduDnHm5iBnbg5y9uZEi0QEKgt8bFxhE0jNyhxqSrOpKvTrtCpqcQsGIMM354enfMIQkR8D000B+bvGmOecbX4XqAfeamYQqLYw1FSRqOFyT4DTNwbtcnOAUzcGudQ1HBsj8bhdbBhPIqU51KzMZvPKHPJ8GckNXqkFkvIJ405E5L3Ah4FXG2MCM3mMJgw1U6OhCOc6hjh53SaQUzcGOHl9kJ7hYGyblbkeakqz2VyWw6aVOWxemaOtEbUkLeoz7onI64FPAa+YabJQajY87jRqV+VSuyo3ts4YQ+fQGKeuD3Ly+oCzDPLzs12EneaI150Wa4FsLrNJpKY0R8dG1LKQki0METkHZALjJ6loNsZ8+E6P0xaGSoSxcISzN21r5MT1AU5cs8lkYNQWJLoEVhf52VKWy+ayHGrLctlSlkO+X7u01OKwqFsYxph1yY5BqXGZ6dO3Rtp7Rzh+bSKJHLjUw3eOXIttU5brYcsqmzxqy+zjV+RkInruaLVIpWTCUCrViQgVzqSLr6+dOHajZzjIiWsDHL/Wz3Hn8scnbzLekC/KymBLWS61q3LYuiqXLWW5lOd7NYmoRUEThlLzqMCfwX3ri7hvfVFs3fBYmBPXBzh+tZ/WawO0Xu3nF+e6iDjjInk+d6wFstVZKgo0iajUowlDqQTzZ6azu7qA3dUTcwCNhiKcvjHIsav9tF7tp/VaP1/8xQVCEZtEcr1utjrdYNvKbRLRlohKNk0YSiWBx51GXUUedRV5sXVj4clJ5NjVfv7u5xdiR2jl+9xsLc9j26pctpbbRFKa49EkohaMJgylUkRmehrbyvPYVj6RROJbIsfa+zl6tZ8vvHA+1p1VnJ3JtlW5zuNsEinMykzWS1BLnCYMpVLYdC2R0VCEE9cHOHqlj6NX+zna3s9PT3fEBtZX5Xmpq7BJpK48j63luWRl6kdd3T19Fym1yHjcaeyszGdnZX5s3dBYmNar/Rxt7+NIu7383rEbgJ1Da11xlk085bnUVeRRU5qjZz5Us6YJQ6klICsznaY1hTStKYyt6x4asy2QK/0cae/jZ6c6+NeD7QBkpLvYvDKH7RV5bHdaMNWFPh0PUbeVkpXec6WV3krdmjF2Rt8jV/ppudLLkXY7LjISstPC53rdbCvPjSWR7RV5Oh6yTCzqSm+l1PwTEcrzfZTn+/iVbfY80eFIlHOdQxy50kfLlT4OX+7j8z87F5vJt6LAy/aK/FgC2VKWg8et82YtV9rCUEpNEgiGab06QMuVXlqu9NFyuY9r/aMAuNOETXFdWTsq87UrawlY9NObz4UmDKUSo2NglMNOK6Tlch9H2/sYds5wmOdzT0og28vzyPXpKXIXE+2SUkrNm5IcDw9tKeWhLXberEjUcLZjkJbLE11ZL5w5Gzu0d22xnx2V+eyozGNHRT4bVmSRnqZHZS122sJQSs2LwdEQR9v7OXy5l8OX+zh8pS92QipfRhrbynPZWZkfSyRFOqCeMrSFoZRaUNkeN/euK+LedXbiRWMMV3pGOHTZjoUcutzLsy9OTHVSWeBjZ6XtxtpZmU/Nymzc2gpJaZowlFIJISJUFvqoLPTx6I5VAIwEI7Re6+dQm22F/PJ8N99usecQ8bhdbFuVx46qvFhhYnG2tkJSiXZJKaWSxhjDtf5RDrX1csjpyjp+rT82a29FgZedlfnsqnJaIaXZOhaSANolpZRKeSLCqjwvq/K8vLmuDLBzZR2/1s+hNtuN9dL5bp5zWiHxYyE7K/PZWZVPgZ4Kd8FowlBKpRSPO41dVQXsqrLnDxmvUD90uS/WEokfC1lT5GdnlW2F7KrKZ11xFi6X1oUkgiYMpVRKi69Qf9hphYwEIxxt7+PQ5T4OtvVOmicr25POjsp8djldWdsr83S23nmif0Wl1KLjzUijcU0hjc5ki8YY2roDHGzr5UBbL4cv9/K5n5zBGHAJbCzNYVdVnm2FVBboKXDnSAe9lVJL0sBoiMNON9ZBJ4mMV6cXZ2fGWiC7qvPZUpZDZvrynSNLB72VUstajsfNKzYU84oNxYCtTj99Y5CDl3s5eKmHg5d7+f5xe86QjHQXdeW57KoqoL5KB9NvJSVbGCLyR8AjQBToAJ4yxly70+O0haGUmo2OgVEOOi2QA229kw7pXVPsZ1dlPvXV+eyqKmBtsX/JdmMt6skHRSTHGDPgXP8NYLMx5sN3epwmDKXU3RgNRTja3s+Bth4OXurl4OVe+gIhAPJ9budIrALqq/PZuip3yUz1vqi7pMaThcMPpF5WU0otOR53Gg2rC2hYbQ/pjUYNF7qGOHCpN9YS+fHJDgAy0lzUrsphd3VB7JDepX7CqZRsYQCIyB8D7wH6gVcaYzpvsd0zwDMAlZWVu9ra2hYuSKXUstM9NDapG+tYez/BSBSwNSH11fnUO62Q1UWLoxsr5bukROTHQOk0d/2uMea5uO0+DXiMMb9/p31ql5RSaqGNhiIcu9rvtEJ6ONA20Y1V6M9gZ1U+u51xkK2rcslIT72pTVI+YcyUiFQB/26Mqb3TtpowlFLJFt+Ntf9SLwfaemjrDgCQme6iriKP+qp8dlcXsLMqn1xv8k82tagThoisN8acda5/DHiFMeaxOz1OE4ZSKhV1DNoJFm0C6eX41X7CUYMIbFyRzS4ngexeXcCqPO+Cx7fYE8Y3gI3Yw2rbgA8bY67e6XGaMJRSi0EgGKblch8HnHGQQ229DI2FASjL9VBfXcDu6nzqqwvYsCKbtATPjbWoE8ZcacJQSi1Gkajh1I0Bpxurh/2Xerg5MAbYubFiLZDqAraVz//hvJowlFJqkTLG0N47woG2HvZdtIPpZ24OAfZw3q3ludRX57PbORorz3d3VemaMJRSagnpCwRtC6SthwOXejna3herSt+wIou/ebKe1UX+Oe17URfuKaWUmizPl8FrNq/gNZtXAPZw3iNX+th/yR7KW5rjSXgMmjCUUmoR8rgnT/G+EFKvgkQppVRK0oShlFJqRjRhKKWUmhFNGEoppWZEE4ZSSqkZ0YShlFJqRjRhKKWUmhFNGEoppWZkSU0NIiKd2Nlt56II6JrHcBYDfc3Lg77mpe9uX2+VMab4ThstqYRxN0TkwEzmUllK9DUvD/qal76Fer3aJaWUUmpGNGEopZSaEU0YE55NdgBJoK95edDXvPQtyOvVMQyllFIzoi0MpZRSM6IJQyml1IxowgBE5PUiclpEzonI/5vseBJNRP5eRDpEpDXZsSwEEakQkZ+JyEkROS4iH092TIkmIh4R2SciR5zX/AfJjmmhiEiaiBwWke8mO5aFICKXROSYiLSISELPUb3sxzBEJA04A7wWaAf2A48bY04kNbAEEpEHgCHgy8aY2mTHk2gishJYaYw5JCLZwEHg0SX+PxbAb4wZEhE38Avg48aY5iSHlnAi8ptAPZBjjHlTsuNJNBG5BNQbYxJeqKgtDGgAzhljLhhjgsDXgEeSHFNCGWNeBHqSHcdCMcZcN8Yccq4PAieBVcmNKrGMNeTcdDvLkv91KCLlwK8Af5fsWJYiTRj2i+NK3O12lviXyXImItXADmBvciNJPKdrpgXoAH5kjFnyrxn4HPDbQDTZgSwgA/xQRA6KyDOJfCJNGCDTrFvyv8SWIxHJAr4BfMIYM5DseBLNGBMxxmwHyoEGEVnS3Y8i8iagwxhzMNmxLLB7jTE7gTcAH3W6nBNCE4ZtUVTE3S4HriUpFpUgTj/+N4B/NMZ8M9nxLCRjTB/wPPD6JIeSaPcCDzt9+l8DXiUiX0luSIlnjLnmXHYA38J2syeEJgw7yL1eRFaLSAbwTuA7SY5JzSNnAPiLwEljzJ8nO56FICLFIpLnXPcCrwFOJTeqxDLGfNoYU26MqfWMZM8AAANqSURBVMZ+jn9qjHl3ksNKKBHxOwdyICJ+4HVAwo5+XPYJwxgTBn4d+AF2MPTrxpjjyY0qsUTkq8BLwEYRaReR9yc7pgS7F3gS+4uzxVnemOygEmwl8DMROYr9UfQjY8yyOMx0mVkB/EJEjgD7gH83xnw/UU+27A+rVUopNTPLvoWhlFJqZjRhKKWUmhFNGEoppWZEE4ZSSqkZ0YShlFIpaD4nCRWRV8YdIdgiIqMi8uhs96MJQ6k7EJGI8yE7IiKHROSeW2znFZEXnAkt5/P5i0UkYYdKqpT1Jeap2NIY8zNjzHan8v9VQAD44Wz3owlDqTsbcT5sdcCngT+9xXbvA75pjInM55MbYzqB6yJy73zuV6W26SYJFZG1IvJ9Z96on4tIzRx2/RjwH8aYwGwfqAlDqdnJAXpvcd8TwHNg560SkZ84LZJjIjLtDMjOuQyKnOv1IvL8Lfb9bWf/anl7FviYMWYX8FvAX81hH+8EvjqXJ0+fy4OUWma8zqyvHmwF9aumbuBMK7PGGHPJWTUKvMUYM+AkhGYR+Y6Ze6XsAeC/zfGxaglwJs+8B/gXO9sNAJnOfW8F/nCah101xjwUt4+VwFbszBazpglDqTsbcfp+EZE9wJdFpHbKl38R0Bd3W4A/cWYOjWKnzF8B3JhjDB1A2Rwfq5YGF9A3/l6M50yoOZNJNd8OfMsYE5prAEqpGTLGvIRNDsVT7hrBtkDGPeFss8v5gN+ccv+4MBOfw+nuJ+6+kbnErJYGZ0r+iyLyNrCTaopI3Sx38zhz7I4CTRhKzYozyJj2/7d3hyoRxEEcx7+DD+ATiM0uXhU0mgxGMXnBLCYfwmI2+wCi6Yq+gILYLwkGMYvIz7B7YPCOVW45w/eTdpdl9r9lh9llZ4DX78eTvAFLVTV56C/TzGb4qKptYHVKyDGw0W7vzbj0Gj12IdX/M6VJ6D5w2DYbfOIX00Hb4WErwO2f12TzQWm2qvoEHie7wGmS6x/OuwAuk4za7xZXNKNRH2g65u4kGVfVDTBM8lxVmzSt119opgAOkmxV1QA4SjJsY58A70nO+71baToThjQnVbUOHCc56CH2HbDbVjLSQvhKSpqTJPc0Myjm/uMecGay0KJZYUiSOrHCkCR1YsKQJHViwpAkdWLCkCR1YsKQJHXyBVZrew+IE3upAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the energy correction as a function of magnetic field amplitude\n",
    "for x in range(degS*degI):\n",
    "    plt.plot(np.array(Blist), Earray[:,x])\n",
    "plt.xlabel(\"B (a.u.)\");\n",
    "plt.ylabel(\"$E^{(1)}_Z$ (a.u.)\");\n",
    "plt.title(\"Zeeman effect in the hydrogen ground states\");\n",
    "#plt.savefig(\"Zeeman-H-au.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-02-24T15:07:15.230086Z",
     "start_time": "2019-02-24T15:07:15.000523Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8XFd98P/PkUYazYyW0S5Zq3dbXmTJtiSHBAKEJKShSYEEEsjGEqBQmteP/kop/bW0fehDf09LoU+fQtM+PA1LCUtpA5RCCSUbjeVFi21Z3i1ZsmRrX0fSbOf541zNjGTJkrXNSPq+X6/7mpl779w5V5qZ75x7zvccpbVGCCGEWEpx0S6AEEKItUeCixBCiCUnwUUIIcSSk+AihBBiyUlwEUIIseQkuAghhFhyElzEslNKbVdK1SulhpVSn1RKOZRSP1JKDSqlvrfMr92klLpziY51p1KqfYmO9aRS6rVYO9Z6ppRqUUrdFe1yrBUSXFaQUup9SqmRGRatlPrDaJdvGf0u8JLWOkVr/dfAu4FcIFNr/dBCD6qU+pxS6ps320drvUtr/dICj6+VUlsWVDixLizkB8d6eV9JcFlBWutvaa2TIxfgGeA68PdRLt5yKgGapj0+p7X2R6k864ZSyhbtMiyUUio+2mUQi6C1liVKC1ABDAN3RqxLA/430AlcBf4bEB+x/QNAM9AP/Awoidimgd8EzlvH/VNgM/A6MAR8F0i09k0Hfgx0W8f6MVAYcayXrOf/yjrWfwBZNzmX+4EGYAD4L2Cvtf4/gQAwDowA3wa8gM96/MF5nNcu4OdAHyYQ/z5w77TjNM5SrhbgLuv+56y/wdetc2oCDszyvFesv+eodfz3AHcC7cCngC7rf/RUxHPswF8AV6xyfhVwzHL8J4HXrP37gcvA261tDwHHp+3/KeBfrfuZwA+t/+kR6//02rT3wcet98Fla91twFFg0Lq9LWL/jdb5DgMvAv8L+GbE9hrrfzoANDL1/Xqr75Pftf5uHcCHrLJusbb9I/AV4CfW3/0uzOfh65j3aSvwB0BcxP8zspyl1vFs8ykb8Jh1zF7gs5HvlRnKfR9w2jrOVeB3ABcwBgSt98gIsAGownzmBqxz/RvCn7sb3lc3+/xY2z5tveYwcBZ4a7S/u+b1/RbtAqzXBXADF4FPT1v/r8DfWW/cHOvL4yPWtgeBC8BOwGZ90P4r4rka86WTivlCngB+AWyyPqSngSesfTOBdwFOIAX4HtaXl7X9Jat82wCH9fgLs5xLJebLthqIB56wPqj2iGN9KGL/6V8Ks56XVbZOzJdrkvW4eqbjzFK20BeGtf+49UURD/x34PBNnhv64rMe3wn4gT8BEqzjeIB0a/uXrL9/hlXOHwH/fZZjP4kJjB+2yvIxzBeuwgSpPmBnxP71wLus+89jgqQL2I354pkeXH5ulcNh3fZjvkxtwCPW40xr/9cxQS4RuB0TtL5pbSvAfPneh7nS8TbrcfYC3if3Atcw700n8A1uDC6DwBus10rCBJYXrL9nKXCO8A+SKf9/Zg4uM5YNKMN8ub/R+nt/0frfzhZcOoE7rPvpQGXEe6J92r77MQHZZpWpGXjmJu+rWT8/wHagDdgQcY6bo/39Na/vuGgXYD0umC+QF6xFRazPxQQER8S6R4BfWvf/ffKDZT2Ow3y5lViPNfCGiO3HiQhewF8CX5qlTPuA/ojHLwF/EPH4N4GfzvLcrwB/Om3dWeBNEce6WXCZ9bys86+f5XWnHGeWfVqYGlxejNhWBozd5LkzBZcxrC8va12X9UWiML9GN0dsO4RVc5jh2E8CFyIeO63Xy4v4m37eur8LEwzs1pePD9gR8dw/48bg8paIx48BR6a9/utWGYoxX6rOiG3fJBxcPg18Y9pzf0b4R8qtvE++RkSwBbZwY3D5esT2eMznoSxi3Ucw7XczvY9KuTG4zFg24A+B5yO2uTA14dmCyxXrtVOnrb+TacFlhuc+A/zLTd5Xs35+rL9RF6YWl3Cz14m1RdpcouPTmF+cT2jrnWQpwfwi7lRKDSilBjC1mJyI7V+O2NaH+VIriDjG9Yj7YzM8TgZQSjmVUn+nlGpVSg1hquvuade5r0Xc90w+dwYlwKcmy2WVrQhziWA+bnZeRZhfn0tl+jkl3WK7RK+e2lY0+XfJxgSI4xHn8VNr/Zxl0Vp7rLuTf+PngEeVUgoTHL6rtZ6wjmfD/Jqd1DrDsSO3b5hhn1bM33cD0Bfx+tOfWwI8NO1/ezuQP9N5cPP3yYZpx26bYZ/IdVmY2lRk2SfLPV+zlW1KWbTWo5ga2Wzeham9tSqlXlZKHZptR6XUNqXUj5VS16zP1p9Z5zKbWT8/WusLmOD0OaBLKfW8Umq+n6uokuCywqxusZ8F3q21Hpi2uQ3zSy1La+22llSt9a6I7R+J2ObWWju01v+1gKJ8ClPlrtZap2IuD4D5Ur9VbZhf2ZHlcmqtv30Lz5/tvNow7UYz0bOsj4YeTPDeFXEOadp02rhlWuvDmF/SdwCPYi4hgWl78GO+fCYVz3SIiPsdmC+wSMWYy2mdQIZSyhmxLfLYbZiaS+T/xqW1/sKtnpP1WoWzvM5M5e7B1NIiyz5ZbjA1xchy591iWUKvb51/5mw7a62Paq0fwPzQ+1fMZcnp5Z30FeAMsNX6bP0+N/9c3fTzo7X+J6317Zi/gwb+fL4nGU0SXFaQUiofc738Ga11/fTtWutOTKPjXyqlUpVScUqpzUqpN1m7fBX4jFJql3W8NKXUQrvypmC+DAeUUhnAHy3wOGB6un1UKVWtDJdS6teUUinzfP7NzuvHQJ5S6hmllF0plaKUqra2XQdKlVLL9T6+jmmvmpPWOoj5O/yVUioHQClVoJS6ZxGv/3VMY7Bfa/2a9ToB4AfA56zaZxnmGv3N/ATYppR6VCllU0q9B3NJ8Mda61bgmHW8ROsX+TsinvtN4B1KqXuUUvFKqSSr+23hjS8zp+8CTymldlpf5jftfm+d63eBz1v/9xLg/7HKBKYB/I1KqWKlVBrwmVsoy/eB+5VStyulEjHtaDO+j6y/y/uUUmlaax+mTSpgbb4OZFqvPynF2mdEKbUD054Wafr7atbPjzI5Ym9RStkx7YVjEa8d0yS4rKwPY9pVvjxDrstXrX0ex1wKOI25zv59rEsQWut/wfxqed6qbp8C3r7AsnwJ08jZAxzGXMJZEK31Mcy5/Y1V5guY6/nzff6s56W1HsY0Ir8Dc4njPPBm66mTCZi9Sqm6hZb/Jj4HPGddqnh4Hvt/GnPuh63zeBFTO1yob2Aun35j2vpPYC7vXMO0U/yfmx1Ea92L6Y30Kcyln98F7tda91i7vA/TPtSL6Z34HUwNGq11G/AA5td3N+ZX9v/LAr47tNb/Dvw18EvM3+l1a9PETZ72W5gayiVM77p/wrTdoLX+uVXWE5j2xR/fQlmaMD3q/glTi+nH9ASczWNAi/V//Sjwfus4ZzA9IC9Z75MNmJ5kj2J6d/29VcZInyPifTXH58cOfAHzOb2GqTn9/nzPM5rU1Ev+QohYoZRyYBpzK7XW51fwdb8DnNFaL6Y2O5/X2Yn5IWHXkvO05kjNRYjY9THg6HIHFqXUQevya5xS6l5MTeVfl+m1fsO6zJSOqa3+SALL2rRqs3eFWMuUUi2YRuAHV+Dl8jDtOJmYS0Mfm6lNcIl8BHMpLwC8jOkeLNYguSwmhBBiycllMSGEEEtu3V4Wy8rK0qWlpdEuhhBCrCrHjx/v0VrfLDkYWMfBpbS0lGPHjkW7GEIIsaoopWYaEeIGcllMCCHEkpPgIoQQYslJcBFCCLHkJLgIIYRYchJchBBCLLmoBRelVJFS6pdKqWalVJNS6ret9RlKqZ8rpc5bt+nWeqWU+mul1AWl1AmlVGXEsZ6w9j+vlJprhFghhBDLLJo1Fz/wKa31TsxMfh+3hg//PeAXWuutmCl6f8/a/+3AVmt5GjNnAhHDxVdj5q7+o8mAJIQQIjqiFly01p1a6zrr/jBmnukCzKB5z1m7PUd4bKUHMFOgamsiJbc1P8o9wM+11n1a637M3OH3ruCpCCFEzPP39TH005/S9cW/WpHXi4kkSqVUKVAB1AK51qRZaK07JydewgSeyClQ2611s62f6XWextR6KC6eafI+IYRYGwLDw3iOHsNTe5jRw7VMnD0LQJzLRcaTT2DLyFjW1496cFFKJQP/jJmdcUipWWcDnWmDvsn6G1dq/SzwLMCBAwdkxE4hxJoRHBvDU1eH53Ato7W1jJ86BcEgym7HUVlB9jPP4DpUQ9KuXSjb8n/1RzW4KKUSMIHlW1rrH1irryul8q1aSz5msiQwNZLIObcLMXODtwN3Tlv/0nKWWwghok17vYydOMHo4Vo8hw/jaWwEnw9sNhx795L10Y/grK7Bsa+cOLt9xcsXteCiTBXlfwPNWusvRmz6IWZO8C9Yty9ErP+EUup5TOP9oBWAfgb8WUQj/t3c2lzaQggR83QgwPjp04wePozncC2eujr02BgoRVJZGRmPP4aruhrn/v3EuVzRLm5Uay5vwMxLfVIp1WCt+31MUPmuUuqDwBXgIWvbT4D7MPNLe4CnALTWfUqpPwWOWvv9ida6b2VOQQghlocOBpk4fyHUZuI5epTg8DAAiVs2437Xu3DVVOM8eJD4tLQol/ZG63aysAMHDmgZFVkIESu01vhaWxk9XMto7WE8tUcI9JnfyQnFxbiqq3DW1OCqqsKWPeeI98tGKXVca31grv2i3qAvhBDrla+zM9RmMlpbi//aNQBsOTkk33E7zqpqXDXVJBTM2AE2pklwEUKIFeLv7cVz5IipnRx+HV/rFQDi3W5TK6mpxlldTWJpKTfpObsqSHARQohlEhgawnPsWKgRfuLcOQDikpNxHjhA+iOP4Kqpwb5tGypubQ31KMFFCCGWSNDjwVNXH2qEH29qMrkmSUk4KytJvf9+XDXVJJWVrUiuSTSt7bMTQohlFPR6GW9sDDXCjzWeMLkmCQkm1+RjHzPBpLycuMTEaBd3RUlwEUKIedJ+v5VrYiUu1tWhx8chLo6kXbvIfPIJnNU1OCsriHM6o13cqJLgIoQQszC5JudNb67JXJOREQDsW7fifvghXDU1OA8cID41NcqljS0SXIQQwqK1xtvSgqe21gST2loC/f0AJJQUk3rffbgO1eCsqsKWmRnl0sY2CS5CiHXN19FhBZLDjNYeCeea5OaS/MY7cNYcwlVdRcKGDVEu6eoiwUUIsa74e3oYra0NjR7su2LlmqSn46ypxlVdjaumhoSSklWfaxJNElyEEGtaYHAQz9GjodrJxPkLgJVrcvAgGe9/H87qGuxbt6y5XJNokuAihFhTgqOjeOrqQomL46dPg9Ym12T/flLf8etmXpOdO9d8rkk0yV9WCLGqBb1exhoazGWuw4cZO3EC/H5ISMBZXk7Wxz+Oq7pqXeaaRJMEFyHEqqL9fsabmkKXuTzH69ATE+Fck6eewllTjbNCck2iSYKLECKm6WCQiXPnwpNkHT1KcHQUAPu2bbjf8zCumkM4Dx4gPiUlyqUVkyS4CCFiitYa7+WW8CRZtbUEBgYASCwtNeNzVVfhrK6WXJMYJsFFCBF1vqtXw5NkHa7F39UFgC0/n+Q3vzk0FH1CXl6USyrmS4KLEGLF+bu7Ga09Eqqd+NraAIjPzAzPuFhTQ0JRkeSarFISXIQQyy4wMMDokSOhxEXvxYsAxKWm4qw6SMbjj+OsrsK+dasEkzVCgosQYskFRkYZqzseGj14vLnZ5Jo4HDj378f9zt/AWVVNUtlOVHx8tIsrloEEFyHEogUnJhirbwi1mYydPAl+PyohAce+fSbX5FANjj17UJJrEhXegJdTPac413+O9+5477K/ngQXIcQt0z4fY6dOhUYPHqurQ3u9Jtdkz24yP/ABXDXVOCoqiHM4ol3cdWlwYpDG7kbqrtdR31XPqZ5TeINeAO4pvYf0pPRlfX0JLkKIOelgkIkzZ8IzLh49RtDjAcC+YwfpjzxiEhcPSK5JtHSOdFLXZQJJXVcdF/ovoNHYlI2dmTt5ZMcjVOZWsi9n37IHFpDgIoSYgdYa76VL4cTFI0cIDA4CkLhxI6kP/DquatM92Ja+/F9UYqpAMMCFgQs0dDVQ11VHXVcd10bNVAFOm5N9Ofu4u+Ru9ufuZ3fWbhy2la89SnARQgDgbW8Pz7hYW4u/uxsA24Z8kt/61nCuSW5ulEu6/oz7xznVcypUK2nsamTYNwxAtiObytxKntz1JBU5FWxL34YtLvpf7dEvgRAiKnzXu/AcqQ3VTnxXrwIQn5WFq6oK5yEr16SwULoHr7C+8T4auhpCweR072n8QT8AW9xbuGfjPVTmVFKRU0FBckFM/n8kuAixTvj7+/EcPTpjromruoqMp57CVVNN4ubNMflltVZprbkyfIW663U0dDdQd72OlqEWABLiEtidtZvHyx6nIqeCipwK0uxp0S3wPElwEWKNCoyMMnb8WKgRfqL5jMk1cTqtXJN34qypJmnHDsk1WUG+gI/Tfaepv15PfVc9Dd0N9I33AZCamMq+nH08uOVBKnMrKcsswx5vj3KJF0aCixBrRHB8nLGGhtBlrrGTJyEQQCUmmlyT3/oErppDOPbsRiUkRLu468Zkl+D6rvpQl+CJwAQARSlF3F5wO/ty9lGZU8nGtI3EqbUxG6YEFyFWKe3zMXbyVGh8rrH6epNrEh+PY88eMj/0IZO4uG8fcUlJ0S7uuqC1pm24LRRIGroauDhoLj9Odgl+ePvDVOaYLsFZjqwol3j5SHARYpXQgQDjzWfCQ9EfP472eEAp7Dt3kP6+9+GsrsJ54CDxya5oF3dd8Aa8nO49HWp8j7zElZKYQnl2Ofdtuo+KnIqodQmOFgkuQsQorTXeixdDMy6OHjlKcDLXZNMm3A8+aBIXDx6UXJMV0jfeR2NXI/Xd9TR2NU7Jep+8xFWeXU5FTgWb3ZvXzCWuhZDgIkSM0Frja28PtZmM1tYS6OkBIKGwkJS73mpmXKyuIiEnJ8qlXfuCOsjlwcuhWkljd2OoF5ctzkZZZhmP7HiEfTn71vwlroWQ4CJEFPmuXw+Nz+U5fBhfRwcA8dlZuA4dCs1tklhYGOWSrn0en4dTPado6G6goauBxu5GhrxDAKTb0ynPKefBLQ9SkVPBrqxdq7YX10qR4CLECvL39+OpPRIaPdh7+TIAcWlpuKqqyPjgB3BVS67JctNa0zHaEQoiDV0NnOs/R0AHANiUtom3lbzN1Eqy91GSWiL/j1skwUWIZRQYGZmSuDhx5gwAcU4njgP7cT/0EK6aauw7dqDi1u/1+eXmDXhp7mueEky6x8zwNg6bg71Ze/nA7g9QkVPB3uy9qyZRMZZJcBFiCQXHxhirrw8lLo6fagrnmlRWkv3bn8RZU4Njt+SaLKcuTxeN3Y00djXS2N3I6d7ToYb3guQCDuYdZF/OPsqzy2NmLK61Rv6iQiyC9noZO3kynLjY0ID2+cBmM7kmT38YV7U1r4ldrtEvB1/Qx9m+syaYWAGlY9S0XSXGJbIraxeP7nyUfdn7KM8pl4b3FSLBRYhboAMBxk83T801GRsL55o89piZJKtyv+SaLJOesZ4ptZKm3qZQxnuOM4d92ft43873UZ5Tzs6MnSTGy8yX0RDV4KKU+hpwP9Cltd5trcsAvgOUAi3Aw1rrfmVa074M3Ad4gCe11nXWc54A/sA67H/TWj+3kuch1i6tNd4LF0KXuTxHjhIcMj2IEjdvDo3P5Tp4kHi3O8qlXXt8AR9n+8O1khPdJ7g6YkZvtsXZKMso46FtD1GeU86+7H3kufKiXGIxKdo1l38E/gb4esS63wN+obX+glLq96zHnwbeDmy1lmrgK0C1FYz+CDgAaOC4UuqHWuv+FTsLsWZorfG1tU3NNentBSChqIiUu99mck2qDkquyTK4NnqNE90nzNJzgtO9p6fUSsqzy3lkxyOUZ5ezM3OndAeOYVENLlrrV5RSpdNWPwDcad1/DngJE1weAL6utdbAYaWUWymVb+37c611H4BS6ufAvcC3l7n4Yo3wXbsWyjUZrT2Mv6MTAFtODq433GbNuFhDYmFBlEu6toz7x2nua+ZE94lQzaTL0wWYtpKyzDIe3v4w5dnllGeXS61klYl2zWUmuVrrTgCtdadSavLnYQHQFrFfu7VutvU3UEo9DTwNUFxcvMTFFquFv69vSuKit7UVgHi3G2d1Na4PfxhndTWJGzdKbsMSmZyzJLJWcq7vHH5tJsAqSC5gf+5+yrPL2Zu1lx0ZO0iIl950q1ksBpfZzPQp1zdZf+NKrZ8FngU4cODAjPuItScwNITn2LHQpa6Jc+cAiHO5cB44gPuR9+KqqcG+bZvkmiyRwYlBTvac5GT3SU70nOBkz0kGJ8y4aE6bk91Zu3li1xOUZ5ezJ3uP9OBag2IxuFxXSuVbtZZ8oMta3w4URexXCHRY6++ctv6lFSiniFFBjwdPXX2oR9d4UxMEgyi7Hef+SlLvewbXoRqSdu1C2WLxI7C6TDa6TwaTkz0nQ2NwKRSb3Zt5a/Fb2ZO1hz1Ze9ji3kJ8nExOttbF4ifrh8ATwBes2xci1n9CKfU8pkF/0ApAPwP+TCk1OSzs3cBnVrjMIoq018vYiROhy1yexkaYzDUpLyfrox/FWV2No2IfcYnSLXUxJucrOdlzMrSc6T0TSlDMcmSxJ2sPD2x5gD1Ze9iVuYvkxOQol1pEQ7S7In8bU+vIUkq1Y3p9fQH4rlLqg8AV4CFr959guiFfwHRFfgpAa92nlPpT4Ki1359MNu6LtcnkmpwOXeby1NWFck2SysrIfOJxnNU1OCsriHNJrsli9I710tTbxInuE5zqOcWp3lOhy1tJ8UmUZZbx6M5H2Z21m71Ze8lz5Uk7lQBAmc5X68+BAwf0sWPHol0MMQ86GGTi/IVw4uLRowSHhwGwb92Cs6oa16EanAcPEp8mY0ItlMfnoam3iaaeJk72nKSptymUUxKn4tji3sLurN2hy1ub3Ztl2JR1SCl1XGt9YK795J0hYo7WGl9razhxsfYIgT5TGU0oKSb17W/HWV2Fq7oaW5Y0BC+EN+DlbN9ZmnqbONVziqbeJi4NXiKog4DpvbU7azeP7HiEXZm7KMssw5ngjHKpxWoiwUXEBF9nZ6jNZLS2Fv+1awDYcnNJvuN2nNU1uKqrSCiQXJNb5Q/6uThwkdO9p0OXts71n8MfNN2AM5Iy2J21m7tL7mZX1i52Z+0mIykjyqUWq50EFxEV/t5ePEeOMPr6YUZrD+NrvQJAfHo6zqoqXB/9iMk1KS2Va/i3IBAM0DrUai5vWZe4zvSdYTwwDkByQjK7MnfxeNnj7Mo0gSTflS9/Y7HkJLiIFTEl1+T1w0ycPw9AXHIyzoMHSX/kEVyHDmHfulVyTeYpqIO0DrVyuvc0Tb1NnO49TXNvMx6/BzAN7jszd/Lube+mLLOM3Vm7KUktWdfzuouVI8FFLItZc02SknBWVpL6jnfgqq6SXJN5CuogLUMtnO49HVrO9J1h1DcKgD3ezo6MHTyw5YFQG8nGtI3S4C6iRt55YkkEvV7GGxtDjfBjjSdMrklCAo7yvWR97GO4aqpJKi+XXJM5+IN+Lg1eorm3mea+Zpp7mznTdyZUI7HH29mevp13bHoHZZlllGWWSc8tEXPk3SgWRPv9jDc1mUb42sN46urR4+MQF0fSrl1kPvkEzpoanBUVxDmll9FsJgITnO8/T3NfM2d6z9Dc18y5/nOhkYAdNgfb07fz4JYH2Zm5k7LMMjalbZJAImKevEPFvOhgkIlz58KJi0ePEhy1Lsls24b74YfM6MEHDxKfmhrl0samwYlBzvWfC9VEmvuauTx4mYAOAJCSkMKOzB08vP1hUyPJKKMktUSGShGrkgQXMSOtNd7LLeHExdpaAgMDgJVrcv/9uGqqcVZVYcvMjHJpY4vWmo7RDs72neVs31nO9J3hTN+Z0NS7ADmOHLZnbOfNRW9mZ+ZOdmTsoDC5UHptiTVDgosI8V29Gk5cPFyLv8uMGWrLyyP5zjvNjIvV1STk50e5pLFj3D/OhYELnOs/Z4JJ/1nO9Z9j2GtGEFAoSlJLKM8u5+HtD7MjYwfbM7bLKMBizZPgso75e3oYrbUSFw/X4msz0+LEZ2SYDPiaQ7hqqkkoLl73v6gnayPn+89zrv9caGkdag1ltTtsDralb+O+jfexLX0bOzJ2sMW9RTLbxbokwWUdCQwO4jl6NNQIP3H+AgBxKSk4q6rIeOz9OGtqTK7JOg4mgxODnO8/z/mB85zvP8+FgQuc7z/PiG8ktE9hciFb07dyT+k9bEvfxvb07RSmFEoOiRAWCS5rWHB0FE9dXagRfvz0adDa5Jrs30/qr/86rpoaksrKUPHrr9HY4/NwceAiFwYuhJf+C3SNdYX2SUlMYat7K/dvup+t6VvZlr6NrelbcSXIaMtC3IwElzUkODHBWENjqBF+7MQJ8PvDuSYf/ziummoce/ei1lGuyahvlMuDl7kwcIFLA5e4MHCBiwMXpzSw2+PtbErbRM2GGra4t7DFvYWt6VvJdeau61qcEAslwWUV034/46dOhRMX6+rRExMm12T3bjKfetLkmlRWEudwRLu4y65/vJ9Lg5fMMnCJy4OXuTh4kWuj10L7JMQlUJpWSnl2Oe/c+k4TSNK3UJhcKF1+hVhCElxWER0MMnH2bHjGxWPHwrkm27eT/t734JzMNUlJiXJpl4cv6KN9uJ2WwRZahsxyefAylwcvMzAxENovKT6JjWkb2Z+7n81pm9mUtonN7s0UphRKAqIQK0A+ZTHsZrkmiSUlpL7jflw1NSbXJGPtDJEe1EG6PF20DrXesLQPt+PX/tC+GUkZlKaWclfJXWxM3cjGtI1scm8i35UvjetCRJEElxgza65Jfj7Jd95pZlysriYhLy/KJV2cQDBA52gnbcNttA230T7cTutQK1eGr9A+3B4aIh5Me0hxajFb3Fu4q+QuSlNLKU0rpTS1lDS7zDwpRCyS4BJl/u5uRmuPhGqaUjQ2AAAgAElEQVQnoVyTrCwznEp1Fa6aGhKKilZVw7LWmsGJQa6OXuXq8FXaR9pDt+3D7XSMdEypgSTGJVKUUkRRahG3bbiNktQSilOLKUkpIdeVK7UQIVYZCS4rLDAwwOiRI3gO1zJaW4v34kUA4lJTcVYdJOOJJ3BVV5G4ZUtMB5NAMEDveC+do51mGemkY6SDztFOOkY76BjpCA0HPynNnkZBcgE7M3dyd+ndFCYXUpRSRHFqMTnOHAkgQqwhElyWWWBklLHjx0KXuiaaz5hcE6cTZ2Ul7nf+Bs7qGpJ27oiZXBNf0EfvWC/XPdfp8nTR5eni+uh1rnmumdvRa3R5uqbUPMAMvJifnE9BcgEHcw9SkFxAQUqBuU0uICVxbXYyEELcSILLEguOjzPW0BBKXBw7eRICAVRCAo6KCrI+8XFchw7h2L17RXNNgjrI0MQQfeN99I730jvWS+94Lz1jPfSM9dA91k2Px9z2jffd8PyEuARynbnkunKpyK0gz5lHniuPDckbyHPlke/Kl+AhhAiR4LJI2udj7OSpcOJifT3a64X4eBy7d5P5oQ/hqq7CUVlJXFLS4l9Pa8b8Y4z4RhiaGGLIa5bBiUEGJgZCt5NL/3g/feN9DE4MhoZ2j2RTNjIcGWQ7sslz5bEnew85jhyynFnkOnPJceaQ48wh3Z4e05fphBCxRYLLLfIPdzHe2omn9ihjh2sZO3YcPTYGQML2rTjf8y7iD+4jrnwXfmcio0E/A0EfvpGLeIe8eANeJgITeANexgPjjPvHmQhMMOYfw+PzMOYfY8w/xqhvFI/fw6hvlFHfKCPeEUZ8I4x4R264HBUpTsXhtrtDS0lqCfty9pFuTyfTkUlGUgYZSRlkJmWS6cgkzZ4mbR1CiCUnweUW/eqeN5FjXTVqz4SmHYpTpXGcLlYMOy8Dl6Hze9B568eOU3E4bA4cNgeuBBdOmxNXgos8Vx7J7mRcCS5SElNISUwhOSGZVHsqqYlmSbOnkWZPIzkhWYKFECLqJLjcIt/byjjjvc6guxtfoodk4DZHOrdnbEJlbMaWtQ1bkpv4uHgS4hKwxdlIiEsgMT6RxLhEEuITsMfbQ0uSLQmHzUGSLYnEuES59CSEWBMkuNyiu/74n82dYBC6m+HyK2a5+Bqc/qXZlr0TNr7RLKW3gSM9egUWQogoUFrraJchKg4cOKCPHTu2dAcM+OFaowk0l16GK4fBPwYoyC+3gs2boLgG7MlL97pCCLGClFLHtdYH5txPgssy8Xvh6jETaFpehbYjEPRBnA0KDoRrNoUHIWHxvciEEGIlSHCZw7IHl+m8Hrjyugk0l1+BjnrQQbAlQVF1uGazoQLi5WqlECI2LVlwUUodAt4P3AHkA2PAKeDfgG9qrQcXX9yVt+LBZbqxARNsLr8Cl1+F6yfN+sQUKLktXLPJ3Q1x0vtLCBEb5htcbvoTWSn170AH8ALweaALSAK2AW8GXlBKfVFr/cPFF3mdcbhh+9vNAjDaAy2vWcHmZTj/M2u/dCi9I1yzydoK0qNMCBHjblpzUUplaa17bnqAeewTi6Jec5nLUIep0Vx+2QScQTNaMsl54VrNxjdCekl0yymEWFeWtM1FKfX/Af+otW6LWPe01vrZxRUzemI+uETSGvovh7s9X34FRrvNNndJuFaz8Q5IWd3zvAghYttSB5cuoAf4uNb6l9a6Oq115aJLGiWrKrhMpzV0NYc7B7S8CuNW01fW9ogcm9vBuXZmqBRCRN+StLlEuAo8AHxPKfV9rfX/AOTCf7QoBbllZqn+CAQDcO1EOMem4Vtw9O8BBXl7wjWbkkNgl5GLhRDLb741l3qtdYVSKgn4CpAM7NFa71juAi6XVV1zmYvfCx11ETk2tRDwgoqHgv3hmk1RteTYCCFuyVLXXI4BaK3HgaeUUh8H9i+ifGI52RLNSADFNcCnwTdmRgyYvIz22l/Bq38B8XYonp5jkxDt0gsh1gBJolyPxofCOTaXXo7IsUmelmOzR3JshBBTLFWey0lg1uijtd67gLItC6XUvcCXgXjgH7TWX4hykWJXUipsu8csAKO94VrN5Vfg/H+Y9Y500ylg45tMsMnaJjk2Qoh5meuy2P3WrcJk5N+3vMVZGKVUPPC/gLcB7cBRpdQPtdano1uyVcKVCbseNAtYOTavhPNsmn9k1ifnTsuxKY1akYUQse2mwUVr3Tp5Xyk1Efk4xlQBF7TWlwCUUs9jerdJcFmI1A1Q/l6zhHJsXg1fRjv5PbOfuzjcXlN6B6TmR7fcQoiYsVZGSCwA2iIetwPV03dSSj0NPA1QXFy8MiVb7ZSCjE1m2f+ECTbdZ81ltEsvQfOPof6bZt+sbRE5NndIjo0Q69hcbS6RSZIOpVQFEfktWuu65SrYLZqpIeCGtiJrRIFnwTToL3eh1iSlIGeHWao+bCZNu34y3O258Xk4+g+YHJvd4faa4kOmrUcIsS7MVXP5y4j714AvRjzWwFuWvEQL0w4URTwuxAy4KZZbXJyZDC2/HN7wSQj44GodtFidA478Pbz+N1aOTeW0HBtHtEsvhFgma6IrslLKBpwD3ooZTeAo8KjWumm256zrrsgryTcG7UfD7TVXj4MOQHxixDw2b4QNlSY/RwgR05ZkbDGl1Dtv9mSt9Q8WULZloZS6D/gSpivy17TWn7/Z/hJcomRi2CR0XnrJBJxrVm/3BJcZnmYy2OTthbj4aJdWCDHNUgWXINBgLTC1bUNrrT+wqFJGkQSXGOHps3JsrN5oPWfN+qS0iHls3gjZOyTHRogYsFTDv7wLeA+wFzNh2Le11heWoHxCGM4MKHvALADD16bOY3Pmx2a9K8dMKTDZ9Tm9VIKNEDFsvgNXujB5I+8BMoHPaq1fXuayLSupuawS/S3hWs3ll2HkulmfVjQ1oTN1Q1SLKcR6sdQDV44Dg8AQUIyZ6liI5ZdeapbKx0yOTc95q1bzMpz9iZleACBzS7jbc+kdZtQBIUTUzNXm8mbgEUwG/IvA81rrNfFzX2oua0AwCNdPWcHmVWj9FXhHzLbcPeFaTcltkmMjxBJZygb9E8BrmLyWKTtrrT+5yHJGjQSXNSjgg4768ACcVw5DYMLk2GyomJpjk+iMdmmFWJWWKrg8yc1HRX5uQaWLARJc1gHfOLQfCefYdNRB0G9ybAqrwsGmYL/k2AgxT0sSXNYyCS7rUGSOTcur0HkCk2PjNMPTTAab/HLJsRFiFks1n8uzwF9rrU/NsM2F6T02obX+1oJLKsRKsafA1reZBWCsH1peC/dGe/GPrP3SoPQN4WCTUybdnoW4RXP1Fvtb4A+VUnuAU0A3pqfYViAV+BoggUWsTo502PkOswAMX4+YNM3qjQbgzLJybKzeaBmbJNgIMYf55rkkAweAfGAMaNZan13msi0ruSwm5jRwZWpC53CnWZ9aODXHJq0guuUUYgVJm8scJLiIW6I19F6Eyy+FZ+kc6zPbMjZPnccmOTuqRRViOUlwmYMEF7Eokzk2k+Oitf4KJobMtpxdU3NsHO7ollWIJSTBZQ4SXMSSCvihs8FcQrv0MrTVgn8cVBzk7wsHm+JDkmMjVrWlynP5htb6MaXUb2utv7ykJYwyCS5iWfknzDw2kzN0th81OTZxCVB4MBxsCg9Kjo1YVZYquJwG3g78ELiTadMJa637FlfM6JHgIlbUxIjJsZmcobOjAdBgc0ybx6Yc4uc75J8QK2+pBq78KvBTYBNwnGnzuVjrhRBzsSfD1rvMAibHpvW/wqMHvPg5a79UKLFybDa9CbJ3mqmkhVhl5tsV+Sta64+tQHlWjNRcREwZ6TKXzyYvo/VdMutDOTbWPDaSYyOiTBr05yDBRcS0gbbwAJyXX4HhDrM+tWBajk1hdMsp1h0JLnOQ4CJWjVCOzcvhEQQ8vWZbxqaIHJs3So6NWHYSXOYgwUWsWsEgdJ0OB5qW1yJybMrCyZylbzBD3AixhCS4zEGCi1gzAn641hjuHHDlMPjHAGVGeJ5srymuMR0LhFgECS5zkOAi1iz/BFw9Hu4c0HYEgj6Is5m8mtI7TE+0woNgs0e7tGKVkeAyBwkuYt3weuDK6+HLaB31oINgSzK1mcmaTf4+ybERc1qqPBchxGqX6IQtbzULwPggtPwq3PX5F39i7ZcybR6bXZJjIxZMgosQ601SGuy4zywAI91T57E591Oz3pExNccmc4vk2Ih5k+AixHqXnA2732kWgMH28Oycl1+G0y+Y9Sn5U3Ns3MXRK7OIedLmIoSYndZmtIDIhE5Pj9mWXhqu1ZTeASm5US2qWBnS5iKEWDylIHOzWQ48ZYJNV3O4VtP0AtR93eybvSMiofN2ybFZ56TmIoRYuGAAOhutqaBfNb3SfB5Mjs3eiBybQ5Jjs0ZIV+Q5SHARYhn4vSbHZrJm034UAl6TY1Ow3wSayXlsEpKiXVqxABJc5iDBRYgV4PWYWTknazYddeEcm6LqcM1mQ4Xk2KwSElzmIMFFiCgYH4TW18PB5vpJsz4xBUpuC7fZ5O6WHJsYJQ36QojYk5QG2+81C8Bob0SOzStw/mdmvSPd9ECbDDZZ2yTHZpWR4CKEiB5XJux60CwAQx1Wjs3LJtg0/9CsT86dmmOTXhq1Iov5keAihIgdqRug/D1m0Rr6W8KdAy69DCe/Z/ZzF1uB5k4zikBKXjRLLWYgwUUIEZuUgoyNZtn/hAk23WfDwab5R1D/TbNv1vapOTbOjOiWXUiDvhBilQoG4NqJcHtN6+vgGwUU5O0J90QrOQT2lGiXds2Q3mJzkOAixBoT8EXk2LxiukAHvKDirRwbq2ZTVAUJjmiXdtWK6eCilHoI+BywE6jSWh+L2PYZ4INAAPik1vpn1vp7gS8D8cA/aK2/YK3fCDwPZAB1wGNaa+9cZZDgIsQa5xuzcmysGTo76kEHIN5uAsxkQmdBJcQnRLu0q0asB5edQBD4O+B3JoOLUqoM+DZQBWwAXgS2WU87B7wNaAeOAo9orU8rpb4L/EBr/bxS6qtAo9b6K3OVQYKLEOvM+JAZnubSy9DyClyzcmwSXFNzbPL2QFx8dMsaw2I6z0Vr3Qygbuy3/gDwvNZ6ArislLqACTQAF7TWl6znPQ88oJRqBt4CPGrt8xymRjRncBFCrDNJqbDtHrOAybFpfS1cs7nwc2s/t+kUsPFNZjpoybFZkFjrLVYAHI543G6tA2ibtr4ayAQGtNb+Gfa/gVLqaeBpgOJimYtCiHXNlQllD5gFYKgzPDvn5VfgzI/NesmxWZBlCy5KqReBmTqff1Zr/cJsT5thnQZmGgdC32T/GWmtnwWeBXNZbLb9hBDrUGo+7H3YLBDOsZkMNjfk2Fjz2KTmR63IsWzZgovW+q4FPK0dKIp4XAh0WPdnWt8DuJVSNqv2Ern/smjuHGLcF0AphQLilEIpcxsfZxabdZtoi8MWp0iwxZEYb5a4OKleC7EqpJeapfLxqTk2La9A848jcmy2ReTY3CE5NpZYuyz2Q+CflFJfxDTobwWOYGooW62eYVeB9wKPaq21UuqXwLsxPcaeAGarFS2J3/p2PRe6Rhb8fFucIikhHrstjqSEeJISzK0jIR6n3YYzIR6nPR5Xog2X3UayPZ5ku43kpASS7TZSk2ykOhJISbKR5kggJSmBeAlYQiwvpSBnh1mqn4Zg0Ay6eellcymt8Xk4+g+YHJvd4VpNyW2mrWcdilZvsd8A/ieQDQwADVrre6xtnwU+APiBZ7TW/26tvw/4EqYr8te01p+31m8i3BW5Hni/1SHgphbaW6z2Ui8ebwCNRmvzgyaoNUGtCQQhoDWBYBBfQOMPaPzBIF5/EG/A3E74J28DjHmDjPsDjHsDeLwBPL4AY14/oxMBPNatNxCc428JKXYbbmcibmcCbmci6c4E0p2JpDsTyXAlkOGyk+FKJCs5kaxkO2mOBKlBCbGUAj7T1XlymJq2IxCYMDk2GyoicmyqIdEZ7dIuSkx3RY4Fq6Ur8oQ/wOhEgJFxP0PjPoat26ExH0PjfgbHzP0Bj5d+T/i23+NleNw/4zHj4xSZLhNoslPMkmPd5qYmkZtqJycliZxUO3abdMkU4pb5xqH9SDih8+pxCPohPjFiHps3woZKsCVGu7S3RILLHFZLcFkMXyBI/6iX3lEvfdZt78gEPSMT9Ax76RmZoHtkgu5hs/iDN74XMlyJ5KYmkZ+WRF5aEhvSkshLc7AhLYkNbgd5aUkkJUgAEuKmJobN8DQtVrDpPAFoK8fmUHh6gfzymM+xkeAyh/UQXG5FMKjp93jpGp7g+tC4tUxwbWica4PWMjRO3+iNgx9kJdspcCdRkO6gwO2gMN1JYXr41mWPtaY9IaLM0wetvwr3Rus5a9YnpU2dxyZ7R8zl2EhwmYMEl4UZ9wW4NjhOx+AYHQPjdA6McXVy6Te3E/6p7UQZrkSK0h0UZjgpznBSYt0WZzrJT3NIhwQhhq9By2tw6SXTQaC/xax3ZU/LsdkY9WAjwWUOElyWh9aa7pEJrvaP0d4/Rlu/x9z2ebjS5+Fq/9iUy2+J8XEUpjsoyXRSkuliY5aLkkwnG7NcFLgd2OJlqluxDvW3hmfovPQyjFwz69OKpgab1A0rXjQJLnOQ4BId/kCQa0PjXOn10NrnobXXQ2vvKC3WrccbCO1ri1MUZ5hAszHLxcZsc7s5O5mcFPtMwwcJsfZoDT3nw7NztrwKY/1mW+aWqTk2rqxlL44ElzlIcIk9k7Welh4PLb2jtPSMcjliibzclmy3WYHGBJvNOclszk6mNMspPdzE2hYMwvVTVrB51bTdeK3cu9zd4WBTcptpw1liElzmIMFldQkGNZ1D41zuHuVSzwiXuke52G1urw6MhfaLU1Cc4WRLjgk4W3NS2JqTzJacZOlYINamgA86GuDySybYtNWCfxxU3LQcm5olybGR4DIHCS5rx5g3wMXuEbN0jXChe4QLXSNc7hnFFwi/vwvcDrbkJLMt1wo6uclszU0hWYKOWEt849B+NCLH5lg4x6bwoAk0Nb+54JEDJLjMQYLL2ucPBGnt83ChywSb89eHOW/dj7zEVuB2sD0vhW25KWzPS2Zbbgqbs5Mlf0esDRMjZh6byWDTewF+9xLY7As6nASXOUhwWb8CQU1bn4dz14etZYRz14e52D0SqunExylKM53syEtle14KO/JS2JGXSmG6Q4bOEaub17Ooy2MxPVmYENEUH6cozXJRmuXi7l3hWSF8gSCtvaOcuTbMWWs51THIT051MvkbzJUYb4JNfio7rdsdeSmkJMk0uWKVWKGxzaTmIsQcRif8nLs+zJlrw5zpHOLMtWGaO4cYihi7rTDdwc78VHbmp1JmLUUZDukuLdYcqbkIsURcdhsVxelUFKeH1mmt6Rwc58y1IZo7hzndOcSZziF+0XydyRzRFLvNCjgplG1IpSw/ja250pYj1gcJLkIsgFKKDW4HG9wO3rIjN7R+zBvg7HVTs2nuHKKpY4jvH29n9HWTHGqLU2zJSTa1mw2p7C5Io2xDKqlyWU2sMRJchFhCjsR49hW52VfkDq0LBjVX+jw0dQxxunOQpo4hXr3Qww/qr4b2Kc5wsrsglV0b0thlBZ2s5IX15hEiFkhwEWKZxUV0IPi1veH51ruGx03A6RiiqcMEnZ+cvBbanpeaFAo4ewrS2F2QRm6qDHsjVgcJLkJESU5KEjnbk3jz9pzQusExXyjYnLo6yKmOIX5xpivUWy0r2c7ugtRQsNlTkEZ+WpIEHBFzJLgIEUPSHAkc2pzJoc2ZoXWjE37OXBviZPsgJ6+awPPq+R4CVs+BTFciewrT2GsFnL2FbqnhiKiT4CJEjHPZbewvyWB/SUZo3Zg3wOnOIU5dHeREu6nlvHKuO9RTLTvFzt6CNBN0Ck3AkTYcsZIkuAixCjkS49lfks7+knD3aI/XT3OnqeGcsILOf54NX1LbkJZkBRs35YVu9hSmkeaQXmpieUhwEWKNcCbeWMMZnfDT1DHEifYBTrQPcqJ9gJ81XQ9tL810Ul7ktgJOGrs2pOFIlDwcsXgSXIRYw1x2G1UbM6jaGA44gx4fJ66aYNPYNkDtpT5eaOgAzNA423JT2FeURnmhm/IiN1tzkmVGUHHLZPgXIQRdQ+M0tFkBp32AxraB0PA2joR49hSkUV6URnmRuaRWmC5D26xXMiryHCS4CDE7rTUtvR5OtA9Qf2WAE+0DnOoYwmtNVZDpSqTcShbdV2RqONJ+sz7I2GJCiAVTSrExy8XGLBcP7CsAwOsPcvbaMA1WzaahbYBfRnQY2JTlMsGm2AScHXmpJNrkctp6JTUXIcSCDY37ONk+SEObqeE0tA3QMzIBQKItjj0FaaHazb4iuZy2FshlsTlIcBFi6Wmt6Rgcp+HKAPVX+mloG+Dk1cHQzJ9ZyXb2FbmpKHZTUeRmb5FbppleZeSymBBixSmlKHA7KHA7QuOo+QLmclr9lX7q2wZouDLAi82mO3Scgm25KdaUBm4qi91sykqW2T7XAKm5CCFW3IDHG7qUVnelf0rvtJQkG/uK3FRaAaeiKJ00p3QWiBVScxFCxCy3M5E7t+dwpzVoZzCoudQzSp11Ka2utZ//+Z/nQ8PZbM52UVGcTmVxOpUlbrbmpBAvtZuYJjUXIURMGpnwc6JtgHor2NS3DdA36gUg2W6jvChtSu0m3ZUY5RKvD1JzEUKsasl2G7dtyeK2LVmA6SzQ2uuh7kp/6HLa3750MTQ69KZsF5XFZry1yuJ0tuZI2000Sc1FCLFqebx+GtsGrYDTT92VcO0mxW5jX7HbupRmajgynfTiSc1FCLHmORNtU+a/mRxZoK61n+NX+qe03SgF23JSqCxJp7LYzf6SdDZmuSTvZplIzUUIsaYNj/tCtZvjraaGM9kzLcOVSGWxm8qSdPYXp7O30C2jQs9Bai5CCAGkJCVw+9Ysbt9q2m6CQc3F7hGOt5pgc/xKPy82dwFgi1Ps2pBqgo215Kc5oln8VUtqLkKIda9v1Etda3+odtPYPsC4z4wqsCEtif2lGewvdnOgNIMdeSnregoCqbkIIcQ8ZbgSuassl7vKcgEzqsDpjqFQ7ebo5T5+1GjmvHEmxrOvyM2BknT2l2ZIR4FZSM1FCCHmMDlm2rGWPupa+znW2k9z51Coo8D2XNNR4EBJOgdKMijKWLsDdMb0wJVKqf8BvAPwAheBp7TWA9a2zwAfBALAJ7XWP7PW3wt8GYgH/kFr/QVr/UbgeSADqAMe01p75yqDBBchxGKMTvhpaBvgWEs/x1r7qL8ywMiE6SiQnWI3NZuSdA6UZrBrQyoJa+RSWqwHl7uB/9Ra+5VSfw6gtf60UqoM+DZQBWwAXgS2WU87B7wNaAeOAo9orU8rpb4L/EBr/bxS6qtAo9b6K3OVQYKLEGIpBYKac9eHOdbaz/GWPo619tPePwZAUkIc5YVuDpSamk1lSfqqnVwtpttctNb/EfHwMPBu6/4DwPNa6wngslLqAibQAFzQWl8CUEo9DzyglGoG3gI8au3zHPA5YM7gIoQQSyk+TrEzP5Wd+ak8VlMCwPWh8VDN5nhrP199+RKB4MVQzs2B0vRQwFlrc93EQoP+B4DvWPcLMMFmUru1DqBt2vpqIBMY0Fr7Z9j/Bkqpp4GnAYqLixddcCGEuJnc1CR+bW9+aPoBj9dPw5UBjlntNj9s6OBbtVesfe0cKMkIBZud+au7V9qyBRel1ItA3gybPqu1fsHa57OAH/jW5NNm2F8DM/2F9U32n5HW+lngWTCXxWYtvBBCLANn4tTx0gJBzdlrwxxv7eNoSz/HWvr4t5Od1r7xVBS7OVCSwUGrV5prFU2stmwl1VrfdbPtSqkngPuBt+pww087UBSxWyHQYd2faX0P4FZK2azaS+T+QggR0+LjFGUbUinbkMpjh0oB6BgYMzWblj6OtfTz1/95Hq2tffNT2V+SzsHSDA6WppOTmhTdE7iJaDXo3wt8EXiT1ro7Yv0u4J8IN+j/AtiKqaGcA94KXMU06D+qtW5SSn0P+OeIBv0TWuu/nasM0qAvhFgNhsZ9Zqy01n6OtfRT39YfSvAsyXRaNRvTK21z9vKPlRbrvcUuAHag11p1WGv9UWvbZzHtMH7gGa31v1vr7wO+hOmK/DWt9eet9ZsId0WuB95vdQi4KQkuQojVyBcI0tQxxLGWPo5cNh0Feq2RoDNciVbNxgSb3RvSSLQtbbtNTAeXWCDBRQixFmitudwzylHrMtqx1n4u94wCpgv0viI3B0szOFCaQWWxm5RFjiYgwWUOElyEEGtV9/AEx1qsTgKtfTR1DBEIauIU7MxP5ZsfrF7wzJ0xnecihBBi+WSn2Hn7nnzevsd0gR6ZMF2gj7T0caZzCLdz+RM4JbgIIcQal2y3TZl2YCWs3gwdIYQQMUuCixBCiCUnwUUIIcSSk+AihBBiyUlwEUIIseQkuAghhFhyElyEEEIsOQkuQgghlty6Hf5FKdUNtC7w6VmY4f7XCjmf2LWWzgXkfGLdfM6nRGudPdeB1m1wWQyl1LH5jK2zWsj5xK61dC4g5xPrlvJ85LKYEEKIJSfBRQghxJKT4LIwz0a7AEtMzid2raVzATmfWLdk5yNtLkIIIZac1FyEEEIsOQkuQgghlpwEl2mUUvcqpc4qpS4opX5vhu12pdR3rO21SqnSiG2fsdafVUrds5LlnslCz0Up9Tal1HGl1Enr9i0rXfaZLOZ/Y20vVkqNKKV+Z6XKfDOLfK/tVUq9rpRqsv5PSStZ9pks4v2WoJR6zjqPZqXUZ1a67DOZx/m8USlVp5TyK6XePW3bE0qp89byxMqVemYLPRel1L6I99kJpdR75v2iWmtZrAWIBy4Cm4BEoBEom7bPbwJfte6/F/iOdb/M2t8ObLSOE79Kz6UC2GDd3w1cXSug+xsAAAT2SURBVM3/m4jt/wx8D/id1Xw+mBlkTwDl1uPMaL7XluB8HgWet+47gRagdBWcTymwF/g68O6I9RnAJes23bqfvkrPZRuw1bq/AegE3PN5Xam5TFUFXNBaX9Jae4HngQem7fMA8Jx1//vAW5VSylr/vNZ6Qmt9GbhgHS9aFnwuWut6rXWHtb4JSFJK2Vek1LNbzP8GpdSDmA950wqVdy6LOZ+7gRNa60YArXWv1jqwQuWezWLORwMupZQNcABeYGhlij2rOc9Ha92itT4BBKc99x7g51rrPq11P/Bz4N6VKPQsFnwuWutzWuvz1v0OoAuYMzsf5LLYdAVAW8TjdmvdjPtorf3AIOaX43yeu5IWcy6R3gXUa60nlqmc87Xg81FKuYBPA3+8AuWcr8X8f7YBWin1M+tSxu+uQHnnspjz+T4wivlVfAX4C61133IXeA6L+Tyvxu+COSmlqjA1n4vz2d92qy+wxqkZ1k3vqz3bPvN57kpazLmYjUrtAv4c80s52hZzPn8M/JXWesSqyMSCxZyPDbgdOAh4gF8opY5rrX+xtEW8JYs5nyoggLnskg68qpR6UWt9aWmLeEsW83lejd8FNz+AUvnAN4AntNbTa2ozkprLVO1AUcTjQqBjtn2sanwa0DfP566kxZwLSqlC4F+Ax7XW8/qlsswWcz7VwP+vlGoBngF+Xyn1ieUu8BwW+157WWvdo7X2AD8BKpe9xDe3mPN5FPip1tqnte4CfgVEe7yuxXyeV+N3wayUUqnAvwF/oLU+PO9XjVYjUywumF+ElzAN8pMNX7um7fNxpjZKfte6v4upDfqXiG6D/mLOxW3t/65o/0+W4nym7fM5YqNBfzH/n3SgDtP4bQNeBH5tFZ/Pp4H/w/9t745ZowiiAI7/H4GQoAiKEWwtxEI0FhZqIZZa6AdQEfMJRLCKhYggWPgZgrWIgoVIRBREEMETjKAEYicIURtJITIWM4EzyZGNjps9+P/gYNl7O3uPW+7dzN7N5G/YW4D3wIGu59MXO8PqG/oL5X3aXrZ3DGkuo8AT4NKGz7uZb2AXH8Ap4CN5XHG67LsOnC7bY+RfHM0Dr4A9fcdOl+M+ACeHNRfgKnkMvNf32DWs+axo4xodKC4VrrVz5B8nvANubXYu/3i9bS3758iF5cpm59Iwn8PkXsEPYBGY6zt2quQ5D1wc1lzKdfZzxWfBZJNzOv2LJKk677lIkqqzuEiSqrO4SJKqs7hIkqqzuEiSqrO4SC2IiF8R0YuIt2XKlqMD4sYj4llEHCzxvYj4GhELZXs2IiYi4lHbOUgb4fQvUjuWUkqTAJGXY7gJHF8jbgq4l/KklMvxM8DDlNLd5aCI+BwRx1JKL/77K5f+gj0XqX3bgG8DnjsLPGjQxv0SK3WSPRepHeMR0SP/S303sGoBtogYJf9r/VOD9l4DN6q+Qqkii4vUjv5hsSPAnYjYn/6cImMn8L1he1/IswhLneSwmNSylNJLciFZuejSErln08RYiZc6yeIitSwi9pGXnl3s35/yqoUjEdGkwOwlT1opdZLDYlI7lu+5QJ5a/kJae2nix+SFwGbXae8EeY0NqZOcFVnqkIg4BFxOKZ1fJ+45cKb0dqTOcVhM6pCU0hvgaUSMDIqJiAngtoVFXWbPRZJUnT0XSVJ1FhdJUnUWF0lSdRYXSVJ1FhdJUnW/AZjHrSNur9moAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for x in range(degS*degI):\n",
    "    plt.plot(np.array(Blist)*B0, Earray[:,x]*EH/h/1e6) # convert units to T and MHz\n",
    "plt.xlabel(\"B (T)\");\n",
    "plt.ylabel(\"f (MHz)\");\n",
    "plt.title(\"Zeeman effect in the hydrogen ground states\");\n",
    "#plt.savefig(\"Zeeman-H-SI.pdf\")"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}