{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0-π Qubit\n",
    "\n",
    "In this notebook, we try to reproduce the result of the [\"Coherence properties of the 0-π qubit\"](https://iopscience-iop-org.stanford.idm.oclc.org/article/10.1088/1367-2630/aab7cd) paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal is to reproduce the enegry spectrum and eigenfunctions that [Groszkowski2018](https://iopscience-iop-org.stanford.idm.oclc.org/article/10.1088/1367-2630/aab7cd) is calculated for the $0-\\pi$ qubit. The diagram of the circuit is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"pics/zeroPiQubit_pic1.png\" width=\"220\" align=\"left\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose the parameter set 3 for the circuit from the following table of the paper:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"pics/zeroPiQubit_pic2.png\" width=\"400\" align=\"left\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Circuit Description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, we import the SQcircuit and the relavant libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import SQcircuit as sq\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the single inductive loop of the circuit via `Loop` class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "loop1 = sq.Loop()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The elements of the circuit can be defined via `Capacitor`, `Inductor`, and `Junction` classes in SQcircuit, and to define the circuit, we use the `Circuit` class. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# define the circuit ’s elements\n",
    "C = sq.Capacitor(0.15, \"GHz\")\n",
    "CJ = sq.Capacitor(10, \"GHz\")\n",
    "JJ = sq.Junction(5, \"GHz\", loops=[loop1])\n",
    "L = sq.Inductor(0.13, \"GHz\", loops=[loop1])\n",
    "\n",
    "# define the circuit\n",
    "elements = {(0, 1): [CJ, JJ],\n",
    "            (0, 2): [L],\n",
    "            (0, 3): [C],\n",
    "            (1, 2): [C],\n",
    "            (1, 3): [L],\n",
    "            (2, 3): [CJ, JJ]}\n",
    "\n",
    "zrpi = sq.Circuit(elements)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By creating an object of `Circuit` class, SQcircuit systematically finds the correct set of transformations and basis to make the circuit ready to be diagonalized.\n",
    "\n",
    "Before setting the truncation numbers for each mode and diagonalizing the Hamiltonian, we can gain more insight into our circuit by calling the `description()` method. This prints out the Hamiltonian and a listing of the modes, whether they are harmonic or charge modes, and the frequency for each harmonic in GHz (the default unit). Moreover, it shows the prefactors in the Josephson junction part of the Hamiltonian $\\tilde{\\textbf{w}}_k$, which helps find the modes decoupled from the nonlinearity of the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\hat{H} =~\\omega_1\\hat a^\\dagger_1\\hat a_1~+~E_{C_{22}}(\\hat{n}_2-n_{g_{2}})^2~~-~E_{J_{1}}\\cos(-\\hat{\\varphi}_1+\\hat{\\varphi}_2+0.5\\varphi_{\\text{ext}_{1}})~-~E_{J_{2}}\\cos(-\\hat{\\varphi}_1-\\hat{\\varphi}_2+0.5\\varphi_{\\text{ext}_{1}})$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$------------------------------------------------------------$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\text{mode}~1:~~~~~~~~~~~\\text{harmonic}~~~~~~~~~~~\\hat{\\varphi}_1~=~\\varphi_{zp_{1}}(\\hat a_1+\\hat a^\\dagger_1)~~~~~~~~~~~\\omega_1/2\\pi~=~3.22489~~~~~~~~~~~\\varphi_{zp_{1}}~=~2.49e+00$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\text{mode}~2:~~~~~~~~~~~\\text{charge}~~~~~~~~~~~~~~~~n_{g_{2}}~=~0$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$------------------------------------------------------------$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\text{parameters}:~~~~~~~~~~~E_{C_{22}}~=~0.296~~~~~~~~~~~E_{J_{1}}~=~5.0~~~~~~~~~~~E_{J_{2}}~=~5.0~~~~~~~~~~~$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\text{loops}:~~~~~~~~~~~~~~~~~~~~\\varphi_{\\text{ext}_{1}}/2\\pi~=~0.0~~~~~~~~~~~$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "zrpi.description()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above output shows that the first mode of the `zrPi` circuit is a harmonic mode with the natural frequency of 3.22 in the frequency unit of SQcircuit (which is GHz in default), and the second mode is a charge mode."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To determine the size of the Hilbert space, we specify the truncation number for each circuit mode via `set_trunc_nums()` method. Note that this is a necessary step before diagonalizing the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "zrpi.set_trunc_nums([35, 6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Circuit Spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generate the spectrum of the circuit, firstly, we need to change and sweep the external flux of `loop1` by the `set_flux()` method. Then, we need to find the eigenfrequencies of the circuit that correspond to that external flux via `diag()` method. The following lines of code find the `spec` a 2D NumPy array so that each column of it contains the eigenfrequencies with respect to its external flux. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# external flux for sweeping over\n",
    "phi = np.linspace(0,1,100)\n",
    "\n",
    "# spectrum of the circuit\n",
    "n_eig=5\n",
    "spec = np.zeros((n_eig, len(phi)))\n",
    "\n",
    "for i in range(len(phi)):\n",
    "    # set the external flux for the loop\n",
    "    loop1.set_flux(phi[i])\n",
    "    \n",
    "    # diagonlize the circuit\n",
    "    spec[:, i], _ = zrpi.diag(n_eig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAELCAYAAADDZxFQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA88ElEQVR4nO3dd3yb933g8c8PgwRBAuAeIqlB7S1Z07Jka9myZMlu2vTipGmaNM100qa967j02lyTXi659tomzXDdNG1GEyeXKXlbsvYe1t6boihuEsRev/vjB1K0DImkRBIc3/fr9byAB3gIfB8SxPf5baW1RgghhLiTJd0BCCGEGJwkQQghhEhJEoQQQoiUJEEIIYRISRKEEEKIlGzpDqCvFBYW6rFjx6Y7DCGEGFIOHz7cqLUuSvXcsEkQY8eO5dChQ+kOQwghhhSl1LW7PSdVTEIIIVKSBCGEECIlSRBCCCFSkgQhhBAiJUkQQgghUpIEIYQQIiVJEEIIIVIaNuMgxPATjyXwNgZpqw/ibwsT8kcJ+WPEwnG01mgNSoEtw4o904rdYSUrJ4Mslx2nOwNXvgNHjh2lVLpPZVCIRuL4mkP42yIE280W8seIReLEwnHisQRYFEopLFaFw2kjM9tOlsuOp9CJpzgLR7Y93achBpAkCDEoaK1prvVTe7GNuitt1F3x0loX4M7lSmx2C7ZMK8qiUAp0QhONJIhF4pBiaROb3YKrwIGnKAtPsZPcEid5pU7yR2WTlZMxMCc3gOLRBC11AVpq/bTUBWitC9BWH6C9OUSwPZryZ6w2C7ZMC1abBbT5W8Rjmkgw9q5js1x2ise4KRnnpnSch7IJHmwZ1v4+LZEmkiBE2sSjCa6faebqiUaun2zC1xIGzJdQyTgP4x8qxlOcRW6xk5y8TBzZ9rt+GWmtiYbjBNujBNsjBNoitLeEaG8O0d4Uoq0+wI2zLcSiic6fyXLZKSjPobAih8JKF0WVLnJLnVgsQ6PEEfBGaLjeTuONdhpv+Gi64aO1PohOJDOlAleeA09xFuPmFOHKd+DKd5DtySDLnYHTlUGm04bFmrqmOZHQhANRAt4IbfWmJNd8y0/9VS/XTjWBBqvdQvmkPMbMKKBqThE5eZkD+BsQ/U0NlxXl5s+fr2WqjcEvkdDUnG3h/KE6Lr/dQCQYw+6wUjk1nzHTCyifnIe70NEv1UI6ofG1hmmp9dNc66fppp+mGz6ab/pN9Qpgy7RSVJlD8Vg3JcnNVdA/8fRGJBSj/lo79Ve91F3xUn/N25lQAVwFDgorcigozyG/LJv8Udl4irOw2fvn6j4SjFF7uY3rp5q4fqqZ1roAKBg1IZeJ84uZMK8ER45URw0FSqnDWuv5KZ+TBCEGgrcxyJm9tZzdU4uvJUyGw0rVnCImLCihYnKeqd5Ik0Q8QcutAA3V7dRfa6fhmpeG677OpOF0Z1A63kNplYey8R6KRrv6NV6tNe1NIWovtXHrUhu1l9torvF1Vre5i7IoGeumaLSL4tEuCitzyHSm98u4tS7AhUN1XDhYR8utABabompOEVOXlFE5JR81REplI5EkCJEWiYTm+qkmTm6vMVUSwOip+UxZUsa42YX9dnXbF+KxBE01PuqueLl1xXxRextDgKlWKRnrpmy8h7IJuZRWuR/oCzqR0DTd8FF7qZXai23UXmzF3xYBwO6wUjrOTUmVh9JxHkrGugf1lbnWmqYaH2d213LuwC3C/hjuQgfTHy1n2pJRgzr2kUoShBhQ4UCUM3tqObHtBt7GEE53BtOWjmLa0lG48h3pDu++BbwR8yV+qY3aC600VPtMfb+CglHZlFaZUkbxGDeekiysKer2tdYEvBEab/i4ddkknrqrXqKhOAA5eZmUTchNJh8P+aNyhkybyJ1i0ThXjjZyckcNNy+0YrVZmLigmFkrKymqdKU7PJEkCUIMiNa6AMffqubMvlvEwnHKJniYubyCqrlFKb8sh7poOE7dlbbOqqBbl9uIJL/oLVZFbomTnDwHKnnqkWCM5pt+wgHTO0gpKKjI6ewNVDYhd0gn0HtpqvFxYnsN5/bVEoskGDUxl9krKxk7u3DIJsDhQhKE6Ddaa2rOt3JsSzVXTzRisSomzS8xV4mjR9ZVok6YrrqN1e001wZorvXjb73dkGzPtJJX6iSvLJuC8hyKx7jIcIysjoQhf5Qzu03psr05hLvQwayVlUxdUjbifheDhSQI0efisQQXD9dzdPN1Gqt9OHLszHisnBmPlpPtka6O4t4S8QSXjzZybEs1ty63kZFlY/rSUcxaWUFO3vAsRQ1WkiBEnwn5o5zaWcOJbTX4W8PklTqZvaqSyYtKZcCUuC+3rrRxbHM1l47Uo5Ri/Lxi5qyupHiMO92hjQj3ShBSphM90tm+sNfUIVdMyWPF705h9FTpwigeTOk4D6Uf8+BtDHJ86w1O777JhYN1lE3wMGf1aMbOknaKdJEShLgrrTU151pM+8LJJtO+sLCUOasqKSjPSXd4YpiKBGOc3n2T4291aadYkWynyJJr2r4mVUyiV2KROOcP1HF8azVNNX6yXHamPyrtC2Jg3dlOYXdYmbqkjJnLK8gtdqY7vGFDEoToEW9TkFM7aji16yZhf4yC8hxmraxg0sKSQT2oTQx/dVe8HHurmkuH60lozdgZBcxcXkGlVHE+MEkQ4q50QnPjXAsntt3g6vFGAMbNLmLWygpGTcxN+xxEQnTlbw1zckcNp3bWEGyPklviZMaj5Ux5uDTt040MVZIgxLuE/FHO7q3l5I4a2uqDOHLsTFs6ihmPlg/bwVpi+IhHE1w8Us+JbTeou+LFlmFh4oISZjxaLr2fekl6MQnANDrXXmzl1M6bXDrSQDyWoGy8hwVPjWPCQ8VY7cNvtLMYnqx2C5MXlTJ5USkN19s5uaOG8wducWZ3LUWjXUxbOopJC0qkUfsBSQliBPC3hTm37xZn9tTSWhcgI8vG5IUlTFtWTmGF9EYSw0M4GOPcvluc3lVDU40fW6aVCfOKmbqkjLLxHqkuvQupYrqHSDBG000/yoJZatGisNiUWWXLbsGWYSXDYb3roiqDVSwa5+rxJs7tq+XaqWZ0QlM2wcPUJaOYML8YuwxqE8OU1pq6q15O77rJxUP1RMNxckucTF5sShxDrQpVa00skiASihGPJojHEsSiCXRCoxPmeUe2ndyS++vZJQniHm5ebOWXf3+k2+NsdgsZThuZTjuObBuObLPucZY7g2xPJjm5mWTnZuLKd5CZbUvL1YpOaGovtXLuQB2XDtcTDsTI9mQweXEpU5eMuu8PkBBDVSQU49KRBs7suUntxTZQUD4pl0kLSxn/UDGZaaqCikXj+JrD+FrD+FtC+FrDBL1RAp1rhUcJB2KE/VEi4dTL6XY1eVEpqz8y7b5iGVQJQilVCXwfKAUSwAta66/dcYwCvgasAwLAh7XW9/wWv98EEfJHqb/qRWvzBZtIaOKxBIlklo6G40TDcSKhOJFAlFDyjxb0maUtg77ou/54tkwr7gIH7sIsPEVZ5Bab9ZDzSrPJzs3o0+ShtabhejsXD9Vz4VAdvpYwtgwLVXOKmLK4jPIpeTIKVQigrSHI+QO3OLvvFt6GIFabhTEzC5g4v4QxMwqwZ/ZtqToSitFaF6DlVoC2hiBt9ea2vSlEwBt51/H2TCtZLjtZrgwcOXYcTjuZThsZWTbsmVbsmVZsGVasdlPDYbFaUAqURZGTm3nfg1cHW4IoA8q01keUUi7gMPAbWuvTXY5ZB3wWkyAWAV/TWi+61+umqw0iEU8Q8Ebxt4bxtYbwNYdpbwrhbQribTTr+HZdBznDYSUvuSRkwagcCirMmsiO7J530dNaU3+tnctv13PxcD3exhAWi6Jyej6TFpQwdlahzIwpxF1oram/2s75g7e4cKieoDeCLcPC2JmFVM0tYsyMgl79/8TjCVpqAzTdaKepxixl23zT944lYVFmrQ9PkRN3oQN3QXJ98DxHZ+1DXyeonhpUCeJdASj1a+AbWus3uzz2L8A2rfWPk/vngOVa69q7vc5gbaTWCY2/Ldx5JdFyK0BzrY+mGj8hX7TzuJz8TIoqXRSPcVM8xtx2XX0rHktw80IrV443cuVoA76WMMqiqJySx/h5xVTNKepVkhFCmNX8ai+0cvFwPZferifYHsVqs1A5NY9xc0yy6Dp7QDyWoPmmn/prXuqvemmo9tF000ciZr5HLTZFXmk2BaOyySvNJq/MSW6JE09R/60P/qAGbYJQSo0FdgAztNbeLo+/BHxFa70rub8F+HOt9aE7fv7jwMcBRo8ePe/atWsDFfoD61hZrKnGR2O1j8bqduqvt9NWH+w8xpFtw5UNdh2ioc1ONKKx2i2MnpZP1dwixs4slKQgRB9JJDS3LrVy+e1GLh9toL3ZLDFbkKfJyFT4whn4WiOda4NnOm0UjXZRVOmicHQOheWuu64kOJgNynEQSqkc4OfA57omh46nU/zIuzKZ1voF4AUwJYg+D7IfKaXI9mSS7clk9LQCAOI+P617DnBz12nqLzXjrYWgo5BAhouitssU+S4wapwTN3NwWueTYctP81kIMXwoBYX2NpyWY4yLHaK+uppb0SKaCqYTsdjICjRQZA2QX5bNqHnjKF6xgIzRo4d199m0JAillB2THP5Ta/2LFIfcACq77FcANwcitoEUa2kh+PZRAocPETh0iNCp0xCLYc/IYMLcuWQ9NhvHjPFkjB1L+FwugSOa4OEjNHz9n0FrlN1O1uzZOBctwrloIVlz5mDJyEj3aQkxJGitiV6/jn/ffgIHDhA4cIBYQwMA1oICih56iNEPzcE5dw46Hid44gSh4ycIHNqK7816fF8B26gynPPn45w3H+f8eWSMG4eyDK0SxL2ko5FaAd8DmrXWn7vLMU8Bn+F2I/XXtdYL7/W6g7UNokMiGCR09iyhU6cJnThB8OhRIskqMWW345g1C+f8+WQ/vJisuXOxZN591tR4ayuBI28TOHiQwP79hM6cMQkjKwvnvHlkP/ww2Y8sIXPy5GF9dSNEb8Wam/Hv2Yt/7x78e/cSu2maNa1FhWQvXIRz4UKcC+abL/q7/O9orYlcuYJ/714CB83FXbzRzGNmcbvJmjWLrNmzcUyfjmP6dOwlxQN2fvdjULVBKKWWAjuBE5hurgCfB0YDaK2fTyaRbwBPYrq5fuTO9oc7DZYEoaNRItU3iFy9QvjCRcLnzhE6f47IlasQNwvaWwsKyJozx1z9z52DY9aseyaE7sTb2ggcOoR/7z78e/cSuXTJvE9hIdlLHiZn6VKyH3kEW0FBX5yiEEOGjkYJvP02/p278O/eTei06SxpcbvJXrQI58OLyV68+J4Jodv30JrI1asEj7xN8OhRgseOEb5wgY7GCmthIY7Jk8mcNMlsVePIGDcOq3twzBk1qBJEfxmIBKG1JuH3E29qIlZfT7SunljdLaI1NUSqbxCtriZy4wbEYp0/Yy8vNx+KKZPJSl5R2EpL+/XKPlpXZ66Sdu/Gv2cP8eZmUArH9OlkL1uK67HHcMycibIOzl4VQjyI6K1b+HbswLdjB4G9+0j4/WCz4Zwzh+ylj5C9ZAmO6dP79fOf8PuTNQanCJ0+Q/j8ecIXL6Ijt8c/WAsLyaisxF5ZQUZFBbbSUuylpdhKSrAVFmLNzR2Q/1FJEPcQqa6m+Qc/gHgCHY9BPEEiHEKHwiRCQRI+P4l2L3FvO/GWlnf8gTtYPB4yKiqwV1aSMXo0GVXjyKyqIqOqCmtOeuc60okEodNn8O/aiW/HToJHj0IigTU3l+xHl+FavpzsZcuwulxpjVOI+6XjcUInTtC+dRu+7dsJnz0LmPaBnGWPkrNsKc7FD2PNyU5vnLEYkevVRK5eIXL5MuErV4jeqCFaXU301i1IJN75A0ph9XiwejxYXC6sbhcqy4nF4UA5MlE2O8pqAYuVrNmz8GzYcF9xSYK4h+DRo1z/g4+B1WqytdWCJdP8ASyZDvOHcbmwuFzY8vOw5uVjzc/HVlzUme3TnQR6I97aim/3bnzbt+PfsZN4a6u5upo3D9eqleSsXElGRUW6wxTinhKBAP49e2jf8ha+7dtNKdlqJWvuHFzLl5Pz2GNkTJgwZNrgdDRKrKGBaF0dsbo6Yk1NxJuaibU0k2jzEm9vJ+H1kgiFSISC6FAYHYtBPI5OJHCvXUvZ3/zP+3pvSRAiJR2PEzx2DN/WrbRv3Urkomm7yJw4kZzVq3CtWo1j+rQh808mhrdYYyPtW7fi2/IW/r170eEwFrebnGXLyFmxgpylj2DNzU13mEOOJAjRI5Fr12h/ayu+t94icPgwJBLYSktxrVqF6/HVOOfPR9lkCg8xcCLV1bS/uZn2zZsJvv02aI29vJyclStxrVqJc948lF0Giz4ISRCi12ItLfi2bad982b8u3ahw2GsubnmH/OJx8leskTGXIh+Eb54Ee8bb9D+5mbCZ84AkDllCq7Vq3GtXiXdt/uYJAjxQBKBAL5du2h/czO+rVtJ+HxYcnLIWbEC95onyF66FItjaM2xLwYPrTXhc+fwvv467W+82dlNO2vuXFyPP47r8dVkVFZ28yrifkmCEH0mEYkQ2LsX7xtv4Nu8hXhbG8rpxLV8Oa4n15Dz6KOSLES3tNaEz57F+9rreF97lei162Cx4FywANcTj+Na/figH2A2XEiCEP1CR6MEDh7E++prtG/eTLylxSSLFStwr32S7GXLHmgAoBhetNaEz1/A++ortL/6mplJwGIhe/EiXE+swfX4ahnMmQaSIES/07EYgQMHTLJ4803ira1YsrNxrV6Fa+1acpYsQUmbxYgUvnwZ7yuv4n3lFSKXL99OCmueNEkhXyadTCdJEGJA6WgU//4D5krxzc0kvF4sHg/uJx7HvW4dzoULZRT3MBe5caMzKYTPngWlcC5YgHvtk7ieeEJKCoOIJAiRNjoSwbdnD96XX8G3ZQuJQABrUSHuNU/ifmodWXPmSI+UYSLW0ID31dfwvvwywWPHAMiaPRv3U0/hWrNG2hQGKUkQYlBIhEL4tu/A+/LL+LZtQ0ci2MvLcT/1FO71T+GYNCndIYpeinu9tL/5Jm0vvURg/wFIJMicOhXPU+twPbmWjIrydIcouiEJQgw6cZ+P9jc3433pJfz79kE8TuakSbjXr8ez/inso0alO0RxF4lwGN/WbXhffgnftu3oaBT76NF41j+F+6mnyBw/Pt0hil6QBCEGtVhjo+nu+NJLZjJBIGv+PDzrN+B+co1MnzAI6HicwMGDtG3aRPvrb5Dw+bAWFeJZtw73+vU4ZsyQqsIhShKEGDIi1dV4X36Ztk0vmQFTdjs5y5bheXoDOcuXyxiLAaS1JnzmDG2bXsL78svE6utNz7QnnsCzYT3ORYuks8EwIAlCDDmdX04bN5kvp4YGLDk55svp6Q04FyyQL6d+Eq2poe2ll2nbtNFM4ChJeliTBCGGNB2PEzhwgLaNm2h/4w0Sfj+2khLc65/C8/TTOCZPTneIQ168rQ3v66/TtnEjwUOHAch66CE8T2/AtWYNtry8NEco+oskCDFsJIJBfFu30rZxE75duyAWI3PSJDxPb8C9fj320tJ0hzhkJCIRfNu34924yfQqi0bJGDcOzzNP416/XtYFGSEkQYhhKdbcjPfVV/Fu3GT63SuFc9EiPBs24HricVklLwWdSBA8csRU3b32GgmvF2thIe51a/FseBrHjOnS2DzCSIIQw17k2jXaNm6i7aVNRK9dR2VkkLNypak3X7p0xE/zEb540TQ2b9pE9OZNVFYWrtWr8Ty9geyHH5Z1PkYwSRBixNBaEzp+nLZfb8T76qvEW1qwejy41j6JZ8MGsubORVks6Q5zQETr6vC+/AptL20ifPqMmQNpyRI8G9bjWr0aS3Z612gWg4MkCDEi6WgU3+7deDe9RPuWLehQCPuoUbjXrx+2I7c7RzZveonA/v2gNY4ZM/BsWI973TpsRUXpDlEMMpIgxIiX8Ptp37KFtk0v4d+zx4zcnjzZTPOxbt2QnhIiEQrh27bdTGGyfbuZwmTMaDxPrce9YT2Z48alO0QxiEmCEKKLWFOTmVSu68jtuXNxr1uH+8k1Q+IqW0ci+PfuxfvKK7Rv3kLC7zeTIK5di2fDBhnZLHpMEoQQdxG5cQPvy6/gffllwufPd65q5l67FtcTjw+qtQp0LIZ//37aX3sN7xtvkmhrw+J243p8NZ7162UadXFfJEEI0QPhixdvL2xz9SpYrWQvWmgWtlm9Ki1rGOhoFP+BA7S//oZZiKmlBYvTSc6qVbjXrSXnkUdGfA8t8WAkQQjRC1prwufO4X3ttXcsjemcPx/X44/jWrWyX2ebTYRC+PfspX3zZnxbuqz7vWKFWfd72TKZ7kL0GUkQQtwns47yedpffwPvG6+buYkAx/Tp5KxcQc7y5TimTXvg+v5YUxO+HTvxbd2Kb9cudCCAxeXCtXIFrieeIPuRRyQpiH4hCUKIPhK+cgXfli20b95iRm9rja24mOxlS8lesoTsxYt7VBWlIxGCx47h37sP3+5dhI6f6HytnFUrca1eTfaCBVJ9JPqdJAgh+kHnVf+2bfj37iXh9QKQUVWFY+pUHFOnYK+oBIspXSS8XkJnzhI6e5bQ6dPoYBAsFhwzZ5Dz2GO4li8nc+pU6X0kBpQkCCH6mY7HCZ0+jX/3HoLHjxM6e4bYzdp3HWdxOsmcMgXHtGlkL16Ec+FCrG53GiIWwrhXgpAJWIToA8pqJWvmTLJmzux8LN7aSrSurnPf4nRiLy8fMVN9iKFPEoQQ/cSamyvLpYohTS5lhBBCpCQJQgghREoDniCUUt9VStUrpU7e5fnlSqk2pdTR5PbXAx2jEEKI9LRB/AfwDeD79zhmp9Z6/cCEI4QQIpUBL0ForXcAzQP9vkIIIXrnnglCKRXp4dbWx3E9rJQ6ppR6VSk1/R7xfVwpdUgpdaihoaGPQxBCiJGtuyqmGLC2m2MU8Ou+CQeAI8AYrbVPKbUO+BUwMdWBWusXgBfADJTrwxiEEGLE6y5B7NRab+/uRZRSu/soHrTW3i73X1FKfUspVai1buyr9xBCCNG9e1Yxaa3XdNxXSs1SSr1r5jCl1HSt9bq+CkgpVaqSk9EopRYmY2zqq9cXQgjRM73pxXQU2KOU2qC1buny+F6gx5PJKKV+DCwHCpVSN4AvAHYArfXzwHuBTymlYkAQeFYPlwmjhBBiCOlNgggAp4C9SqkntdZXk4/3aupJrfX7u3n+G5husAOiIdDA1uqt73rcoiwolLlVqnPfqqxYlAWLsmC1WLGq5GaxYlM2bBabuW+xYVM27Ba72ax2bJY79pVNZu4UYhiJJ+JEE9HbW/z2/VgidnvT5jau48QTceI6TiwRI6ETxHQMrTVxHUdrTUInSOgEGt15e+c18xj3GBaVLerz8+lNgkhorT+hlPorTJJYr7U+DAzpq/vq9mq+tO9LaXv/DEsGGVaz2S12Mq2Znfsd9zOtmWRaM3FYHWRYM3DYHJ37mTZz67A5bt/aHGTZsnBYk7fJ/SxbFnaLXZKSGBHiiTiheIhgLEgwFiQUC3XedjweioU690OxEOF4uHO7cz8cDxOJRzpvI/EIkUSk8340ESWu42k516fHP532BAGA1vpLSqlrwJtKqQ/1eUQDbGbhTLb+l9slCK11Z4buuE2QIJFImNtkNo/ruLlNZv+OK4COLa7j77hq6LgfTUQ7P0wdH7BoPHr7g5fc7/wwxsJ4w953flBjYULxENFEtNfna1XWzmThtDvNrc1Jlt3cOm1OnHZzm23Pfsf9ji3HnkNORk7nvkXJjC3iwWitCcfD+KI+fBEf/pgff8SPL+rDH/UTiAYIxALmfizQuR+MBt+5Hwt2Ph5JRHodh0VZOi+0Oi7MMqwZnRdn7gw3dqu98zm7xW4u6JIXenarvbOWoOOiz265XXtgs9g6axc67luUpbMWoqNGwmKxYMHSWVvRUYOhlOq87fp/l2Htn4WlepMgOi87tdbfV0rVAi8CWX0e1QCyW+0UZhWmO4z7Ek/EzZVO/PZVUDgWNldGKa6c3rVFg53/UG2hNm7Fbt3+J4wGiOlYj+LoSBquDFfnrTvTjcvuwpXhwpPpwZ3hxp3pJjczF0+GB0+mh9zMXOxWez//lsRASegE7ZF22sJtZou0dd73Rrx4I17aI+3v2nxRH76oj1iiZ5+3jouYbHu2ubixZeFxeCizld2++Ek+3rF1lqq7lKY7SuFZ9tv3bRap9u2qNwni97ruaK3fVEqtAN7TtyGJnrJarDgt5p+lP0TiEfxR/7u3mN9c5UX9nVd8Hf/o7ZF2GoONXGm7gjfixRf1kdCJu75Htj2b3Mxc8h355DnyyMvMIz8rnwJHAQVZBRRmFVLoKKQwqxBPpkf+eQdYKBaiIdhAU7CJxmAjjcFGmkPN79haQi20hltpDbfe82+dZcvCneHGlWEuHIqcRYzzjOvc77jQ6FpKddqdZNtul14dNoeUWAdQjxOE1voXKR47Dhzv04jEoNHRFpLnyLvv10joBP6oH2/Ee/vKMrm1hFs6b1tCLTQEGjjXfI6mUFPKq0m7xU5RVhHFzmKKncWUZJdQ4iyhLLvMbDllFDgKJIn0UCAaoNZf27nd8t+izl9HfaDebMF62iPtKX82NzOXPEce+Y58qjxV5DnyyM3MNZvD3Loz3LdLj8mqGTG0dJsglFKf7+4YrfWX+yYcMdxYlKXzCrE8p7xHP6O1pj3a3nnV2nFbH6ynIdBAfaCe8y3n2Vmzk2As+I6fzbRmUpZdRnlOORWuCipdlVS4KhjjGkOlu5JMa2Z/nOaglNAJ6gP1XPNeo7q9mur2am6036DGV8NN301awi3vON6iLBRmFVLqLGWsZywLyxZS7CymwFFAkbOIwqxCChwF5DpysVvky34k6EkJ4vE79h8Buo6c1oAkCNFnlFKdV53jPOPuepzWGm/Eyy3/LW75b1Hjq6HWX0uNr4YaXw3HG4+/4wpYoSjNLmWcZxxVnqrO2wm5E8h15A7AmfWPaDzKFe8VLrVe4krbFS63XeZK2xWue68Tioc6j7NZbCZx5lQwvWA6o3JGUZZdxqicUZRml1KYVYjNIotMittUb8egKaVatNb3X+fQT+bPn68PHTqU7jDEINMWbqO6vZpr3mtc917nqvcqV71XudJ25R2lj8KsQibmTmRy/mQm5U1icv5kqjxVg+4LsyHQwJnmM5xvOc/55vOcaznHNe+1zu6VFmWhPKeccZ5xjHGPYYxrDKPdoxntHk2psxSrxZrmMxCDjVLqsNZ6fqrn7ufTP6THPYiRxZNpekzNKJzxjscTOkGdv45LbZe41HqJCy0XuNB6gR+d+VFn98gMSwaT8iYxrWAa0wunM71gOuNzxw9Y0qgP1HOq8RSnmsx2pukMTaHbs86U55QzMW8iq0avYnzueCbkTmCsZ+yIqkYT/et+ShDNWuv8fornvkkJQvSFWCLGNe81zjaf5UzTGc42n+V002nao6aqymF1MK1gGjMLZzKraBazi2ZTkl3ywO8bjAU53XSaYw3HONFwguONx6kP1ANm7EpVbhXT8qcxtWAqU/KnMClvEq4M1wO/rxD3KkFIghCiGwmdoLq9mpONJznZeJLjjcc523S2s6RR4ixhVtEsZhXOYkbhDKYVTLtn1+N4Is619mvmtRqOc6LxBOebz3eOO6nIqWBm0UxmFpptcv5ksmxDeriRGMQeKEGk6MX0P4C/7frAYOjFJAlCDKRoPMq5lnMcazjGsfpjHG88To2vpvP5suwyqjxVlGaXdna79UV8XG67zNW2q53JJduezfSC6cwums3sotnMLJpJvmPQXX+JYexBE8S7Z7J7J621Xnm/wfUVSRAi3ZqCTZxsPMm5lnNcbrvM5dbLNARvr3TosDqoyq2iylPF+NzxzCycyVj3WGk4FmnVp1VMg5UkCCGE6L17JQgZsy6EECKlbhOEUmqFUuprXfaDSql4l215fwYohBAiPXrSofvTwI+67IeB1cn7i4HPANv6NiwhhBDp1pMEMQ/4cJd9rbXeDaCUOoJJIEIIIYaZnrRBFGqt/V32OxcJ0loHgeI+j0oIIUTa9SRBBJVSFR07WutNHfeTjwdT/pQQQoghrScJYifw3F2eey75vBBCiGGmJ20Q/xvYpZTKxywxWgNUAM8Cvwss7b/whBBCpEu3CUJrfVgp9TTwLeBjmNlcFXAZ+A2t9eH+DVEIIUQ69GjeYq31m8BEpdREoAho1Fqf79fIhBBCpFWvJrbXWl8ALvRTLEIIIQaRezZSK6VO9ORFlFJH+yQaIYQQg0Z3JYgJSqn3Y9oc7mVs34QjhBBisOguQdQBPVnr4VYfxCKEEGIQuWeC0FqPHaA4hBBCDDIy3bcQQoiUetWLqYMyayiOBwqAJuCSHi4rDwkhhADuI0EopVYC/whcAVqBPGCsUuqPtdZv9W14Qggh0uV+ShBfAlZorZs7HkhOw7EJkAQhhBDDRHfjICpSPQx473isne67wgohhBhCuitBnAI8AEqp5VrrbZg5mfYrpXYALUA+ZsK+f+rJGyqlvgusB+q11jNSPK+ArwHrgADwYa31kZ68thBCiL7TXS+mgFLKmby/EUBr/UNgJfAqZtqNV4BVycd74j+AJ+/x/FpgYnL7OPDtHr6uEEKIPtRdCeLnmNLCRsCqlJqgtb6otW4D3rifN9Ra71BKjb3HIc8A30/2itqnlMpVSpVprWvv5/2EEELcn+4Gyn1GKfUsprrHAZxWSgWAt4HDye2I1vpcH8ZUDlR32b+RfEwShBBCDKBuB8pprV/UWn8IuAm4gDXAT4Fc4M+A430cU6rG7pRjLJRSH1dKHVJKHWpoaOjjMIQQYmTrcTdXrXVl8u7+5AaAUiqjj2O6AVR22a/AJKdUMb0AvAAwf/58GagnhBB96IGn2tBaR/oikC42Ah9SxmKgTdofhBBi4N3XVBsPQin1Y2A5UKiUugF8AbADaK2fx/SKWgdcxHRz/chAxyiEECINCUJr/f5untfAcwMUjhBCiLuQ2VyFEEKk1OsEoZT6olLK0R/BCCGEGDzupwTxl4Cz26OEEEIMafeTIGRSPiGEGAGkDUIIIURKkiCEEEKkJAlCCCFESveTIGRKCyGEGAGkkVoIIURK95MgvgeE+joQIYQQg0uvp9rQWsvcSEIIMQJII7UQQoiUJEEIIYRI6Z4JQilVkeKxNclbe38FJYQQIv26a4M4BXgAlFLLtdbbuD0P0w6llEoe8zZmbeo9/RWoEEKIgdVdFVNAKdWREDYCaK1/mbx9GFgKfB3wA+/rryCFEEIMvO5KED8H9iulNgJWpdQErfXFjie11jHgWHITQggxjNwzQWitP6OUehazBKgDOK2UCmCqlA4ntyNa63P9HqkQQogB1W0vJq31i1rrDwE3ARewBvgpkAv8GXC8PwMUQgiRHj0eKKe1rkze3Z/cAFBKZfR1UEIIIdKv1yOp76S1jvRFIGlTcxh+/AGw2MBiNbcZTrBnQ0Y2ZOVCVj44CyCnGFxl4C4Dz2hw5oOSqamEEH0kGoK2avDWgLcW2msh0ATBFgg0Q8SX3AIQj0AiDokYTHsa1n61z8N54AQx5DlyYdKa27/oeARiIfNHCLVCyxXzhwm1vvtnM1yQNxYKqqBwMhRNhuKpUDgJrDJMRAiRgtbQfgvqTkHDWWg8B40XoOWqSQh3smWZi9GsfHC4wVkIuU6wZt6+sC2Z0S+hSoIoGA9Pf7374+Ix8NebP6C3FlqvQ+s1aL4Ct07CmU2gE+ZYawYUTYGy2VA+z2zFUyVpCDESeW+amoqaw3Dzbbh1wpQKOjgLzEXl+JXmgjN3DLhHmS2nBDJz0ha6JIiestpu/9HKUzwfC0PTJag/DbeOQ+1xOPsyvP0D87zdCRXzYfTDMHoxVC42VVlCiOFDa1MquLYbru8zW1u1ec5ig5LpMHkdlM4y94umQHZBemO+B6X18Fj/Z/78+frQoUPpDuOdtDbFxprDUH0AqveZqwedAIsdKhZA1WMwfhWUP2SKikKIoaWtBi5tgcvb4MpOU9MAkFMKYx6GykVQPh9KZ4LdkdZQU1FKHdZaz0/5nCSIARZuN1cVV3bA1Z1w8yigTVvI+BUwcQ1MfByyC9McqBAipXjU/A+ffw0ubjYlBjAJYdyjZhu71FQXDYFOLPdKEFLFNNAyXSYBTHzc7Aea4fJWuPiW+bCd+iWgTOliyjqYsgEKJ6Q1ZCFGvJAXLr5pqo0vbIZwm2lrHLME5n7Q1AIUTx0SCaE3pAQxmCQScOsYnH8Dzr0CtUfN44WTYdozZiuZPuw+hEIMSoFm8394+tem+igeMT2IJj8Jk9ZC1fK0NiD3FaliGqpaq80H9Mwm0+ilE1AwAab/Jsx8r+lWK4ToOyGvKSWc/Lkp2SdiZszTtKdhynqoXDjs2golQQwHvno4+5KpgrqyE9Cm7/PM3zabJ1XXKiFEt2JhOP86nPipKb3HwyYpTP8Ns416aFiX2iVBDDftt+DUr+Dkz+DGQUCZRrHZz5pqqExXuiMUYnDTGqr3w9EfwelfQagNsothxm/CjPeaLunDOCl0JQliOGu6BCd+BsdfhObLZrzF1Kdh7u/AmKVgkVVlhejUdgOO/hiO/ajL/8sGmPVfYNxyM95phJEEMRJobcZaHPsRnPyl6WWRO8b0sJjzAfC8a/VYIUaGWDg5aPWHcOktQMPYZTD7/aZtYYSXuCVBjDTRIJx5yYzivrIdUDBhNcz7PZj0pEz5IUaGhnNw+Htw7McQbAZ3hSlZz/mAGaMgABkHMfLYs2DWb5ut5aq5cnr7h/CTD5q5XeZ+EB76Pcgbk+5Ihehb0ZDplnr43+H6XjNjwZSn4KHfhaoVw64HUn9LSwlCKfUk8DXACnxHa/2VO55fDvwauJJ86Bda6y/e6zWlBNGNeMwMxDv8H3DhdVMlNWEVLPgDmPiE/OOIoa3pEhz6Lhz9TzM1dv54U2Ke/QHIKUp3dIPaoCpBKKWswDeBx4EbwEGl1Eat9ek7Dt2ptV4/0PENW1abGeAz+UnTUHfk+6b4/eNnwVMJ8z5sShXyzySGinjMTHdx8DtmzILFZkoL8z9q2hikg8YDS0cV00Lgotb6MoBS6kXgGeDOBCH6i6cCVnweHv1TMxDv4L/BW1+C7V+F6e+BBR8bUd38xBDjb4Qj34ND/25mSnWXw4q/hIc+BK7SdEc3rKQjQZQD1V32bwCLUhz3sFLqGGYt7P+mtT515wFKqY8DHwcYPXp0P4Q6zFntt6fwaDhvrsSO/giO/wRGzYWFnzD9wm2Z6Y5UCLOWwv4XzCjneNhMirfmy2b67BHYPXUgDHgbhFLqt4E1Wus/SO7/LrBQa/3ZLse4gYTW2qeUWgd8TWs98V6vK20QfSTsM2Mq9v8LNJ43c8/M/31Y8FG5OhMDLx6Ds5tg3/Nmunx7Nsx5vynlFk9Jd3TDwqBqg8CUGCq77FdgSgmdtNbeLvdfUUp9SylVqLVuHKAYR67MHNNwPf+jpl53/7/Ajr+DXf9oShOLP2VKF0L0p2CLaSM78K/gvWHG9Kz5Msz5HbNOvBgQ6UgQB4GJSqlxQA3wLPCBrgcopUqBOq21VkotBCxA07teSfQfpcwSiONXmh4iB14wXWWP/wTGPAKLPw2T10rvJ9G3mi7Bvm+b3kjRgGlsXvd/zPgd+awNuAFPEFrrmFLqM8DrmG6u39Van1JKfTL5/PPAe4FPKaViQBB4Vg+XEX1DUcF4WPtV07B95AemVPGT34H8KpMo5nwAMrLTHaUYqrSGa3tg7zdNpwmr3UxAufhTZhU2kTYyklr0Xke98J5vQM0hsxrego+aRm1XSbqjE0NFPAZnNsKef4abRyAr33yOFnxMPkcDSKbaEP3n+n7Y83Uz143VDrPeB0s+K2tViLuL+E115d5vQus1M6jt4efM3EgZznRHN+IMtkZqMZyMXgSj/9PUHe/9pqk7fvsHpuvhkj80i7YLAeBrMG1ZB//VNEJXLrrdTVUGtQ1KUoIQfcvfaHqeHPiX218CS/8YJq6RL4GRquWqqUZ6+4dmZtUpT5mLh9Gphj+JgSZVTGLgdVQj7PkGtF2Hoqmw9HMw47dkNtmRou6U6R598hegLGZBqyV/CEWT0h2Z6EIShEifeNR8Qez+J6g/DbmjzZfE3A+aWWfF8HN9P+z6BzNPUkaOmefr4efAPSrdkYkUJEGI9EskzCyyO/8Bbhwwyzsu+YwZpT3CF2wZFrQ2Ayt3/F+4tsv0SFr8KVj4McjKS3d04h4kQYjBQ2u4ugt2/l/zheLINV8kiz4hXyRDkdampLDj76DmMLjKTAlx3u/J2JghQnoxicFDKRi3zGw3DsPOv4dt/9u0VSz8mKmKyC5Md5SiO4kEnPm1KTHUnTBTYaz/JzNoUiZ3HDakBCHS79YJU6I49SvTLrHgo+YqNKc43ZGJOyXicOqXpsTQcBYKJsKy/woz3yudD4YoqWISQ0PDeVOiOPH/wJp5O1HIqNr0S8RNZ4PtX4WmC1A0xawnMv09MkfSECcJQgwtjRdNojj+U7BmSKJIpzsTQ/F0eOzPYOrTMq5lmJAEIYampkuw4+/NDLIdieKRz8myqAOhoypp+1fNuiDF02H5n8OUDZIYhhlJEGJoa7pk6ryP/wRsDlj4cXjkj8CZn+7Ihp+OxudtXzFtDEVTYflfSIlhGJMEIYaHxgvmivbEz8wArMWfMr2eZAGZB6c1nHsVtn7Z9EoqnGwSw7TfkMQwzEmCEMNL/RnTNfb0r8HhMe0Tiz5pVsMTvdMxwO2tvzXjGPKrYPl/N1OiSOPziCAJQgxPtcfNFe/5VyG7CJb+iRmZbXekO7Kh4fo+2PJFuLYbPJWm8Xn2B8Aqw6NGEkkQYnirPghvfRGu7AB3hWlMlS+6u6s9Dm99CS68YaY8efRPzchnGeA2IkmCECPD5W3mirjmMBRMgJX/w9ShK5XuyAaHpkuw9X/ByZ+bqrlHPmemOJEpMUY0mWpDjAxVy2HcY2Z1u7f+Fv7fh6FsDqz+Aoxfmebg0qj9lmncP/J901146Z+YXmDSuC+6IQlCDC9KwdT1MHmt6Ra79cvwg/eY5LH6f8KouemOcOCE2mD312DftyEeMdNuP/qn4CpNd2RiiJAEIYYni9VMHDfjt+Dgv5lxFC8sh+m/aaqeCsanO8L+EwvDwe+Ycw62mN/Bir8c3ucs+oUkCDG82TLh4U/D3N8xy17u/Sac2Wh6Oz36Z8NrVHYiASd+Cm/9L7OK3/iVsOoLMGpOuiMTQ5Q0UouRpb0Otn8FDn/PzBy75A/NYLuhPobi4hZ48wtmkFvZbHj8i6ZaTYhu3KuRWoZIipHFVQLr/xGeOwDjV8C2L8M/PwSHvgvxWLqj673aY/D9Z+CHvwlhL/zWv8HHtklyEH1CEoQYmQonwPt+CB99E/LGwUt/DN9+GM6+YkYXD3at1+EXH4d/edSMa3jyK/CZg2ZdBpkaQ/QRaYMQI1vlQvj910zX2M1fgBffD2MegSe+BOXz0h3duwVbYdc/wL7nTY+tpX9sxjNIl1XRDyRBCNHRNXbSGjjyPTOT6b+uNL1/Vv4V5I9Ld4QQiyR7Jv0fkyRmP2t6Y3kq0h2ZGMYkQQjRwWqHBX8As94Hu79uej2d3mimF3/0v6VnenGtzboMW/4GWq5C1QrTAF02a+BjESOOVFYKcadMF6z8S/jDt82V+v5vw9fnmEFn0dDAxXFtD3xnFfzsI2DPhg/+HD70K0kOYsBIghDibtxl8Mw34JO7oGIhvPnX8I35cOxFM+agv9SfhR+/H/59LXhr4Zlvwid3woTV/feeQqQgCUKI7pRMhw/+DD70a1PN9MtPmN5DFzb3bY8n703Y+FnTm+rqLtP+8dnDMPeDsjaDSAtpgxCip6qWmzEGp35hZo39z9+CMUvNZICVC+//dQPNpmfSgX81a0Ev/ISZMym7oK8iF+K+SIIQojcsFjPWYOrTcPg/zHxH//Y4TFxjFtypSDkgNbVAs5lIb//zEG437R3L/wLyxvZX9EL0iky1IcSDiPjNF/yeb0Cw2ZQyHvkjM+343aqFWq7BoX8zkwhGfDB1Ayz/PJRMG9DQhQBZMEiI/hf2mek69vwz+Oshp9SMo6hafjtRNF+GEz+D6n2gLOb5Zf8ViqemNXQxsg26BKGUehL4GmAFvqO1/sodz6vk8+uAAPBhrfWRe72mJAgxKESDcP41kwguvGHWYeiqaCrM+m2Y8V7IG5OeGIXoYlCtKKeUsgLfBB4HbgAHlVIbtdanuxy2FpiY3BYB307eCjG42bNg+nvMFmyBhvO3n3PmQ+HE9MUmRC+lo5F6IXBRa30ZQCn1IvAM0DVBPAN8X5vizT6lVK5SqkxrXdsfAf3NplOcvuntj5cWoovG5CZE35o2ys0XNkzv89dNxziIcqC6y/6N5GO9PQal1MeVUoeUUocaGhr6PFAhhBjJ0lGCUCkeu7MhpCfHoLV+AXgBTBvE/QbUH5lXCCGGunSUIG4AlV32K4Cb93GMEEKIfpSOBHEQmKiUGqeUygCeBTbeccxG4EPKWAy09Vf7gxBCiNQGvIpJax1TSn0GeB3TzfW7WutTSqlPJp9/HngF08X1Iqab60cGOk4hhBjp0jLVhtb6FUwS6PrY813ua+C5gY5LCCHEbTKbqxBCiJQkQQghhEhJEoQQQoiUJEEIIYRIadjM5qqUagCuPcBLFDKy5kEYaecLcs4jhZxz74zRWhelemLYJIgHpZQ6dLcZDYejkXa+IOc8Usg59x2pYhJCCJGSJAghhBApSYK47YV0BzDARtr5gpzzSCHn3EekDUIIIURKUoIQQgiRkiQIIYQQKY2oBKGUelIpdU4pdVEp9RcpnldKqa8nnz+ulHooHXH2pR6c8+8kz/W4UmqPUmp2OuLsS92dc5fjFiil4kqp9w5kfP2hJ+eslFqulDqqlDqllNo+0DH2tR58tj1KqU1KqWPJcx7Ss0Irpb6rlKpXSp28y/N9//2ltR4RG2Zq8UtAFZABHAOm3XHMOuBVzIp2i4H96Y57AM55CZCXvL92JJxzl+Pewswq/N50xz0Af+dczLrvo5P7xemOewDO+fPAV5P3i4BmICPdsT/AOT8KPAScvMvzff79NZJKEAuBi1rry1rrCPAi8MwdxzwDfF8b+4BcpVTZQAfah7o9Z631Hq11S3J3H2b1vqGsJ39ngM8CPwfqBzK4ftKTc/4A8Aut9XUArfVQP++enLMGXEopBeRgEkRsYMPsO1rrHZhzuJs+//4aSQmiHKjusn8j+VhvjxlKens+H8VcgQxl3Z6zUqoceA/wPMNDT/7Ok4A8pdQ2pdRhpdSHBiy6/tGTc/4GMBWzXPEJ4I+01omBCS8t+vz7Ky0LBqWJSvHYnX18e3LMUNLj81FKrcAkiKX9GlH/68k5/xPw51rruLm4HPJ6cs42YB6wCsgC9iql9mmtz/d3cP2kJ+e8BjgKrATGA28qpXZqrb39HFu69Pn310hKEDeAyi77FZgri94eM5T06HyUUrOA7wBrtdZNAxRbf+nJOc8HXkwmh0JgnVIqprX+1YBE2Pd6+tlu1Fr7Ab9SagcwGxiqCaIn5/wR4CvaVNBfVEpdAaYABwYmxAHX599fI6mK6SAwUSk1TimVATwLbLzjmI3Ah5K9ARYDbVrr2oEOtA91e85KqdHAL4DfHcJXk111e85a63Fa67Fa67HAz4BPD+HkAD37bP8aWKaUsimlnMAi4MwAx9mXenLO1zElJpRSJcBk4PKARjmw+vz7a8SUILTWMaXUZ4DXMT0gvqu1PqWU+mTy+ecxPVrWAReBAOYKZMjq4Tn/NVAAfCt5RR3TQ3gmzB6e87DSk3PWWp9RSr0GHAcSwHe01im7Sw4FPfw7fwn4D6XUCUz1y59rrYfsNOBKqR8Dy4FCpdQN4AuAHfrv+0um2hBCCJHSSKpiEkII0QuSIIQQQqQkCUIIIURKkiCEEEKkJAlCCCFESpIghBBCpCQJQoi7UEr9VX9Oi52clnnRQL6nEL0xYgbKCdETSikL8Ang08A0wKKUug78FPhbrXVrH73PTMycSAcG6j2F6C0ZKCdEF0qpbwHvA34fM/f+Y8AfY+aqsgCLktNLP+j7/BVQprX+9EC9pxC9JVVMQiQl587/JPAnWutfA1EArfXbwG8CczDThPeF9wC/GuD3FKJXJEEIcds4zJw9h+98Qmt9DWjCrGDWK0qpOcnJ0zr2Ryffa2tv3lMp9UGl1N7ktqq3cQjRW5IghLjtCmb+/Het5Zv8Ui/g/mYDfS9Q3GX/N4BXtNbRnr6nUioX+BPMZG3rgX9USlnvIxYhekwShBBJyamRXwD+QSm1DjNLKEqpGZhpwY8Av0w+tk4p9ZZSardS6u+VUhal1D6l1FilVIVSakfydhXwHPBnSqlvJt/qPcCvevmei4CdWutwcs2Oq5hFcIToN9KLSYh3eg4zJfbfYdYPsGCmUX4R06MoopSaDPx34AmtdVAptQmztOPngX8A8oCPaq1vADeSC9Us01prpVQBsIB3Lu3ak/csAFq6/EwLpnQhRL+RBCFEF1rrOPAtzPoYfwWs1lo/dsdh78OsRPdqcg2NAkw7wjFMFdBzWusLAEqpfKBJ3+4uuAHYprX29fI9mzCJp0Ne8jEh+o10cxWil5RSfwsc11r/NLlvBzIxVUHfBz6ste5YyWwl8LTW+nPJ/V8BL2mtv9PL98wF3gIeBrKB7cCcZHIRol9IG4QQvfcC8JxSaptSajOmO+pPgC9prX+AqVbakDz2MrBaKbUlOSBuD8n2h95IDpb7J2Ab8DKmW6wkB9GvpAQhhBAiJSlBCCGESEkShBBCiJQkQQghhEhJEoQQQoiUJEEIIYRISRKEEEKIlCRBCCGESEkShBBCiJT+P3L0blk1XIgVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "for i in range(n_eig):\n",
    "    plt.plot(phi, spec[i,:]- spec[0,:])\n",
    "    \n",
    "plt.xlabel(r\"$\\Phi_{ext}/\\Phi_0$\", fontsize=13)\n",
    "plt.ylabel(r\"$f_i-f_0$[GHz]\", fontsize=13)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell shows the spectrum from the figure 4 of the paper, which is the same spectrum that SQcircuit calculated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"pics/zeroPiQubit_pic3.png\" width=\"250\" align = \"left\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eigenfunctions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate the eigenfunction in the phase space by `eig_phase_coord()` method for the ground and first excited state. To calculate the eigenfunction at $\\Phi_{ext} = 0.9\\Phi_0$ similar to paper, we set back the flux of our loop to $0.9\\Phi_0$ and diagonalize the `zrPi` again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "loop1.set_flux(0.9)\n",
    "zrpi.set_trunc_nums([35, 11])\n",
    "_, _ = zrpi.diag(n_eig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a range for each mode\n",
    "phi = np.pi*np.linspace(-1,1,100)\n",
    "theta = np.pi*np.linspace(-0.5,1.5,100)\n",
    "\n",
    "# the ground state\n",
    "state0 = zrpi.eig_phase_coord(0, grid=[phi, theta])\n",
    "    \n",
    "# the first excited state\n",
    "state1 = zrpi.eig_phase_coord(1, grid=[phi, theta])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\theta$')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUIAAACeCAYAAABHJgN8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOfUlEQVR4nO3df6gl513H8fcnN/sj2R/ZJLvSJtl1S0nFoG2DS1ArpWijaQ2JFYSWWgTBpWAwpRbbtGKpkj+kUKKSP1xIwGK0CGlpixHbqqH2j9TcDZs0zSayFEO2jcn+cLO72Ww2N/v1j3Pm7rO7M3dm7p4zc+Y8nxdcmDtnzpkvy90P32eeeeYoIjAzy9llfRdgZtY3B6GZZc9BaGbZcxCaWfYchGaWPQehmWXv8r4LKLN169bYuXNn32WY2ZzZu3fv4YjYduH+mQzCnTt3sri42HcZZjZnJD1ftt9DYzPL3kx2hJavSa90kjTRz7P55CC0mTCtpZ7p5zoUrYqD0HpTFn7T7AgdilbFQWidqgq6Yn/d6xcqC7Qm4VfsdyAaOAitA03CrSwIm7yvLgirgq7Y7y7RwEFoU9Im/NLts2fP1h5bJg2xyy5rfjOEh84GHQahpAVgEfhxRNze1Xmtf3XhB+cCMA3CJqFYBFZV+FXtd9BZqsuO8G5gP7C5w3Nax+o6t7LwS7fffPPN5X3pdlUoFkGXBt7CwkLpudNjyup0d5ivToJQ0g3AbwL3Ap/s4pzWnbrZ37ouEM6F3tLS0vK+dDsNxbIgTMOvycRKm6G2J1bmX1crS+4D/gQ4W3OcmVnnpt4RSrodeDki9kp63wrH7QZ2A+zYsWPaZVmHmnSERfd35syZ5X1vvPHGRa9f+BlFR7hmzZrSc1dNohT723R5Hi7Pry6Gxu8B7pD0QWA9sFnS30fE76YHRcQeYA/Arl27/I1SM67N/YB1w2EoD8LTp08vb6ehmH5eMSROPyuVht9qg7DsGIfifJl6EEbEPcA9AOOO8FMXhqDNtyZBWARgGn6nTp1a3n799ddLP+Pyy0d/wldcccXyvjSY0muHZdcRJxFoDsXh832E1libVR9thsNwLgjT8Dtx4sTy9muvvVb6GcWQON2XBl46ZE63i2OqZpJXO1niUBymToMwIh4FHu3ynHbp2qz/Leu0qrrAdLhbdHxpEB4/fnx5++TJk6WfsW7duovOV3SJAGvXri3dLo5ps3qlbch5tnk43BFaqUt9+EFVR1h2XRDOdXxp4B07dmx5+5VXXil9XzEkTsOmCEeAK6+8svR9RTjX3ahteXAQ2qrU3SfY5rognLs2mHaBR48eLd1OO8kNGzYA5w+H0/DbtGlT6bnLVrJM64ZrB+zscxDaskk8Aqtu1rhqaFx0hOl1wSNHjixvHzp0aHk7DdDNm0cLldIu8KqrrlreTidZ0vPVdYRtOOiGz0Foja3mWiE0u0ZYBGFVR/jSSy8tb6dBWARd2gVu3br1otehfKXKpK8ROhSHyUGYuWk9GTpV9SCFuo4wDcK0Iyzr8tIu8NVXXy09tmxo3OaxX22Drex9DsrZ5CC0WqtZS9x2aFwEVtVkyeHDhy86Nv3stJNMZ57T7rHNZEldYDnQ5ouD0Kai7X2EdbfPpKGY3nRd3AaTBmjVipSyJ9hM8rooOBSHyl/naWbZc0doU9emO0yHvWl3mHZ86SqTYpIk3ZcOh+ueaTiJWWMbPgehlZrErSR1r9c9j7Bsdrhqu+pJNU2eeF2n7p5CGz4Hoc2Eug6tzeRM2eearcRBaL0pe1ZgurojXTNct52uLKl79NaF22YOQivV5rH2Ve8r21f1oNQi0NIHI6xfv355O71hOv2MYq1xurIkfcpMVRDW1VZ2jMNzfjkIbSqahF/Z47LS5woW64jh3FI6OD8si7XE6bFVoVjVNV4qB+TwOQhtVVb7RUhVQ98ivNJAS1eLXHPNNcvb6X2CV199NXD+wxXSMK0KwrIuz4GWLweh1Spba1v2erpdFTBVD00thsEbN25c3peGX9VqkS1btgDnh2Y6jE67wzR4i0CuGw5XWe0j/m02OQgzt9prgW0+t244DOe6uHQInAZh1Zc6FccXgQjnd5XpMLpsaNzmGuEkOBxnk4PQGmsSmm06wjSkyoLw2muvLT1fGoRFB1kMkdN9cP6ES9oRlj2qvy4U3QXOLwehLZvETHHZE1fadITptb60I0w/L71huhgGV3WEdbPJkwguh9/wea2xmWXPHaGVqusO656xN4mOML0umA5r046wGPqmkyWTuEY4iWHyao61fjgIbeqqQjENwrLbZ9LAS48t+xa7qttnyr65Lq2jze0zbQPNATgcDkKr1eb2mbJ9VZMlZfcRpre+pOdLA63se43T91UF4Wo7Qpt/DkJrbBI3UVfdUF0EVhpiqbQjTM9dhFs6O9xmZckkZoIdmsPnILRVafM9HG2uF6bdXvq+uiBMX68aDrdZWeLwy8vUg1DSduDLwFuAs8CeiPiraZ/XZkeb7rAqVNNAS48pPiN9vcn64tU8SMHhN7+66AiXgD+OiCckbQL2Svp2RDzTwbmtA3X3EVaFW1V3WHZs1UNVi89I31+13WYpYFUdNp+mHoQR8SLw4nj7hKT9wPWAgzATVUFS9wSYqiAsO6aq02zybEKzTq8RStoJ3Ax8v8vzWnfafBl6XRCmr9cFYZPhd11HWPa5lofOglDSRuBh4BMRcbzk9d3AboAdO3Z0VZZNSZt1yel2GlxNwrTNpMckb5K2+dJJEEpawygEH4qIr5YdExF7gD0Au3bt8hdNzKmqgCwCsMm3yrWZpfa9gdbE1Ncaa/TX9wCwPyK+NO3zmZm11UVH+B7gY8APJO0b7/tsRDzSwbltBjSZQS6bbW7yeSvta/J57hINupk1/h7gvzYD2t2IPYlzXMoxlg+vLLHetHnCzaTPZ5aqDcLxLS9/CLwdOArsA74ZEc9PtTLLSl1ItVnbbNZWk8mSrwPPAvcDtwLvAr4r6X5J61Z8p9mESFrxx+xSNAnChYh4ICL+DTgaEX/AqDv8H8a3u5iZDVmTIPyOpLvG2wEQEUsR8UXgl6ZWmZlZR5pMlnwSuEfSInDdeAXIKUYheGSaxZmZdaG2I4yIsxFxL/BeRkvg3gL8AvA08IHplmdmNn2Nb5+JiFPAN8Y/ZmZzw1/naWbZcxCaWfYchGaWPQehmWXPQWhm2XMQmln2HIRmlj0HoZllz0FoZtlzEJpZ9hyEZpY9B6GZZc9BaGbZcxCaWfYchGaWPQehmWXPQWhm2eskCCXdJuk5SQckfaaLc5qZNTX1IJS0wOg7kT8A3AR8RNJN0z6vmVlTXXSEtwAHIuJHEXEG+ApwZwfnNTNrpIsgvB54Ifn94HifmdlM6CIIVbIvLjpI2i1pUdLioUOHOijLzGykiyA8CGxPfr8B+MmFB0XEnojYFRG7tm3b1kFZZmYjXQTh48CNkt4maS3wYfzdyGY2Qxp/wftqRcSSpLuAfwUWgAcj4ofTPq+ZWVNTD0KAiHgEeKSLc5mZteWVJWaWPQehmWXPQWhm2XMQmln2HIRmlj0HoZllz0FoZtlzEJpZ9hyEZpY9B6GZZc9BaGbZcxCaWfYchGaWPQehmWXPQWhm2XMQmln2HIRmlj0HoZllz0FoZtlzEJpZ9hyEZpY9B6GZZc9BaGbZU0T0XcNFJB0Cnq94eStwuMNyLpXrnS7XO13zVu9PR8S2C3fOZBCuRNJiROzqu46mXO90ud7pyqVeD43NLHsOQjPL3hCDcE/fBbTkeqfL9U5XFvUO7hqhmdmkDbEjNDObqEEGoaS/kPSUpH2SviXpur5rWomkL0p6dlzz1yRt6bumlUj6HUk/lHRW0szOGEq6TdJzkg5I+kzf9axE0oOSXpb0dN+1NCFpu6T/kLR//Ldwd981rUTSekn/JenJcb1faPX+IQ6NJW2OiOPj7T8CboqIj/dcViVJvw78e0QsSfpLgIj4dM9lVZL0s8BZ4G+BT0XEYs8lXUTSAvDfwK3AQeBx4CMR8UyvhVWQ9F7gJPDliPi5vuupI+mtwFsj4glJm4C9wG/N8L+vgA0RcVLSGuB7wN0R8ViT9w+yIyxCcGwDMNNpHhHfioil8a+PATf0WU+diNgfEc/1XUeNW4ADEfGjiDgDfAW4s+eaKkXEd4GjfdfRVES8GBFPjLdPAPuB6/utqlqMnBz/umb80zgXBhmEAJLulfQC8FHgz/qup4XfB/6l7yLmwPXAC8nvB5nh/6hDJmkncDPw/Z5LWZGkBUn7gJeBb0dE43pnNgglfUfS0yU/dwJExOciYjvwEHBXv9XW1zs+5nPAEqOae9Wk3hmnkn0zPTIYIkkbgYeBT1wwEps5EfFmRLyb0YjrFkmNL0FcPrWqLlFEvL/hof8A/DPw+SmWU6uuXkm/B9wO/FrMwIXZFv++s+ogsD35/QbgJz3VMpfG19oeBh6KiK/2XU9TEXFM0qPAbUCjyamZ7QhXIunG5Nc7gGf7qqUJSbcBnwbuiIhTfdczJx4HbpT0NklrgQ8D3+i5prkxnnx4ANgfEV/qu546krYVd2NIugJ4Py1yYaizxg8DP8NoZvN54OMR8eN+q6om6QCwDjgy3vXYjM9yfwj4G2AbcAzYFxG/0WtRJSR9ELgPWAAejIh7+62omqR/BN7H6OkoLwGfj4gHei1qBZJ+BfhP4AeM/p8BfDYiHumvqmqS3gn8HaO/hcuAf4qIP2/8/iEGoZnZJA1yaGxmNkkOQjPLnoPQzLLnIDSz7DkIzSx7DkIzy56D0Myy5yC0QZP025L+uu86bNgchDZ0NwNP9F2EDdvMPnTBbCWS3gHcD/wicETSloi4r9+qbKi8xM4GR9I6Rs/G+xjwdeCXgWeA6yLidJ+12TB5aGxDdCvwJKPHbh2PiP8FTjNacG/WmoPQhujdjJ6K8i7gKUk/BZyIiFd7rcoGy9cIbYiOA+9k9LTvp4AvMLpeaLYqvkZogyPpGuBrwM8D/8foi5v+dBae/G3D5CC0wZL0JPCrEXGk9mCzFfgaoQ3SeOZ4k0PQJsEdoZllzx2hmWXPQWhm2XMQmln2HIRmlj0HoZllz0FoZtlzEJpZ9hyEZpa9/wcAXFPe3I8CJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5, 2))\n",
    "plt.pcolor(phi, theta, np.abs(state0)**2,cmap=\"binary\",shading='auto')\n",
    "plt.xlabel(r'$\\phi$')\n",
    "plt.ylabel(r'$\\theta$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\theta$')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUIAAACeCAYAAABHJgN8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQXElEQVR4nO3dfaxl1VkG8Oe5l2GAmYHh45IUZvA2DRqNtpBM8KPGGAQltQE1McGPxsRE0kQijRprxdhUQ6JpUlHDH04KiUaUmNCmVDH9sCW1f1C5Mw4IDDWkKWHaKndmGGaGYYBhXv84Z515585aZ+195uyvs55fcpN9993n7PeSuQ/v2mvvdWhmEBEp2VLXBYiIdE1BKCLFUxCKSPEUhCJSPAWhiBRPQSgixbug6wJirrrqKltdXe26DBFZMHv27DloZisb9/cyCFdXV7G2ttZ1GdKxNu5xJdn4OaQ/SL4U29/LIBQBFFLSHl0jFJHiKQhFpHgKQhEpnoJQRIqnIBSR4ikIRaR4rQUhyWWS/0XyX9o6p4hIFW12hPcA2N/i+UREKmklCEnuAPDzAD7dxvlEROpoqyO8H8AfADjd0vlERCprPAhJfhDAK2a2J3PcXSTXSK6tr683XZaIyEQbHeH7AdxO8tsAHgFwM8l/2HiQme02s11mtmtl5ZzFIUREGtN4EJrZx8xsh5mtArgTwFfM7NebPq+ISFW6j1BEitfqMlxm9gSAJ9o8p4hIjjpCESmeglBEiqcgFJHiKQhFpHgKQhEpnoJQRIqnIBSR4ikIRaR4CkIRKZ6CUESKpyAUkeIpCEWkeApCESmeglBEiqcgFJHiKQhFpHgKQhEpnoJQRIqnIBSR4ikIRaR4CkIRKZ6CUESKpyAUkeIpCEWkeApCESle40FIcifJr5LcT/I5kvc0fU4RkTouaOEcpwD8npntJbkNwB6SXzKz51s4t4hIVuMdoZl9z8z2jrePAdgP4NqmzysiUlWr1whJrgK4EcA32jyviMg0rQUhya0AHgXwETM7Gvn5XSTXSK6tr6+3VZaISDtBSHITRiH4sJl9JnaMme02s11mtmtlZaWNskREALQza0wADwLYb2afavp8IiJ1tdERvh/AhwDcTHLf+OsDLZxXRKSSxm+fMbOvA2DT5xERmZWeLBGR4mU7wvEtL78N4D0ADgPYB+DzZvZSo5WJiLSkSkf4OQAvAHgAwK0A3gfgayQfILm5yeJEYsxspi+RlCpBuGxmD5rZvwM4bGa/hVF3+G0Au5ssTsrRRqApKCWlShB+meTd420DADM7ZWafBPDjjVUmC6nPIdTn2qRZVWaNfxfAx0iuAbiG5F0ATmAUgoeaLE6Gq60Aaeo8o9tf0+fwP5fhywahmZ0GcB/JvwRwC4AbAFwO4FkA9zZanQxCnTBqIyBnPUcq/GKhl/u5DEvl+wjN7ASAx8ZfUrhc2DQZeE0Fb51wU/gtljbWI5QFkQqV2P4uh8Z16vRSHeHS0tLUn6feQ4ZDQSgzyYVNX4KwynbgQywWaH5fCMdp51YoDoeeLBGR4qkjlKxcd3f69Ompx1bpDmftnmIdqN9O1ZbrCGMdn9+nofFiURBKZXWGnD6AqggBUjdIYuerux34oPPby8vLyXqn7Q+1KRz7T0EoUbFwS3Vase1c9wXEr8nlrtOl3tvX8M4772S3Y3Wmwi8WaKnuUdcIh0lBKJVVGX6GsKnSEcbCpEoQeuE8qcB7++23J9unTp2Kbm+sAQAuuODcP41Ubbn7DxWO/acglIk6M8G5DszvSwVBrJPy+6qETSwIffi99dZb0e1wTKojvPDCC7FRqmPMzSAr/PpPQShRsw4/Q8BUCUIfJmE7FUypIIwFrw/CN998c7J98uTJc/b738nX498vnDtW78Y6UzVLvykIJSvWEfqg8MPMEEJ+XyrcYjO6sSFpqh5fhz+f7/x8+L3++uuT7TfeeOOc38Of2+8PNfuf++3Y9US/XefRPemGgrBwVZ7CyHWEsaGo35cKwk2bNk2tzQdMriNMDYdPnDgx2T527NhkO4Sif52vxwdrqMMPlzdvPrMUZ53bhxR+/aQglIncLTGp4bAPnjDk9Pv8e6W6p9yssd+OBXJqOOyD8OjRMx+n/dprr51Tpw86//uF/RdffPFk30UXXRR9Xa4jlH5SEEpWrCOMDYeBM0PRVBCmhr6hU6wyWRKrI1YDcPZwOIQfABw6NFpBzodmKtBCAG7dunWyb8uWLdF6/HZqEkX6R0EoUXVmZmMTEn6ffy8//IxNnKSGw3WGxj4Ijx8/Ptl+9dVXJ9sHDx4EcOZaIXD2cNefLwTgZZddNtnnz5eaGIrt0/XCftL/skSkeOoIZaLOY3O5WVrfafn38l1XbOIkNlu7cbvO0Nh3hEeOHJlsh47Q/9xf9/ND+Msvv/yc963SEWqZ/+FQEEpWbGhcJwi91PXCECx+6Jy6xhabtEkN1f01Qh+E6+vrAM6eSU5Nhlx99dUA0kGYukaoGeThaCUISd4G4K8ALAP4tJn9eRvnlbgqncqsHWEIIR8aqfP57jC8X+x54I1is8a+nlQQ+smSw4cPAzg7HC+55JLJ9rZt2ybboWv0v5M/36zXCKU/Gg9Ckss485nIBwA8RfIxM3u+6XNL3qwryuS6w9RkiZ8MiT0HnOqucrfPpILZd6a++wsBGAJxY81XXnnlZDvcguPfN3XDuIbGw9RGR3gTgBfN7FsAQPIRAHcAUBAOxCxPlvjQ8FKLIIRwSw0tc6vdVJnRjt1c7e8t9O/rrx2GMPW/U6pzVfgNUxtBeC2Al933BwD8aAvnlTlLdTuxUIyt7rLx2FiYVOkIZ33qxYdiCDffMfprkv7Y2PPTCr/F0kYQxi6KnPMvZ/x5yXcBwHXXXdd0TTKDOsPoeawYXaemOkuEAfHATl33iy0tpiHwYmnjPsIDAHa673cA+O7Gg8xst5ntMrNdKysrLZQlbSGZ/Wr73EtLS1haWoru27i/zpcMUxsd4VMArif5bgDfAXAngF9t4bwyZ1VWaI49IZJadCG2XeVZ49h2lWX2/a054XE6P3PtH7Hzx4bbfLTE1uJqPAjN7BTJuwF8AaPbZx4ys+eaPq/MT2xBhNTafCE0fJCk7h2MrelX5VnjWAinVpf24ebvEwzPCvtrgf75YX8rTQhL/zul1iNUQA5TK/cRmtnjAB5v41wyu1m7rlinFVvheeP+WChWWeQ0dkwq/PyN0T7oLr30UgBnX/dL3UcYAtS/b5VnohWKw6EnSwpU5UORYsenui4fhKF7St1nlxt+Vln52QvHp4LQh1sIPwC44oorzjmfP9YvsBAWXfDD6FRnm/tgeOknLbogIsVTRyhRsSFe7FogEF+5OXWrSWpyItYRVrkOF+sI/Tn8cHj79u2T7fDkiD/WX0MMCy0AZzpCP8xOPROtofEwKQglKneNsMo1uSA1NPYhFIKl7tA41JGqwS+m6oMwPDfsh8P+dWHoDJy5XphanSY3NPYUjv2kIJSJ1B9pCKHUrSg+0MJTGKknQWITK35/blZ543uHY/z7poLQd3nhGqb/uf89YkHoO8ZURxirM/XfVaHYHwpCiTrfjjD2mcXTXpcbGqeGnLHbdVJB6CdAQk2ppfp9aMY6wtStNBoaD5OCUKJyQZjrDlNBmHpdriNMrT4T6wh9Z5eaCQ71pT7Fzs8wh+uMsaG8f6+NdcpwKAgLV+VWmtjQOHX7THhdlSdLYp8TnOoIczX7cEwFoT8mnM/f5uPr8ZMssSD0x+YmSxSO/acglKjcH3SqswtSHWGqq4w9mpcKmNj9jlWeJokN0VMf8O6HwbkbquvMGisU+0lBKFm5oXFsJelUEOaeEKn7KXaxfamP5fTvEcI7NswGzu7+QijqEbvFpSCUidwwORVivpOqMxyMBUiVUIntT3WoqSAMNae61djsdt0gVCgOh4JQsmIh5cU6tCqf35sbfs87SFLhHTtf7BqohsOLS0EoWeEPucoQN/bzOueo+7pUOOeODUPieV+/lGHSs8YiUjx1hBIVG9rmnjzxx9Rdyr7OkDL38ZipoWruc0Zyi83WXStRneJwKAilsln/sOu8btZzVJmlrvMccG7xV4XcYlEQSlbuGuG012x8XdVznY8qEzU5uWeG1QUuFgWhVDaPgJl3HfN8Xe52HQXe4lIQykz6Eoo5uc60zr2KWkVmcSkI5bzNGjY58x5Sz7OTVPgtFgWhzFXseuL5vlddTV2TVPgtLgWhNKJKaMwjLGc9dxvvIcOhIJTO5MKmqaCsU4OUgW38Y6uL5DqAlxI/vgrAwRbLOV+qt1mqt1mLVu/3mdnKxp29DMJpSK6Z2a6u66hK9TZL9TarlHr1rLGIFE9BKCLFG2IQ7u66gJpUb7NUb7OKqHdw1whFROZtiB2hiMhcDTIISf4ZyWdI7iP5RZLXdF3TNCQ/SfKFcc2fJbm965qmIfnLJJ8jeZpkb2cMSd5G8pskXyT5h13XMw3Jh0i+QvLZrmupguROkl8luX/8b+GermuahuRFJP+T5NPjej9R6/VDHBqTvNTMjo63fwfAD5nZhzsuK4nkzwL4ipmdIvkXAGBmH+24rCSSPwjgNIC/BfD7ZrbWcUnnILkM4H8A3ArgAICnAPyKmT3faWEJJH8KwHEAf29mP9x1PTkk3wXgXWa2l+Q2AHsA/EKP//sSwBYzO05yE4CvA7jHzJ6s8vpBdoQhBMe2AOh1mpvZF80sfJL4kwB2dFlPjpntN7Nvdl1Hxk0AXjSzb5nZWwAeAXBHxzUlmdnXABzuuo6qzOx7ZrZ3vH0MwH4A13ZbVZqNHB9/u2n8VTkXBhmEAEDyPpIvA/g1AH/SdT01/CaAf+u6iAVwLYCX3fcH0OM/1CEjuQrgRgDf6LiUqUguk9wH4BUAXzKzyvX2NghJfpnks5GvOwDAzO41s50AHgZwd7fV5usdH3MvgFMY1dypKvX2XOwh4V6PDIaI5FYAjwL4yIaRWO+Y2TtmdgNGI66bSFa+BNHbRRfM7JaKh/4jgH8F8PEGy8nK1UvyNwB8EMDPWA8uzNb479tXBwDsdN/vAPDdjmpZSONrbY8CeNjMPtN1PVWZ2RGSTwC4DUClyanedoTTkLzefXs7gBe6qqUKkrcB+CiA283sRNf1LIinAFxP8t0kLwRwJ4DHOq5pYYwnHx4EsN/MPtV1PTkkV8LdGCQvBnALauTCUGeNHwXwAxjNbL4E4MNm9p1uq0oj+SKAzQAOjXc92fNZ7l8E8DcAVgAcAbDPzH6u06IiSH4AwP0AlgE8ZGb3dVtRGsl/AvDTGK2O8n8APm5mD3Za1BQkfxLAfwD4b4z+zgDgj8zs8e6qSiP5XgB/h9G/hSUA/2xmf1r59UMMQhGReRrk0FhEZJ4UhCJSPAWhiBRPQSgixVMQikjxFIQiUjwFoYgUT0Eog0byl0j+ddd1yLApCGXobgSwt+siZNh6u+iCyDQkvx/AAwB+DMAhktvN7P5uq5Kh0iN2MjgkN2O0Nt6HAHwOwE8AeB7ANWZ2ssvaZJg0NJYhuhXA0xgtu3XUzP4XwEmMHrgXqU1BKEN0A0arorwPwDMkrwZwzMxe77QqGSxdI5QhOgrgvRit9v0MgE9gdL1QZCa6RiiDQ/IKAJ8F8CMAXsXog5v+uA8rf8swKQhlsEg+DeBmMzuUPVhkCl0jlEEazxxvUwjKPKgjFJHiqSMUkeIpCEWkeApCESmeglBEiqcgFJHiKQhFpHgKQhEpnoJQRIr3/y2kWRbMLYHXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5, 2))\n",
    "plt.pcolor(phi, theta, np.abs(state1)**2,cmap=\"binary\",shading='auto')\n",
    "plt.xlabel(r'$\\phi$')\n",
    "plt.ylabel(r'$\\theta$')"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}